The coupled strength and toughness of interconnected and interpenetrating multi-material gyroids

The growth of additive manufacturing technologies has spurred interest in examining multi-material micro-architected materials for filling the so-called white spaces in the Ashby strength versus toughness plots. We investigate this problem using interconnected and interpenetrating double gyroids comprising ductile and brittle phases as an exemplar. Both strength and toughness at the initiation of crack growth are shown to vary non-monotonically with the volume fraction of the two phases and multi-material double gyroids significantly outperform their single material counterparts. However, we establish that at a given relative density, the strength and toughness cannot be simultaneously enhanced for architecture designs, which include varying gyroid orientations, phase volume fractions, and the unit cell length scales of the two phases. Intriguingly, even crack flank bridging by the ductile phase during crack growth is insufficient to overcome this inherent property of the interpenetrating gyroids. Our conclusion is that multi-material interpenetrating micro-architected solids are unlikely to outperform single material non-interpenetrating lattices from a strength–toughness perspective but rather become optimal when multi-functionality is required. The integration of materials and architectural features at multiple scales into structural mechanics gave us structural designs such as the Eiffel Tower. The explosion of additive manufacturing methods has opened new avenues for the invention of multi-material micro-architected materials that simultaneously possess high strength and toughness at a low density, and thereby can fill the so-called “white spaces” in the Ashby strength–toughness space. The idea is to construct three-dimensional materials with a network of crack arrestors like in rip-stop nylon and break the link between toughness and strength. We use interconnected and interpenetrating double gyroids comprising ductile and brittle phases as an exemplar to investigate the opportunities of such designs. Intriguingly, from a perspective based solely on strength and toughness, we show that multi-material micro-architectures cannot outperform their single material counterparts at a given relative density. In fact, in most designs the coupling between the two phases is non-synergistic. Rather, we argue that multi-material designs such as those used in rip-stop nylon are driven by multi-functional considerations beyond mechanical properties.


Introduction
Lattice materials are artificially structured materials that derive their properties primarily from their engineered topology. These materials have sparked great interest in the additive manufacturing (AM) community because AM is often the only feasible production route for these structures. In turn, lattice materials greatly expand the range of effective properties that are achievable. [1][2][3][4][5][6][7] While most three-dimensional (3D) printers can only print a single material there is a growing interest in multi-material lattices that is going hand-in-hand with the expanding AM capability to manufacture such materials. Advances in lattice materials have traditionally been driven by new unit cell topologies [6][7][8][9][10][11][12] or by manufacturing advances that leverage unique materials and nanoscale effects to maximise performance. 6,10-13

Impact statement
The integration of materials and architectural features at multiple scales into structural mechanics gave us structural designs such as the Eiffel Tower. The explosion of additive manufacturing methods has opened new avenues for the invention of multi-material micro-architected materials that simultaneously possess high strength and toughness at a low density, and thereby can fill the so-called "white spaces" in the Ashby strength-toughness space. The idea is to construct three-dimensional materials with a network of crack arrestors like in rip-stop nylon and break the link between toughness and strength. We use interconnected and interpenetrating double gyroids comprising ductile and brittle phases as an exemplar to investigate the opportunities of such designs. Intriguingly, from a perspective based solely on strength and toughness, we show that multi-material micro-architectures cannot outperform their single material counterparts at a given relative density. In fact, in most designs the coupling between the two phases is non-synergistic. Rather, we argue that multi-material designs such as those used in rip-stop nylon are driven by multi-functional considerations beyond mechanical properties.

ThE cOUpLEd sTrEngTh And TOUghnEss Of inTErcOnnEcTEd And inTErpEnETrATing MULTi-MATEriAL gYrOids
While novel unit cells remain undiscovered, incorporating multiple materials to access effective properties that fall within the "white-space" of Ashby material maps 14 is of growing interest. The idea of using disparate materials in composites to access toughening mechanisms is well-established. [15][16][17] It is even used in lattice-type materials such as rip-stop nylon where crack arresters in the form of a small volume fraction of ductile fibers are added to blunt an advancing crack. The goal of this work is to explore the use of multi-material lattices to achieve combinations of strength and toughness that do not exist in current engineering materials by attempting to break the inverse relation between strength and toughness 18 that exists across a wide range of engineering materials.
A large number of studies have been reported in the fracture of single-phase lattices. [19][20][21][22][23][24][25][26] Most of these are for twodimensional (2D) [19][20][21][22][23] lattices with very few investigations on 3D 24,26 lattices. The key finding of all these studies is that the mode-I fracture toughness K IC of a lattice with relative density ρ scales as K IC ∝ σ f ρ n √ ℓ , where σ f is the strength of the parent solid, ℓ is the strut length and n is a topology-dependent constant. These scaling laws have been experimentally verified for 2D 20 and 3D 24 lattices. The scaling with strut length implies that macroscale lattices are tougher compared to their nano-/microlattice counterparts, and in fact that the toughness of the lattice can far exceed its parent solid material toughness. 20 This feature has often led researchers to refer to these lattice materials as structural metamaterials. In fact, material property bounds do not exclude the possibility to invent a structural metamaterial with the strength and toughness of steel at the density of water. 14 This is a holy grail in structural metamaterials that has remained elusive primarily due to the inverse relationship between strength and toughness that is ubiquitous in most classes of materials, but not due to a theoretical requirement. The idea of using multi-material lattices [27][28][29] to obtain synergistic mechanical properties is gaining traction with interpenetrating lattices being the most important class of these multi-material systems.
We restrict attention to a class of lattices known as interpenetrating lattices (IPLs) 28 where two lattices made from different materials interweave through a given volume. If these lattices are disconnected, then there is no load transfer between the lattices and synergies in strength and especially toughness are unlikely to be achieved. Rather, we shall consider a class of 3D IPLs which are connected at a subset of their nodeswhile this restricts available topologies it is hoped it opens up synergistic interactions where properties of the interpenetrating lattice exceed the sum of the individual responses. The gyroid topology is a promising candidate whose mechanics has been extensively studied. 6,[30][31][32] Here, we explore the potential of chiral double gyroid lattices as IPLs to attain unique combinations of strength and toughness (Figure 1a). Using two lattices of opposite chirality permits the construction of an interpenetrating topology with the two opposite chirality gyroids connected to each other at a subset of their nodes (Supplementary Material Figure S1)-this is referred to as the interconnected double gyroid. The basic hypothesis being investigated is as follows (Figure 1a). Consider a single gyroid made from a ductile but low-strength material. This lattice is then expected to have a high toughness but low strength-both properties inherited from its parent material. Next, consider a similar gyroid but made from a brittle but high-strength parent material which then correspondingly has a low toughness but high strength. Our aim is to explore the potential of interconnected double gyroids comprising brittle (but high-strength) and ductile (but low-strength) phases to achieve properties that outperform a simple rule of mixtures ( Figure 1a) for fixed total volume fraction of the solid phases. In this pursuit, we also explore the role of different toughening mechanisms in such IPLs (Figure 1b). Although, a solely brittle phase lattice experiences catastrophic failure upon the initiation of crack growth, intrinsic toughening mechanisms are operative ahead of the crack tip in the ductile phase through the formation of a plastic zone. With the increase in crack growth, the plastically strained material unloads in the crack wake resulting in hysteresis and dissipation. At the same time, IPLs also offer extrinsic toughening in the crack wake through bridging by the unbroken struts of the ductile phase of the crack growing in the brittle phase. We hypothesize that both these mechanisms in IPLs can potentially result in a strong R-curve response, which might delay the onset of catastrophic fracture.
Let ρ denote the relative density (volume of solid material to volume of the effective smeared-out continuum) of the double gyroid comprising interconnected and interpenetrating single gyroids made from a brittle and ductile material with relative densities ρ B and ρ D , respectively (Figure 1c), so that ρ = ρ B + ρ D and ρ ≡ ρ B /ρ is the volume fraction of the brittle phase. We traverse the map in Figure 1a from left to right by varying ρ from 0 to 1, while keeping ρ fixed. In addition to the strategy of varying the relative proportions of the ductile and brittle phases, two additional schemes are investigated to enhance synergistic effects within the IPL. At the mm and larger length scales, the strut length ℓ and/ or unit cell length L = 2 √ 2ℓ ( Figure 1c) does not affect its strength, but its toughness (i.e., the critical strain energy release rate required to initiate a pre-existing crack, J IC ) is known to scale linearly with L. 14 We thus investigate a matrix of combinations ( Figure 1d) where we vary ρ and l ≡ ℓ B /ℓ D , where ℓ B and ℓ D are the length of the struts of the brittle and ductile single gyroid lattices, respectively. With this scheme of interpenetration, new sets of nodes with a connectivity of 5 or 6 are formed (they comprise ∼14% of the total nodes) by intersection of struts from the two opposite chirality gyroids ( Figure S1 and "Methods" section). Depending on the choice of ρ the struts of the two opposite chirality gyroids can have different diameters resulting in complex node geometries (Figure S1). We limit our attention to IPLs where each phase has ρ ≤ 0.125. In this range, the node volumes are negligible compared to the strut volumes and it thus suffices to model the struts as shear flexible Timoshenko beams in the finite element (FE) calculations ("Methods" section). In addition, we also exploit the cubic symmetry of the gyroid topology to determine the effect of the orientation of the lattice with respect to the crack (Figure 1e

Results
We present results for interconnected double gyroid lattices with parent material properties representative of steel. Thus, both the brittle and ductile phases have a Young's modulus GPa , respectively, and corresponding Poisson's ratios ν B = ν D = 0.33 . The brittle phase is modeled as a high-strength steel that is elastic-brittle with a tensile strength B = 2.5GPa such that it fails without undergoing plasticity. On the other hand, the ductile phase is modeled as an elastic perfectly plastic solid with a yield strength D = 1GPa and a ductility (i.e., the plastic strain at onset of stress softening) ε D . Results are discussed for the choices of ε D = 0.1 and 0.4 . Details of the numerical finite strain FE calculations, including boundary conditions and material failure models are provided in the "Methods" section. Unless otherwise specified, we restrict attention to the ρ = 0.1 lattice with l = 1 and discuss the effect of l in the final section of the manuscript. We emphasize that the focus of our study is on strength and toughness where the elastic properties of the parent phases have a minor influence. Thus, in the bulk of the manuscript, the analysis is restricted to the case where E B = E D and ν B = ν D . Further numerical calculations shown in Figure S2 confirm that varying the ratio R E ≡ E B /E D has little effect on the strength and toughnesses of the double gyroids.

Effective elastic-plastic properties of the interconnected double gyroid
We consider configurations shown in Figure 2a where the orientation and ratio ρ of the brittle to ductile phases are varied. Predictions of the normalized uniaxial compressive stress σ yy ≡ σ yy /� D versus strain ε yy are included in Figure 2b-c for the 0 o and 45 o orientation (with a failure strain ε D = 0.1 ) where the (x, y, z) co-ordinate system is defined in Figure 2a and σ yy is the stress conjugated to ε yy . The single gyroid has a nodal connectivity of 3, while the interconnected double gyroid topology has a mixed connectivity with some nodes retaining the connectivity of 3 of the parent single gyroid while it increases to 6 for ~ 14% of the nodes that are common between the two interpenetrating gyroids ("Methods" section and Figure S1). However, this increase in nodal connectivity is insufficient to change the deformation mechanism and the interconnected double gyroid remains a bending-governed lattice much like the single or disconnected double gyroid. 30 The Young's modulus of the interconnected double gyroid E DG thus scales with ρ 2 and the predictions of the normalized modulus E DG ≡ E DG / E B ρ 2 are included in Figure 2d as a function of ρ for the 0 o and 45 o orientations with E DG denoting the effective Young's modulus in the y-direction. The 45 o orientation is stiffer as it has more struts directly aligned along the loading y-direction. Also, since the modulus of the parent material is the same for both phases, E DG is equal for ρ = 0 and 1. More importantly, Figure 2d shows that the modulus is a minimum at ρ = 0.5 . To understand this, we observe that symmetry considerations dictate that under compressive loading in the y-direction, the interconnected opposite chirality single gyroids deform independently. Thus, the modulus of the double gyroid is the sum of the moduli of the individual single gyroids: where the scaling constants are obtained from the FE calculations as α ≈ 0.340 and 0.545 for the 0 o and 45 o orientations, respectively. The predicted dependence (1) is included in Figure 2d ( R E = 1 ) and not only accurately captures the FE predictions but also gives a clear insight into functional dependence of E DG on ρ . As ρ is increased from 0 (i.e., material is removed from a single to create the second interpenetrating phase), the bending stiffness is reduced. The non-linear scaling of the bending stiffness with the strut diameters gives rise to the ρ 2 scaling which in turn results in E DG being a minimum at ρ = 0.5 . A similar analysis for R E = 2 is presented in Figure S2. We note in passing that if the IPL is derived from stretching-dominated lattices (where the modulus scales linearly with relative density) with identical parent moduli (as in the case under consideration here), E DG is then independent of ρ.
Predictions of the full compressive stress versus strain curves (Figure 2b-c) show that while the brittle ρ = 1 lattice has an elastic-brittle response, varying levels of ductility are observed at lower ρ . The ρ = 0 case, which is a ductile single gyroid, displays a nearly perfectly plastic response until the onset of softening with the ductility of the 0 o orientation greater than the 45 o orientation. The double gyroids at intermediate values of ρ have load drops corresponding to failure of struts in the brittle gyroid with a complete loss of load carrying capacity occurring only when the ductile lattice ultimately fails. Note that this ultimate failure strain of the double gyroid lattice should not be confused with ε D , which is the tensile plastic failure strain (i.e., ductility) of the parent ductile phase strut. We define DG as the peak compressive strength of the double gyroid, and the bending-dominated response of the gyroid implies that DG ∝ ρ 1.5 . 30 Predictions of the dependence of the normalized strength � DG ≡ � DG / � D ρ 1.5 on ρ are included in Figure 2e Figure 2e and accurately capture the FE simulation data. In this scaling relationship, the constants β B and β D are obtained by calibrations for the single gyroid cases with extreme ρ values of 0 and 1, we have confirmed that (2) captures the strength accurately for double gyroids of relative densities in the range 0.05 ≤ ρ ≤ 0.125.

Toughness at the initiation of crack growth
We adopt a thin strip geometry to analyze fracture of the double gyroid lattice. This geometry comprises an infinitely long strip of height H with a semi-infinite crack and has been adopted extensively to investigate non-linear fracture. 33 The method allows the easy imposition of an J-integral that is independent of the crack length while not requiring smallscale yielding conditions to be satisfied. This method thereby significantly reduces the size of the computational domain ("Methods" section). Briefly, loading is imposed by applying a displacement rate u ∞ on the top surface of the strip as shown in Figure 3a and a simple energy balance provides a direct relation between the imposed u ∞ and J . A zoom-in of the crack tip extending within the 0 o and 45 o orientation lattices is shown in Figure 3b. Upon loading, crack growth occurs by the evolution of damage resulting in failure of the struts of the ductile and brittle lattices. The distributions of the von Mises stress σ e (normalized as σ e ≡ σ e /� D ) at the first instant of strut failure for the 0 o orientation double gyroid with ρ = 0.75 and ε D = 0.4 are included in Figure 3c. The associated plastic zone (the region with a von Mises effective plastic strain, ε p > 0 ) in the ductile phase is also marked in Figure 3c. These stress and strain distributions are suggestive of a mode-I crack tip field, and it is instructive to now define a toughness. However, the definition of crack extension and resulting toughness is not unambiguous in such IPLs, as while one lattice might fail the other could remain intact. Nevertheless, the definition of a crack flank is clear in that it is a traction-free surface, and this fracture is defined when new crack surface is formed with no bridging struts. In the brittle/ductile IPLs considered here, failure in the brittle lattice outruns failure of the ductile lattice such that the unbroken ductile struts form a bridging zone (see Figure 3d). We define the initiation of crack growth at the instant of the failure of the first bridging strut in the ductile lattice, and the value of J at this instant is labeled as J IC . Calculations using the more commonly used boundary layer method 21 were also performed as a check and gave identical results.
Predictions of the variation of the normalized toughness  (1) and (2) for modulus and strength, respectively. An extrapolation of numerical studies on two-dimensional lattices 21 would suggest that a scaling law of the form ρ is a non-dimensional function of R σ ,ρ and n an exponent that is only topology dependent. The numerical results in Figure 3e-f suggest that a scaling law of this form does not describe the non-linear fracture of the three-dimensional gyroid lattices-we attribute this to a combination of significant finite deformation effects at the crack tip even at the onset of fracture and bridging effects that are at play for the interpenetrating multi-material lattices. In fact, the scaling of toughness with ρ is highly non-linear as illustrated in Fig Figure 3g-h, and this is especially so for the ε D = 0.1 case. To understand this, recall that toughness scales approximately with the energy absorption per unit volume of the parent materials of the double gyroid. The ductile lattice has a high ductility but low strength, while the reverse is true for the brittle lattice. For a low ε D parent ductile solid, the ductility cannot always compensate for its lower strength. For example, increasing ρ from 0 to 0.25 for the ε D = 0.1 case in Figure 3h results in a decrease in toughness as replacing the ductile material by the higher strength brittle material cannot compensate for loss of energy absorption of the ductile lattice. However, in the range 0.5 ≤ρ ≤ 1 , the toughness increases with the increase in ρ as the additional load carrying capacity provided by the higher strength brittle lattice relieves stresses in the ductile phase and thereby delays ThE cOUpLEd sTrEngTh And TOUghnEss Of inTErcOnnEcTEd And inTErpEnETrATing MULTi-MATEriAL gYrOids fracture. This complex coupling is not amenable to simple analytical scaling laws necessitating detailed numerical calculations as presented here.

Coupling of compressive strength and initiation toughness of the multi-material lattices
Is there a synergistic coupling between the brittle and ductile phases of the interconnected and interpenetrating gyroids that offers the possibility of enhancing toughness without a loss of strength, or equally the ability to enhance strength without a loss of toughness? Predictions from   Figure 4b and qualitatively all conclusions remain unchanged. A question that arises at this point is whether single material interconnected and interpenetrating double gyroids are superior to the multi-material systems previously considered. We include in Figure 4a, b corresponding predictions for the ρ = 0.1 brittle interconnected and interpenetrating double gyroid (i.e., interpenetrating double gyroids each of relative density 0.05 but both made from the brittle parent material) and the corresponding ductile single material double gyroid for both the ε D = 0.1 and 0.4 parent materials. Clearly, higher toughness and/or strength is accessible via the multi-material route (i.e., interpenetrating multi-material lattices outperform single material counterparts). This is the key benefit of multi-material lattices but as previously clarified, the multi-material lattices cannot synergistically enhance strength and toughness together.

Does bridging during crack growth allow the development of synergistic coupling between the phases?
Crack bridging is known to be a potent toughening mechanism in a vast class of composite materials 15 and interconnected and interpenetrating ductile and brittle lattices allow for the possibility to exploit such crack bridging effects in 3D architected lattice materials. Although such toughening mechanisms have been well explored in a wide range of fully dense multi-phase systems, including composites with soft and stiff phases, [34][35][36] fiber and particulate composites 15,37 and biological systems like bone, 18,38 they remain unexplored for 3D truss-based architected metamaterials. We now exploit our strip geometry numerical setup (Figure 3a) to calculate resistance or J-R curves (i.e., evolution of the imposed J with crack extension a ). Recall that we define the crack flanks as being traction-free and thus the new crack tip is located at position where the lattice is contiguous ahead of it: the crack extension a defined as the difference between the current and original crack locations as measured in the x-direction. Predictions of the distribution of the normalized von Mises stress σ e and corresponding plastic zone (defined as the region with ε p > 0 ) at the crack tip along with the plastic wake in the ductile phase are included in Figure 5a for the ρ = 0.75 lattice with ε D = 0.4 in the 0 o orientation. The results are shown for a crack extension �a/L = 100 (i.e., after 100 unit cells of the ductile lattice are broken). As anticipated, the stress fields are similar to Figure 3c but translated to the new tip, but the crack extension leaves a plastic wake with plasticity in the zone spanning ~6 unit cells on each side perpendicular to the crack flank. Moreover, similar to Figure 3d, the crack in the brittle lattice outruns the ductile lattice with the ductile struts bridging the crack flanks in the brittle lattice. Thus, two toughening mechanisms 18 are seen to operate (also cf. Figure 1b): (1) toughening due to the plasticity in the wake of the crack-a mechanism seen in continuum plastic solids and (2) crack tip bridging as seen in fiber and other composite materials. Previous work 21,25 on bending-dominated lattices suggests that the contribution to toughening from the former mechanism will be low due to their small plastic zone sizes. Predictions of the normalized J-R curves, viz., J ≡ J /(� D L) versus �a ≡ �a/L are included in Figure 5b, c for the 0 o orientation case for ε D = 0.1 and 0.4, respectively, across the range of ρ . After the initiation of crack growth, J increases with the increase in a in all cases although the rise is steepest for the low ρ cases with ε D = 0.4 where the ductile phase dominates.
To illustrate the development of bridging by the ductile phase with crack propagation, we include in Figure 5e Figure 5d. The trends are qualitatively similar to Figure 4a suggesting that the basic conclusion that multi-material interpenetrating double gyroid designs do not allow for a synergistic improvement of the strength-toughness coupling remains unchanged.

Designs with length scale mismatches between the two phases
Toughness scales with the strut length of the lattices, and a pertinent question is whether including a contrast in length scales between the interpenetrating lattices induces the synergy between strength and toughness that we are seeking. The calculations previously discussed were restricted to l ≡ ℓ B /ℓ D = 1 where the ductile and brittle lattices of opposite chirality were otherwise topologically identical. We now consider two additional cases of l = 0.5 and 2.0 such that the period of the brittle lattice is half and twice that of the ductile lattice, respectively. The corresponding unit cells of such lattices for a matrix of ρ and l values are sketched in Figure 6a; see the "Methods" section for a detailed description of the procedure to generate these lattices and consequence on nodal connectivity. Three-dimensional views of the crack extending in such lattices with ρ = 0.75 are shown in Figure 6b for l = 0.5 and 2.0. The sketches in Figure 6a clearly illustrate that the unit cell size L of the interpenetrating lattice is set by the lattice with the larger length scale, so that . This is more clearly seen in full 3D images of additively manufactured interpenetrating double gyroid lattices shown in Figure 6c-g. For visual clarity, the ductile and brittle phases are dyed with black and yellow tracer particles, respectively. These lattices were printed using microstereolithography that has recently been adapted to allow simultaneous printing of two different materials, 39 and extended to a large area projection system by continuously moving projection so curing occurs in sub-sections. 40 In Figure 6c, we include an optical image as well as an X-ray computed tomographic (XCT) scan of the l = 0.5 lattice with a zoom-in showing nodes with a connectivity of 5. Optical images of four lattices for a selection of ρ and l values are included in Figure 6d ThE cOUpLEd sTrEngTh And TOUghnEss Of inTErcOnnEcTEd And inTErpEnETrATing MULTi-MATEriAL gYrOids included in Figure 6h-i for the choices ε D = 0.1 and 0.4, respectively, and length scale contrasts of l = 0.5 , 1.0 and 2.0. The overall qualitative behavior remains unchanged, viz., allowing different length scales for the brittle and ductile lattices does not allow the creation of materials that can simultaneously enhance strength and toughness. This conclusion carries forward to the crack growth resistance ( Figure  S3) as well as the 45 o orientation (Figure S4). We emphasize that the non-dimensional toughnesses in Figure 6h, i need to be interpreted with care. For example, if the l = 2.0 double gyroid is generated by keeping the ductile lattice identical to the l = 1 case but doubling the periodic cell size (i.e., L ) of the brittle lattice, then the l = 2 lattice has an actual toughness J IC that is twice that of an l = 2 lattice generated by keeping the brittle lattice size equal to the l = 1 case and halving the size of the ductile lattice. However, the normalized J IC results (Figure 6h-i) remain unaffected in either case.

Outlook
The geometry of interpenetrating networks or lattices is well established within the framework of Euclidian geometry, viz., dual polyhedral, as reciprocal pairs of polyhedral where the vertices of one polyhedron correspond to the faces of the other, and both polyhedral share the same symmetries. 41 Crystallographers using these geometries refer to the resulting constructions as Wigner-Seitz cells. 42 This idea has also recently been adopted by the metamaterials community with nomenclature adopted from crystallography as well as classical geometry. 43 The concept of dual or interpenetrating lattices is thus the latest adaptation of classical ideas in crystallography to invent new metamaterial topologies although previous work 27, 28 in such interpenetrating lattices has focussed on disconnected phases. While the concept is simple, the explosion of additive manufacturing technologies has opened routes to manufacture such materials and capitalise on the design freedom to achieve unique properties. Additive manufacturing techniques such as projection microstereolithography 44 and two-photon lithography 45 have led to realization of polymeric, 46 metallic, 47 and ceramic 48 metamaterials with micron and sub-micron unit cell sizes. We have included in Figure 1f-g as well as in Figure 6c-g optical and XCT images of multi-material double gyroids printed via microstereolithography with the two interpenetrating but interconnected phases comprising the brittle TMPTA * lattice and ductile lattice made from a mixture of BPAEDA † and TMPTA (9:1 v/v). These 3D printed multimaterial lattices serve as a proof of concept for the ability of the state-of-the-art manufacturing methods to fabricate such complex multi-material systems, wherein the interconnected nodes between struts of different materials could have widely varying diameters ( Figure S1 and Table SI). The scale-up of this method to manufacture metamaterials with 10s of millions of unit cells 40 implies that such interpenetrating lattices are now viable candidates as large-scale engineering materials.
The key advantage of the double gyroid lattice is that the constituent interpenetrating networks are single gyroids of opposite chirality. Moreover, it is straightforward to induce connectivity between the two networks by translation of one with respect to the other ("Methods" section and Figure S1). The interconnected and interpenetrating gyroids of opposite chirality considered here are a bending-dominated topology. 30 The increase in nodal connectivity by inducing interconnectivity within them is insufficient to switch their behavior to stretching-dominated. This is because the fraction of nodes with a connectivity greater than 3 is still low (~ 14% for l = 1 and ~ 11% for l = 0.5 and 2) in the range of lattice architectures considered ("Methods" section). A direct consequence of this low average connectivity and consequently the bendingdominated behavior is that strength and toughness to vary in a non-monotonic manner with the volume fraction of the two phases (Figures 2 and 3). This low average connectivity of the interconnected and interpenetrating gyroids is a critical limitation as further increase in the average connectivity is not possible as that results in an overlap of struts of the two phases making such IPLs infeasible to manufacture.
It is conceivable that other topologies might exist where the interconnectivity is high which in turn might result in a synergy between strength and toughness being achieved. One possibility is the construction of stretching-dominated interpenetrating and interconnected architectures such as the socalled compound trusses. 43 These will undoubtedly enhance the strength/toughness-to-weight ratios of the interpenetrating lattices. However, these lattices are also unlikely to succeed in achieving the goal of inventing topologies that simultaneously enhance strength and toughness at a fixed relative density. The reason for this in the context of lattice materials can be understood as follows. Enhancing the ratio of the strong (but brittle) to weaker (but ductile) phase will increase the strength of the lattice. This higher strength comes at the cost of increased strut stresses, including higher crack tip stresses that increase the propensity for fracture and thus reduce toughness. This basic mechanics is inherent to lattice-based metamaterials including interpenetrating lattices.
Notwithstanding, different variants of IPLs are extensively used. For example, rip-stop nylon has crack arresters in the form of a small proportion of ductile threads interwoven into the fabric that comprises otherwise brittle threads. These ductile threads blunt an advancing crack to enhance toughness. However, this has been achieved by increasing the weight of the fabric by including the ductile threads. Our investigations suggest that at a constant relative density (i.e., removing a certain proportion of the brittle threads to add ductile threads) this increased toughness will come at the price of reduced strength. This focus on mechanical properties, however, neglects that non-mechanical design considerations are often critical in many material selections. IPLs will find applications in the multi-functional property space where properties in addition to strength * Trimethylolpropane triacrylate. † Bisphenol A ethoxylate diacrylate. and toughness become important. For example, heat shields for hypersonic aircraft not only require excellent mechanical properties but also zero or low thermal expansion coefficients 49 to minimise the generation of thermal stresses. Similarly, cathodes in Li-ion batteries ideally should not swell upon lithiation and must possess a high compliance to reduce their propensity to cracking. 50 IPLs are already known to provide some of these multi-functional properties, and it is possible to imagine numerous unusual mechanical, thermal, chemical or electrical functions that such multi-material lattices can fulfil.

Construction of double gyroids
The gyroids are lattices with a cubic symmetry. Let e i where (i = 1, . . . , 3) denote the orthogonal axes aligned with the cubic directions of the gyroid. Then the topology of a single gyroid with period L is well-described 51,52 by the function with (x, y, z) denoting co-ordinates along (e 1 , e 2 , e 3 ) , respectively, and t 0 a scaling parameter that sets the relative density of the gyroid. The double gyroid comprising two interpenetrating but disconnected gyroids of opposite chiralities is then constructed by infilling the spaces bound by F − t 0 ≥ 0 and F + t 0 ≤ 0 to obtain the unit cell shown in Figure S5a. For the range of relative densities we considered here, it suffices to approximate the full 3D geometry by a network of 3D beams with a circular cross section of radius r and strut length ℓ = L/ 2 √ 2 , as shown along different projected views in Figure S5b. The nodes (i.e., connection points of the beams) of this network are located at the geometric centre of the nodes of the full 3D topology. The relative density of the ductile gyroid is then with ρ B = ρ Dρ / 1 −ρ and r 2 B = r 2 Dl 2ρ / 1 −ρ and ℓ B =lℓ D . These expressions hold whether the interpenetrating gyroids are connected or disconnected as the proportion of the overlapping nodal volume for different ρ cases considered in this work ( 0.05 ≤ ρ ≤ 0.125 ) is low.
To construct the interconnected interpenetrating double gyroids, first consider the l = 1 case such that ℓ B = ℓ D = ℓ . The interpenetrating gyroids are connected by translating one of the gyroids by L/4 ê 1 +ê 2 +ê 3 where ê i are unit vectors along e i . Then a fraction of the nodes of the two gyroids overlap as shown in Figure S1a and the two gyroids are fixed to each other at these nodes. This interconnected double gyroid retains cubic symmetry but has a mixed connectivity: a unit cell comprises 14 nodes such that two of these nodes have a connectivity of six and involve three struts each from the two opposite chirality gyroids, while all the remaining nodes have , a connectivity of three and involve struts belonging to a single gyroid; see Figure S1a. The generation of the l � = 1 interconnected double gyroids is more complex as we now illustrate for the l = 2 case. We first construct the surface of the large cell gyroid (the brittle single gyroid in this case) by using the equation (3) and draw the network of beams that form the struts of this gyroid. Next, the ductile single gyroid is drawn with L = 2 √ 2ℓ B /l in (3) and skeletonised to form the network of struts. The smaller unit cell lattice (ductile gyroid in this case) is then translated by L/16 5ê 1 −ê 2 + 11ê 3 to form the unit cell of the interconnected double gyroid (see Figure S1b). In this unit cell, (1) two nodes are generated by overlapping of the nodes of the original ductile and brittle gyroids so that these nodes now have a connectivity of six with three struts from each of the single gyroids; (2) six nodes formed as a result of the intersection of the struts of the ductile gyroid with the midpoint of struts of the brittle gyroid to form new nodes: these nodes have a connectivity of five and involve three struts from the ductile gyroid and two struts from the bisected strut of the brittle gyroid; and (3) the remaining nodes remain unchanged from their parent gyroid connectivity. Construction of the l = 0.5 gyroid follows an identical procedure except that now the brittle gyroid is translated by L/16 5ê 1 −ê 2 + 11ê 3 to form an interconnected and interpenetrating double gyroid where the connectivity is as described above with the ductile and brittle gyroids swapped (Figure S1c). Table SI reports the strut radii ratios of the ductile and brittle phases for all the interpenetrating and interconnected gyroids considered in this study. The previously discussed procedure generates the gyroids in the 0 o orientation and the 45 o orientation is obtained by rotating these gyroids by 45 o about any of the e i -axes (the cubic symmetry means that all the e i -axes are topologically identical).

Finite element (FE) calculations and material model
All calculations were performed in the finite deformation setting using the general-purpose FE code ABAQUS. The struts were modeled by discretising each strut using eight 3D Timoshenko beam elements (B31 in the ABAQUS notation). Both the ductile and brittle phases were assumed to be made from steel, which is modeled as an isotropic J2 flow theory solid with a Young's modulus E B = E D = 200GPa (i.e., R E = 1 ) and Poisson's ratio ν B = ν D = 0.33 . The brittle phase was taken to be representative of a Maraging steel with a yield strength B = 2.5GPa with a softening response immediately after yield. The ductile phase represents a TRIP or low carbon steel with a yield strength D = 1GPa (i.e., R σ = 2.5 ) followed by a perfectly plastic response until the initiation of softening at a tensile plastic failure strain ε D . Beyond, the onset of softening a linear softening response was assumed for both the ductile and brittle phases until the load carrying capacity of the material vanishes. This fracture energy after the onset of softening for both phases was set to be G f = 25kJm −2 . Thus, the plastic strain at a complete loss of load carrying capacity is 2G f /(� D e) + ε D and 2G f /(� B e) for the ductile and ThE cOUpLEd sTrEngTh And TOUghnEss Of inTErcOnnEcTEd And inTErpEnETrATing MULTi-MATEriAL gYrOids brittle phases, respectively, where e is the element size used to discretise the beam elements. To reduce the effect of node designs on the predictions (the details of the node designs are not modeled in the beam description of the lattices which is a good approximation for the low relative density lattices considered here), we only include the softening response in elements at the centre of each strut. Thus, while the beams have homogenous properties, failure is restricted to occur only at the centre of each strut. The uniaxial compression responses were calculated by considering a single unit cell and imposing periodic boundary conditions where the strain ε yy (Figure 2a) was specified and the stress components conjugated to all the remaining strain components were set to zero.

Calculation of fracture toughness and R curves
The initiation toughness and R-curves were calculated using thin strip specimens in the x-y plane with the cubic e 3 -direction aligned with the z-axis (Figure 3a). Interpenetrating gyroid lattices with one unit cell thickness in the z-direction were sandwiched between the platens with the crack extending along the x-direction and the y-direction normal to the crack plane. The specimens of thickness H comprised 101 unit cells in the y-direction and 600 unit cells in the x-direction with the initial crack extending 100 unit cells along the x-direction. The left and right surfaces of the specimen were traction-free, and all degrees of freedom (DOFs) were constrained to zero on the bottom surface. A displacement rate u ∞ is imposed uniformly on the top surface in the y-direction with all other DOFs set to zero. Periodicity was imposed for all DOFs between the front and the back face nodes of the sandwiched gyroid. For an imposed u ∞ , the J-integral then follows directly as J = WH where W is the strain energy per unit volume in the double gyroid strip far upstream from the crack tip at the imposed displacement u ∞ . 33 Application of an increasing displacement u ∞ (or J ) results in failure of struts at the crack tip and the consequent propagation of the crack. The initiation toughness is set by the value of J when the first ductile lattice strut fails (for l = 0.5 and 1 ) or when the first brittle strut fails (for l = 2 ), while the R curves follow directly from the measurement of a as discussed in the main text. This method requires an independent estimate of W as a function of u ∞ /H . This calibration is obtained as follows. We take a strip specimen as above of height H but width w = 10L with no crack. All DOFs other than the displacement in the y-direction are constrained to zero on the left and right surface nodes of the specimen. By imposing an increasing displacement u ∞ , these calculations directly provide W (u ∞ /H ) as W ≡ 1/(wHL) F ∞ is the force conjugated to the imposed u ∞ .