Encoding kirigami bi-materials to morph on target in response to temperature

Shape morphing in response to an environmental stimulus, such as temperature, light, and chemical cues, is currently pursued in synthetic analogs for manifold applications in engineering, architecture, and beyond. Existing strategies mostly resort to active, namely smart or field responsive, materials, which undergo a change of their physical properties when subjected to an external stimulus. Their ability for shape morphing is intrinsic to the atomic/molecular structure as well as the mechanochemical interactions of their constituents. Programming shape changes with active materials require manipulation of their composition through chemical synthesis. Here, we demonstrate that a pair of off-the-shelf passive solids, such as wood and silicone rubber, can be topologically arranged in a kirigami bi-material to shape-morph on target in response to a temperature stimulus. A coherent framework is introduced to enable the optimal orchestration of bi-material units that can engage temperature to collectively deploy into a geometrically rich set of periodic and aperiodic shapes that can shape-match a predefined target. The results highlight reversible morphing by mechanics and geometry, thus contributing to relax the dependence of current strategies on material chemistry and fabrication.


S1. Testing apparatus
Supplementary Fig. S1 shows the schematic of the experimental set-up built to test the thermal deformation of the fabricated metamaterials. The heating chamber consisted of two 200-Watt strip heaters (McMaster-Carr, USA) placed underneath the testing plate, where the temperature was adjusted through a proportionintegration-differentiation (PID) controller (CN7800, OMEGA, USA & Canada). A data acquisition system (cDAQ-9174, National Instruments, USA) was used to collect the temperature values from the thermocouples placed at different locations in the heating chamber with the goal of assessing any instance of temperature heterogeneity throughout the heating chamber. A borosilicate glass cover was placed on the top of the testing system to provide a thermal insulation shield. Specimens were tested both in air and oil, the latter to provide a uniform source of heating and eliminate frictional dissipation between the specimen and the testing plate.
Digital Image Correlation was used to capture the full-field displacement and effective thermal strain of the specimen at increasing levels of temperature. Supplementary Fig. S2 illustrates a representative set of testing on a periodic metamaterial with compound units. Any rigid body movement of the specimen was prevented by anchoring one point (highlighted in red in Supplementary Fig. S2a) of the specimen to the testing plate. Two pairs of black speckles were applied on the rigid frame to trace their movement.
Temperature gradually increased from 20 to 120 °C. Sample deformation was first captured by a digital camera (EOS Rebel T6i, Canon USA), and then processed with a correlation algorithm (Vic-2D, Correlated Solutions Inc., USA) that provided the full-field displacement and strain data between pairs of black speckles in the deformed specimen. The effective thermal strains ( * ) was then obtained from the relative displacement normalized by the initial distance between pairs of black speckles. As an illustrative example, regions with speckles dispensed on the wooden frame (see insert on the top of (a)).

S2. Characterization of constituent properties: thermal expansion coefficient and Young's modulus
A thermomechanical analyzer (TMA Q400, TA Instruments Inc., USA) was used to measure the coefficient of thermal expansion (CTE) of the constituent solids, silicone rubber and hardwood. For the latter, the CTE was measured in both directions, parallel and perpendicular to the wood grain. Supplementary Fig. S3a reports the CTE for both materials with measured values ( = 215 × 10 −6 −1 , // = 5 × 10 −6 −1 , and ⊥ = 32 × 10 −6 −1 ) almost constant throughout the temperature range under investigation (20 to 120 ℃).
To assess the tensile elastic modulus of the constituents, we tested a set of laser-cut dumbbell-shaped specimens ( Supplementary Fig. S3b) with an Instron tensile tester (5982 Series Universal Testing Systems, Instron Inc., USA). For silicone rubber, 3 uniaxial tensile tests were performed under displacement control, from which the most representative nominal stress-strain curve was obtained. To capture the experimental measures, a first order Ogden model was adopted (Supplementary Fig. S3c) with strain energy function given by: ( 1 , 2 , 3 ) = ( 1 + 2 + 3 − 3) where and are the material constants chosen to fit the experimental data (in this case = 1.03 and = 6.20), and 1 , 2 , and 3 denotes the principal stretches. In addition, the elastic modulus of silicone rubber was extracted from the curve at 25% strain ( 25% = 3.7 MPa) and used in our analyses.
For the elastic response of wood, uniaxial tensile tests were performed on specimens with orientation either parallel or perpendicular to the wood grain direction ( Supplementary Fig. S3b). Ten dumbbell-shaped specimens were laser cut from a 1/8-inch-thick sheet of hardwood. Supplementary  (d) Stress-strain curves obtained from testing hardwood specimens oriented in the direction parallel and perpendicular to the wood grain.

S3-A. Theoretical model
This section examines the mechanics of the functional building block shown in Fig.1a subject to a uniform thermal field. We theoretically study the temperature-driven deformation for both the unidirectional and rotational floppy modes. Theoretical expressions for the elastic deflection of both U-BB and R-BB are obtained as a function of temperature, , under the following assumptions: • Both constituent solids have constant CTE within the temperature range here investigated.
• Hooke's law of elasticity holds for the high CTE core.
• The low CTE frame is assumed rigid, therefore the deformation induced by the thermal mismatch of the two materials is neglected.
• The length of the connection between the upper and lower portions of the high CTE core is considered significantly smaller than the BB length ( ≪ ).

i. Unidirectional Building Block (U-BB)
Supplementary Fig. S4a shows U-BB in its undeformed state: two distinct portions joined centrally via a connection of length along the horizontal axis of symmetry. Each is shaped with a low CTE (grey) and high CTE (blue) material, where and width ℎ describe the length and width of the core, which is cut along semi ellipses with major axis of length 2 and minor axis of length 2 . The thermal mismatch between the low CTE confinement (grey) and the high CTE expansion (blue) ( Supplementary Fig. S4b) governs the unidirectional response that is endowed by the BB topology. A temperature increase expands the core, which is constricted horizontally by the surrounding frame, and thus forced along its floppy direction to space the upper and lower portion of the BB by ∆ℎ, here assumed anchored at its center. Due to the negligible thermal expansion of the low CTE solid, an effective boundary condition is enforced to replace the frame action onto the core. For the sake of symmetry, a quarter of the core is examined with one end clamped and the other free to move transversely via a guided support ( Supplementary Fig. S4c). In addition, we simplify the semielliptical quarter of the groove by straightening its profile so as to form an inclined beam, which is in turn attached to the preceding horizontal part. The analysis of the building block subjected to uniform temperature is now reduced to the solution of a statically indeterminate problem of a Timoshenko beam-column (Supplementary Fig. S4d) with an equivalent CTE, effective = core − frame , compensating the neglect of the thermal expansion of the frame.
Supplementary Fig. S4: (a) Geometry of unidirectional building block, U-BB, with its dimensional parameters. (b) Deformed U-BB with expansion ∆ℎ depicted at a given temperature, . (c) Condensed model of U-BB, where the nominal geometry is reduced to that of an effective core with effective = core − frame subject to ∆ and pertinent boundary conditions. (d) Further model reduction to a statically determinate Timoshenko beam-column subjected to reaction forces, and , both dependent on T.
By releasing redundant constraints on the clamped end of the beam-column, the model in Supplementary and bending moment ( ( )) that both depend upon the level of temperature, as well as its deflection can be obtained through the solution of the following set of second-order differential equations: where is the transverse deflection along the y axis of the beam, and it represents the difference between the deformed and undeformed axis of the beam-column. and are the Young's modulus and the second moment of area of the beam cross-section. The following notation applies: which is substituted into Eqns.
(2) and yields: The general solutions of Eqns. (4) have the form: { ( ) = 1 sin( ) + 2 cos( ) + , 0 ≤ ≤ 1 ( ) = 3 sin( ) + 4 cos( ) + ( 1 − ) + , 1 < ≤ 2 (5) The constants of integration 1 and 2 can be determined from the boundary conditions. Since the beam is clamped at the left end, we impose: To determine the constants of integration 3 and 4 , continuity on deflection and slope is enforced at the junction between the two portions of the beam, thus yielding: Substituting into Eqns. (5) the values of the integration constants from Eqns. (6) and (7)  With the above equations of the deflection curve, we can determine the force and moment by imposing pertinent compatibility conditions. The first specifies that at the right end support of the beam: Since a compact closed-form for ( ) and ( ) cannot be directly retrieved from the above equations, we numerically solve Eqns. (13) and (10) and use the Newton-Raphson method to find their roots. The deflection curve ( ) of the beam at a given temperature can then be obtained by substituting the roots of ( ) and ( ) into Eqns. (8). The expansion of the building block ∆ℎ at a given temperature is thus determined by calculating the transverse displacement of the beam at point A, expressed by: ii. Deformation mode of U-BB The above theory can now be used to investigate the deformation mechanism of U-BB for an increasing level of temperature. To do so, we examine a representative beam-column with given parameters /ℎ =9, 2 / =0.1, and 2 /ℎ =0.2 and plot the reaction and bending moment at the beam end. Supplementary Fig.   S5a shows the results with the axial force and bending moment plotted as a function of the temperature range here investigated. To rule out size dependency in the plot, and are normalized with the cross-sectional dimensions ℎ/2 and , the latter representing the out-of-plane thickness of BB.
Two sequential regimes of deformation can be observed, each controlled by temperature. For low values of temperature, the axial force dominates the bending moment and increases linearly with temperature; here axial compression governs the U-BB response. With a further increase of temperature, the deformation mode switches. Bending increases steeply and non-linearly, and the axial force flattens at a plateau with a resulting bending domination at high temperature. This is also apparent in Supplementary Fig. S5b, where a transition from a stretching to a bending mode appears in an intermediate zone (grey) around T=60 °C.
Below this zone, i.e. in the low temperature regime, the change in slope of ∆ℎ/ℎ shows a moderate expansion induced by extension. In contrast above the transition zone, U-BB expansion is controlled by internal bending of the core with a rate of ∆ℎ/ℎ that increase swiftly and non-linearly with temperature.
Supplementary showing BB sensitivity to temperature variation in the low and high temperature regime, respectively below and above the transition zone.

iii.Rotational Building Block (R-BB)
R-BB differs from U-BB for a mere symmetry breaking as well as end closure of the core. They both seal R-BB with a rotational character. Supplementary Fig. S6a shows a shift in the position of the connection of the core at an offset from the left end of R-BB. The underlying mechanism of thermal deformation does not differ from that of U-BB, i.e. the trigger is the CTE mismatch of the constituents, but the topological difference between the two (U versus R) is responsible for each floppy mode ( Supplementary   Fig. S6b).
First, Eqn. (11), which for U-BB retains beam length compatibility upon heating, still holds for R-BB.
Hence, substituting the expression of ( ) − ( 0 ), thermal , and compression into Eqn. (11) results in: In addition, the following two apply to the right end support of the beam: The fourth condition stems from the deformation symmetry along the OA line that exists in the deformed R-BB ( Supplementary Fig. S6b). This translates into an equivalent condition of slope between point A (offset of e from O), and the symmetry axis (OA), here expressed as: Also in this case, we resort to the Newton-Raphson method to numerically solve the reaction force , , , from the system of Eqns. (18)-(21). After the deflection curve of the beam at a given temperature is obtained by substituting their roots into Eqns. (15), we can express the opening angle, , of R-BB at a given temperature through the slope of point A on the deflection curve as:

iv.Deformation mode of R-BB
Similar to the U-BB investigation on the mechanism of deformation, here we study the deformation behaviour of R-BB. As a representative model, a beam-column is examined with /ℎ =9, 2 / =0.1, 2 /ℎ

S3-B. Experimental validation of BB theoretical models
To verify the validity of the solutions obtained from the previous theoretical analyses, we fabricated a number of building blocks with a specific set of geometric parameters and tested them in a heating chamber (see S1 for details). A representative array of results from these experiments is shown in Supplementary   Fig. S8 for both U-BB and R-BB. The experimentally measured values of Δℎ/ℎ and of actual BB geometry confirms the mechanism of deformations predicted through the mechanics theory on simplified column-beam analogs. In Supplementary Fig. S8, the experimental results show a good agreement with the theoretical prediction, especially in high temperature region (max error 15%). The discrepancy can be attributed to frictional dissipation and the simplified geometry assumed for the BB. In the low temperature regime, the former plays a major role, since the testing plate upon which the BB rests during the experiment is not perfectly smooth and hence the frictional force can significantly dominate the low magnitude of the internal forces generated in BB (Supplementary Fig. S7a and b). In the high temperature regime, on the other hand, the simplification of the semielliptical groove of BB core to a straight inclined beam is the main source of the deviation.
Supplementary Fig. S8: Experimental validation of theoretical predictions for both BBs: (a) normalized expansion ∆ℎ/ℎ versus temperature for U-BB tested with two representative sets of parameters (blue and red); (b) opening angle as a function of temperature for R-BB tested with two representative sets of parameters. Horizontal and vertical bars indicate the standard deviation around the mean from a pool of measures taken for both temperature and deformation, the former measured at three distinct sites in the heating camber, the latter obtained from three repeated measurements.

S3-C. Computational analysis of BB response
The theory presented in S3-A and experimentally validated in S3-B applies to a BB with core of small size ligament (d<<l) and groove geometry abridged with a straight axis. Here we relax these assumptions and examine a BB with actual (un-simplified) core geometry made of two distinct passive materials: silicone rubber and wood (see S2 for material properties). Section S3-C-i explains the method used to generate the deformation property capturing the role of the ligament size (d/l) besides the core size (h/l).
Section S3-C-ii investigates the minor role of other geometric parameters defining U-BB and R-BB, in particular, the groove size, the core connection offset and the size of the end closure.
All the computational analyses here performed account for both material and geometric non-linearity of the constituents and are conducted by combining Python scripts and Abaqus (Dassault Systemes Simulia Corp, France). Eight-node, quadratic plane-stress elements (element type CPS8) are used to discretize the BB models and a mesh size equal to /4 is adopted after performing a convergence analysis on the mesh size. The mechanical behavior of silicone rubber core is captured by a first order Ogden model. The effective properties of the wooden frame are assumed as transversely isotropic, and the values are obtained from experimental measures on dog-bone specimen of the solid material (see S2 for material properties).
The cut motifs introduced into the computational models are represented with seam cracks along which duplicate element nodes overlap. A contact law with hard contact for the normal behavior, and frictionless for the tangential behavior, was assigned to the model.

i. Deformation-property profile
The deformation property profile in Fig.1a is generated through a response surface that fits a set of numerical predictions with model details given above. In general, the relationship between the response of a system (y) and a set of predictor variables ( 1 , 2 , … , ) can be mathematically expressed as a multiple linear regression model which can be written as where y is the true response variable and the parameters ( = 0,1, … , ), are the regression coefficients.
This model spans a k-dimensional space defined by the regressor { }. The parameter is the error of the regression model.
To provide a continuous approximation of the true response, the material and geometry spaces of the deformation-property profile of BB are generated from two sets of simulations. These are conducted in the admissible domain of two sets of variables: ( / , ∆ ) and ( /ℎ, / ), and for each of them a response surface is obtained. In general, the response of an N-order model as a function of two variables, 1 and 2 , is given by: where y is the true value of the response, in this case, ∆ℎ/ℎ for U-BB and for R-BB; 1 and 2 are / , ∆ for the material space, and /ℎ , / for the geometric space. If we let , = ( 2 +3 −2 ) 2 ⁄ , and To estimate the regression coefficients , we use the method of the least squares. Writing Eqn. (23) in matrix notation gives: where: where the index n represents the number of sampling points in the design of experiments.
Since the goal is to find the regression coefficient vector that minimizes the error vector , we can write the set of least squares function as: Eqn. (27) can be further developed into: with the least-squares estimators 0 , 1 , … , satisfying the condition: Eqn. (29) simplifies to the normal equations in matrix form: Solving the normal equations gives the least-squares estimator b of the regression coefficients : Hence, the fitted regression model is: The residuals are: The coefficient of determination is: From the above we can respectively write the sum of squares of the regression , the sum of squares of the residual , and the total sum of squares as: To reduce the error of the approximation, we aim at ensuring the coefficient of determination, 2 , above 0.99, a value that indicates here an acceptable level of accuracy. As per choice of the order of the response surface function (Eqn. (24)), we adopt = 4 after performing a systematic analysis on the role of N.
While Fig.1a shows the response surface of ∆ℎ/ℎ for U-BB, i.e. the property-deformation profile, here we report, as example, the expression of the surface response of for R-BB in its geometric space: where the regression coefficients , is estimated through Eqn. (31). A similar approximation can be expressed in the complementary material space for R-BB as a function of the pair of material properties ( / , ∆ ). Below are two plots constituting the property-deformation profile of R-BB response at = 120 °C (Supplementary Fig. S9a and b). Obtained through Eqn. (38), the spectra provide the range of rotational deformation that a BB can achieve through tuning the most influential geometric and material attributes that control R-BB response to temperature. The maps provide guidelines for attributes selection of BB, in particular, the geometric space on the right-hand side is the foundation of our morphing scheme.
Supplementary Fig. S9: Property-deformation profile of R-BB ii. Sensitivity to least influential geometric parameters

a. Role of groove size ( / and / ) in U-BB
While the surface response in Fig.1 shows the role of the most influential geometric parameters, /ℎ and ⁄ , on ∆ℎ/ℎ, Supplementary Fig. S10 reports complementary results capturing the influence of the semielliptical groove of the core: (2 / ) and (2 /ℎ). For a BB with /ℎ equal to 9, the curves show a minor role of 2 / and 2 /ℎ in the high temperature regime. On the other hand, in the transition zone, i.e. mid-range temperature, more notable differences appear in the shape of the response curve. In particular,

S4. Investigation on stress state and bond strength of silicone-wood interface
Supplementary Fig. S12a shows that the constituent materials after the curing process bond at four sites, thus creating four silicone/wood interfaces. It is at these locations that the BB can potentially become weak during deformation. Here we analyze the stress state of their interfaces resulting from an increase of temperature and then measure the bonding strength through a set of pull-out tests.
Section S3-A provides a full analysis of both the temperature-driven deformation and the internal forces generated in a given BB. The results show the existence of an axial compressive force and a bending moment that act on the bonding interface ( Supplementary Fig. S12a). At each of these edges, the distribution of the normal stresses can be simply described as: where is the normal stress exerted by the core to the surrounding low CTE frame, ∈ [−ℎ/4, ℎ/4], and is a given temperature. A positive value of represents a tensile stress, whereas the negative counterpart indicates compression. and are the area and second moment of area of the cross-section. From Eqn.
(39), we gather that the lower half of the cross-section ( ∈ [−ℎ/4, 0]) is subject to either compressive or tensile stress, depending on the magnitude of and , which is governed by temperature. Only the tensile portion can lead to interfacial debonding, which is controlled by the maximum tensile stress that occurs at the furthest edge ( = −ℎ/4) of the cross-section. This critical value varies with the BB core size.

S5. Time response upon heating
As described in S3-A, the deformation of the BB is caused by a mere mismatch of thermal expansion of the constituent materials. The materials do not undergo any atomic or molecular changes, as in the case of smart materials, such as shape memory polymers. On the other hand, the deformation of our BB is determined by the topological layout and the internal forces that are generated in the BB members, with magnitude dependent on temperature. Once the BB has reached the target temperature, the deformation is instantaneous, as opposed to the time response observed in smart materials. With our metamaterials, the only time span involved is that required to heat the sample to the target temperature, which depends only on the heating strategy and the experimental setup.
Our analysis here, therefore, pertains to our experimental system (see S1), which consisted of a heating chamber with multiple thermocouples monitoring the maintenance of 120 ℃, our target temperature. We conducted a series of tests on a BB placed in the middle of the chamber already heated at 120 ℃ and we recorded with a digital camera (EOS Rebel T6i, Canon USA) the time span required for the BB to deform.
Two medium types were examined: fan-propelled air and oil bath, the former requiring 300±17 seconds to complete the deformation, and the latter 115±18 seconds.
To corroborate our experimental observations, a numerical analysis of the transient heat was performed with the following values assigned to the thermal conductivity, , and specific heat, , of the constituents:

S6. Morphing from a pre-assigned sequence of building blocks
S3 has examined the response of an individual BB through theory and numerical analysis. Here we focus on a monolithic assembly of BBs, each with its own set of geometric attributes. The goal is to express the collective deformation of the metamaterial phenotype at a given . As a paradigmatic example, we consider a series of rotational units sharing low CTE edges and all having equal width . Supplementary Fig. S14 shows its deformed configuration, the phenotype at a given . The central curve (dash type) is specified as representative of the deflected axis, and assumed as an arc spline consisting of a number of 1 continuous arcs and straight segments, all passing through the interface mid-points, , between adjacent units and normal to their contiguous edges.
Supplementary Fig. S14: Schematic of a chain of building blocks in series, with an arc connecting a pair of consecutive mid-points at the BB interface. The local coordinate system (X′Y′) of a generic unit is referenced to the global system (XY) on the first outer frame of the first unit.
With reference to the generic building block , we express the coordinates of its point in the local coordinate system, indicated as the X′Y′ plane, and in the global system residing at 0 and visualized as the XY plane in the first unit. The two systems are rotated by and translated by a vector pointing from point to −1 . The coordinates of in the X′Y′ plane can be represented by the two-dimensional array: where the integer ∈ [1, ] and the expression of , which is the distance between two consecutive midpoints of the low CTE edges, is given by: where ℎ , , , are the main geometric parameters controlling the R-BB response, and , 0 are minor parameters here assumed unchanging for all the units due to their minor influence (Supplementary Fig.   S11). represents the opening angle of building block , whose value is governed by ℎ ⁄ and ⁄ , as illustrated in (Supplementary Fig. S9b). A positive value of represents a clockwise rotation, and negative means counter-clockwise.
The following affine transformation can be used to express in the global coordinate system (XY): where 0 is null. As assumed above, the central axis of the deformed configuration of the BB assembly can be represented by a 1 continuous arc spline passing through all the mid-points . The expression of a generic arc in the arc spline can be expressed from point −1 and to point , as: with being the radius, and ( , ) the coordinates of the arc center, whose expressions are:

S7. Morphing on-target through encoded sequence of building blocks
Morphing on target poses the search of a BB sequence that allows the metamaterial phenotype to accurately conform to a target domain. We describe here the target as a domain with central axis representable with an arc spline and two symmetric boundaries that are continuous with varying distance along the central axis. Both are shown in blue on the upper part of Supplementary Fig. S15, the former in dash-dot and the latter with solid line. In red, at the lower part of the figure, is a metamaterial phenotype with randomly assigned BB sequence expressed at a given temperature.
Supplementary Fig. S15: Target domain (blue -top) with varying width and off-target metamaterial phenotype (red -below), both sharing the initial location where the global coordinate system is anchored.

S7-A. Description of target shape
The central axis of the target shape is here assumed to trace an arc spline mathematically expressed piecewisely in an implicit form as: where denotes the number of arcs and straight lines, and the intervals are the coordinates of the blending points.
For the symmetric boundaries of the target shape, the function describing the width change along the central axis is expressed as:

S7-B. Problem formulation
The phenotype will match the target shape if two conditions are guaranteed. The first pertains to the central axis of the metamaterial phenotype. For a given temperature, a generic phenotype with randomly assigned BB sequence will most likely be off-target, i.e. its central axis would appear far from that of the target. This The second condition relates to the boundary of the target domain. Again here, an arbitrary string of BBs will likely lead to a phenotype with non-conforming boundaries, i.e. incompatible from those of the target ( Supplementary Fig. S15). We thus express the gap between the width of each BB and the width of the boundary at a given position as: where = + 2 0 is the width for BB in the off-target phenotype, and ( , ) is the width of the target shape corresponding to point .
In general, the characteristic distances of an off-target phenotype from the target domain, i.e. Eqns. (50) and (51), should be both minimized simultaneously. This can be expressed through a least-squares model that minimizes the sum of the squares of both residuals. In this case, the overall objective function can be formulated as:
The first four are treated as variables, whereas and 0 are assumed constant because (as S3-C demonstrates) the BB response is less sensitive to changes in their value. We can thus collect the most influential descriptors of the BB geometry in the following vector of variables: = {ℎ , , , } T ( = 1,2, . . ) (53)

S7-D. Design constraints
To achieve the kinematic compatibility between neighbouring units, the edge of each building block should be normal to the tangent direction of the target axis, an arc spline, on point . These constraints can be mathematically expressed by the system of equalities: where denotes the normal vector of the target curve on point , expressed as: with / and / representing partial derivatives of the target arc spline; ( , ) denote the coordinates of point in the global coordinate system. is a vector in the axis direction of the localized coordinate system residing at point −1 , and it can be obtained from the following transformation: The unit vector ′ is expressed as [1,0] T in the local coordinate system and the matrix (Eqn. (44)) defines the relative rotation.
The opening angle for each R-BB is governed by both ℎ ⁄ and ⁄ , and Supplementary Fig. S9b presents the maximum range of that R-BB can generate in the geometric space. To produce admissible deformation, the following constraints should not be violated: where /ℎ and /ℎ are the lower and upper bounds on the ratio of ℎ ⁄ , respectively, and / and / are lower and upper bounds on ⁄ .
For practical considerations, side constraints are also set for each design variable, i.e. ℎ ≤ ℎ ≤ ℎ , ≤ ≤ , ≤ ≤ , ≤ ≤ , where the superscripts L and U represent the lower and upper bounds on the design variables, respectively.

S7-E. Optimization problem formulation
The problem of morphing on-target is thus expressed as the minimization of the sum of the squares of both gaps (Eqns. (50) and (51)) subject to the following set of equality and inequality constraints:

S7-F. Sensitivity analysis
The derivatives of the objective function with respect to the design variables ℎ , , and ( = 1,2, … , ) can be obtained through After calculating the gradients above, Eqns. (62), and (65)-(69) are used to guide the search direction of the gradient-based optimization scheme described in the following section.

S7-G. Optimization method
Eqn. (61) describes a nonlinear optimization problem with multiple equality and inequality constraints.
To solve it, we opt for the Powell-Hestenes-Rockafellar (PHR) method ① , an Augmented Lagrangian algorithm well suited for minimization problems with equality and inequality constraints. The Lagrange function defined for the primal problem is reframed with the definition of a set of penalty functions, which ① Rockafellar, R. T. The multiplier method of Hestenes and Powell applied to convex programming. Journal of Optimization Theory and applications 12, 555-562 (1973).
allow transforming the primal problem into a series of unconstrained sub-problems. The Augmented Lagrangian function corresponding to the primal problem in Eqns. (61) is constructed as: where is a vector representing the design variables (Eqns. (53)), and are the Lagrange multipliers, and the penalty multiplier; represents the -ℎ inequality constraint, including ℎ ⁄ , ℎ ⁄ , ⁄ , and ⁄ . In the -ℎ step of the iterative process for solving the primal problem, starting from the design point, , the design point for the ( + 1)-ℎ step is obtained by minimizing the following unconstrained subproblem: To solve Eqn. (71), we resort to the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method ① that updates the design variables such that the updating scheme, in the -ℎ step of the solution process of the sub-problem, is given by: Using the BFGS method, the updating formula for the approximation of the Hessian matrix is of the form where ∆ = +1 − and Δ = +1 − . To determine the appropriate step size during the search for an improved point, whereby the design objective moves towards a feasible descent direction, we define and minimize the following descent function: A sequence of trial step size, , is defined as = (0.5) = 0,1,2,3,4, … At the − ℎ iteration, we determine an acceptable size as = , with as the smallest integer to satisfy the descent condition The constant is determined using the search direction = ‖ ‖ 2 , where is a specified constant between 0 and 1.
After solving the subproblem Eqn. (71), the design point +1 is obtained for the ( + 1)-ℎ step in the iterating process solving the primal problem, while the multipliers are updated through the following rule: is a specified constant between 0 and 1, and is larger than 1.
The design problem in Eqns. (61) is solved iteratively until the following set of convergence criteria is satisfied: where 1 , 2 and 3 are convergence control parameters taking respectively the value of 0.001, 0.001, and 0.01.

S8. Implementation of morphing on target scheme
This section shows the implementation of the inverse problem scheme explained in S7 to the solution of the last two illustrative demonstrations of morphing on target. While they both share a primary objective, the central axis of the target domain (i.e. an "M"), their boundaries are distinct. The first assumes a uniform, yet unspecified value of the width, which is to find, and the second case features a varying width to match.
In both cases, the constituent solids are silicone rubber and wood with CTE and elastic properties reported in section S2.

S8-A. Demonstration 1: uniform width
As a morphing target, we chose the profile of an "M", and use an arc spline composing smooth, 1 continuous, arcs to prescribe its central axis ( Supplementary Fig. S16). Because = 7.4 is a vertical axis of symmetry, only half "M" is here examined. The width of the "M" boundary is uniform along the "M", meaning that 1 + 2 0 = 2 + 2 0 = ⋯ = + 2 0 for all BBs. The central axis of the half domain is here mathematically expressed piece-wisely by the following set of primitives ( Supplementary Fig. S16 Due to the randomly assigned values of the design variables, the central axis of the metamaterial phenotype is off from the target. After application of the tailoring scheme, though, the phenotype is realigned to the target after 80 iterations. Supplementary Fig. S18 shows the iteration history for the Supplementary Fig. S19 shows, in a stem plot, the optimized values of the design variables along with their relevant ratios, which define the metamaterial genotype at the initial temperature ( = 20 ℃). The values shown in sequence for the 50 BBs are descriptive of each BB geometry as well as their sequence.
As described in the caption of Fig.1, each building block is described with , which can be further condense to Supplementary Fig. S20 shows the resulting configuration at = 120 ℃ along with the phenotype axis and the target axis. The similarity between the phenotype axis and the target axis is evaluated by calculating the coefficient of determination 2 , which is mathematically defined as where ( , ) and ( , ) represent the coordinates of the BB interface point i of the phenotype on target and the corresponding reference point on the target axis respectively. ̅ and ̅ denote the mean value of the relevant coordinates for BB interface points. The coefficient of determination indicates the precision of the fit for both -and -values and ranges from 0 to 1. The better the match between phenotype and target, the closer the value of 2 to 1. In our example, the value of 2 is 0.994, meaning a remarkably good agreement between target and predicted phenotype axis. A minor discrepancy appears (Supplementary As per Eqn. (53), we look for the optimum values of the design variables ℎ , , , ( = 1,2, … , ).
Since both the axis and width functions of the target are specified, the problem is formulated for a given The optimization is implemented under the following initial values of the design variables: ℎ = 0.1, = 0.1, = 0.7, = 0.35, = 1,2, … ,46. The tailoring process gradually converges to an optimal solution after 128 iterations, and Supplementary Fig. S22 shows the iteration history for the objective function. At convergence, the objective value decreases from the initial 2024 to the final 1.82 × 10 −5 , an extremely small value indicating the achievement of the gap closure, where the phenotype axis and its width conform seamlessly to the axis and boundary of the target.
Supplementary Fig. S22: Convergence plot for the objective function.