Optimal face-to-face coupling for fast self-folding kirigami

Kirigami-inspired designs can enable self-folding three-dimensional materials from flat, two-dimensional sheets. Hierarchical designs of connected levels increase the diversity of possible target structures, yet they can lead to longer folding times in the presence of fluctuations. Here, we study the effect of rotational coupling between levels on the self-folding of two-level kirigami designs driven by thermal noise in a fluid. Naturally present due to hydrodynamic resistance, we find that optimization of this coupling as control parameter can significantly improve a structure's self-folding rate and yield.

Thermal fluctuations are critical to how microscopic systems explore their configuration space and converge to the desired target structures [25,[36][37][38]. Due to these fluctuations, the folding trajectories are stochastic and the final configuration might not coincide with the desired one. The trajectories depend on the choice of the initial template [25], on the properties of the material [21], and on the experimental conditions [36]. Folding yield has therefore been considered an important parameter to optimize [22,25]. Nonetheless, folding time is an equally key parameter beyond yield for real-life applications of microscopic kirigami designs, yet its optimization is much less understood. An accurate prediction of the folding time was possible for simple single level structures (e.g. pyramids) [37,38]. However, kirigami designs often present multiple interdependent levels connected via joints, whose folding could be strongly affected by level-to-level correlations, e.g., due to materials or environmental properties [34,35,39,40].
Here, we demonstrate how rotational coupling between levels emerges naturally in microscopic hierarchical kirigami templates folding in a fluid. We numerically show how this coupling can be tailored to minimize the folding time and rationalize the choice of the optimal coupling parameter by mapping our results into a first passage problem.
As model target structures, we consider double pyramids fixed on a substrate with three sets of two hinged lateral faces of constant height h ( Fig. 1a-b). Two target angles define the three-dimensional geometry of the structures (Fig. 1a): φ l ∈ (0, π) (defined with respect to the substrate) and φ u ∈ (−π, π) (defined with respect to the plane containing the lower face). Based on the choice of angles, the two-dimensional template (gray area, Fig. 1b) spans from an hourglass shape to a diamond and is obtained by cutting the edges of the lateral faces of the target structure and unfolding them. The upper level faces are hinged to the respective lower level faces, which are in turn tethered to the substrate via a second hinge. Two time-dependent angles, θ l (t) and θ u (t), describe the motion of the two hinges during the folding process driven in the fluid by thermal fluctuations (here, in water at room temperature).
A model for the coupled motion of two faces joined by a hinge can be derived using the approach of 'resistive-force theory' [42,43]. Under this approach, the hydrodynamic drag on each face combines with constraints of vanishing net torque and force (owing to the absence of inertia) to yield a coupled equation of motion, while secondary hydrodynamic interactions between the faces are neglected (Supplemental Material for full derivation [41]). Thus, where F i l and F i u are the magnitudes of the forces acting in the normal direction to the lower and upper faces ( Fig.  1c and Fig. S1), respectively, C < 0 is the normal hydrodynamic resistive force coefficient for each face, which we assume here to be equal (Supplemental Material for the more general case [41]), andα < 0 andβ < 0 provide negative feedbacks (i.e. contrary motions) between the faces. Equation 1 shows how the motion of the two faces joined via a hinge is thus naturally coupled due to the hydrodynamic resistance of the fluid. In fact, for small upper angles θ i u , the feedback terms in Eq. 1 become roughly constant and are dominated byα, which is 15 times larger thanβ in this limit. Indeed,α remains roughly an order of magnitude larger thanβ over a wide range of θ i u (Supplemental Material [41]), reflecting the fact that the upper face is much more strongly affected by the motion of the lower than the lower is by the upper.
For microscopic structures, the folding is driven by thermal fluctuations, and the driving forces (F i l and F i u ) change rapidly, being well described by a stochastic process uncorrelated in space and time [44]. The variance of this process, and hence the typical rotational diffusion coefficients D θ for each face, can thus be derived assuming equipartition (Supplemental Material [41]). Assuming a hinged circular disc of radius h/ √ π, we obtain D θ = 3kBT 8µh 3 , where k B is the Boltzmann constant, T the thermostat temperature, and µ the fluid viscosity (Supplemental Material [41]). For h being O(µm) to O(10 µm), D θ varies from O(10 rad 2 s −1 ) to O(10 −3 rad 2 s −1 ).
In order to primarily explore the role of the dominant lower-to-upper coupling of the faces, and taking the sheets as being driven by thermal fluctuations alone, we consider here a reduced model of the forṁ where D θ l and D θu are the rotational diffusion coefficients of the lower and upper face, respectively, and η i θ l and η i θ l are independent white noise processes [44]. Here, the weaker coupling is neglected and the dominant coupling is parameterized by a single dimensionless constant α; a reference value of α is obtained in the limit θ i u → 0, which gives α H ≡α(θ i u → 0) = −5/2. This follows directly from the model in Eq. 1 but, more generally, we leave α as a free parameter to allow for this feedback to be 'engineered' or controlled by other constraints not included in the simple model that underlies Eq. 1. Indeed, even within that modeling framework it is possible to change the strength and direction of coupling by allowing variability in the properties of the two faces and the hinge point between them (Supplemental Material [41]).
For each set of coupled faces i (Fig. 1), we solve these differential equations numerically for D θ l = D θu = 0.64 rad 2 s −1 using a first-order (Euler) integration scheme with a short timestep (∆τ = 10 µs) [44]. We implement the interaction with the substrate as reflective boundary conditions [44] and we correct for collisions among the faces by detecting them with the Gilbert-Johnson-Keerthi (GJK) algorithm [45]. We set a cutoff time τ cut = 2 · 10 8 ∆τ for the simulations, when we consider misfolded any structure which is not folded completely. This cutoff time was chosen to guarantee that the rate at which new structures fold is near zero above it. Sample trajectories are shown in Fig. 1d-e for the faces of the upper and lower levels when α = α H , and show how the system converges to the target structure through a series of four irreversible binding events between faces of the lower level first followed by the upper level, where each event is defined by two faces being at the respective target angle (±π/180) at the same time. Figure 2 shows that there is scope to engineer the coupling to optimize the folding process by tuning the valueand sign -of the coupling parameter α. This could be achieved by, e.g., engineering the hydrodynamic resistance of the sheets or exploiting active or driven hinges (Supplemental Material [41]). These observations are qualitatively independent of the exact choice of the diffusion coefficient (Fig. S3). For example, Fig. 2a shows that, for the structure in Fig. 1 with an obtuse lower level's target angle and a negative upper level's target angle, the folding rate is enhanced by a factor of ≈ 1.3 at a slightly negative α (α = −1) when compared to the rate without coupling (α = 0) and marginally (by a factor of ≈ 1.1) when compared to the rate at α H . Instead, for a different target structure such as a diamond with an acute lower level's target angle and a positive upper level's target angle (Fig. 2b), positive α values can enhance folding by a factor up to ≈ 6 when compared to the rate at α H and up to ≈ 4.5 when compared to the no-coupling case. Interestingly, for both structures, the α value that optimizes the rate (α opt ) is closely related to the value that optimizes yield (defined as the percentage of fully folded structures within the cutoff time τ cut ), with a yield of ≈ 70% (against the highest yield of ≈ 74%) and ≈ 61% (against the highest yield of ≈ 62%) at α opt in Fig.  2a-b, respectively. Figure 2 suggests that, by correctly designing the coupling between faces, we can optimize both folding rate and yield of the kirigami at once. This observation can be generalized to a variety of structures obtained by varying the two target angles φ l and φ u (Figs. 3 and S4). Figure 3a shows some examples of structures. To a first analysis, the value of φ u is the main decisive factor discerning whether a negative (for φ u < 0) or positive α (for φ u > 0) optimizes the structure folding rate (Fig. 3b). In general, structures with upper levels opening away from the center (e.g. II, IV and VI in Fig. 3b) tend to benefit from a negative coupling to their respective lower levels, while structures whose upper levels point towards their centers (e.g. I in Fig. 3b) tend to fold faster with a positive coupling. The absence of coupling instead tends to be optimal for structures (e.g. III in Fig. 3b and Fig. S4a) where the lower and upper levels roughly lie on the same plane (φ u ≈ 0). In fact, Fig. 3c shows that the benefits in folding rate due to the presence of coupling between levels increases the further a target structure is from this condition, as k α ≈ k 0 around φ u ≈ 0. This result highlights the importance of engineering the coupling parameter even more, as the folding of target structures in this range of relatively smaller upper angles (Fig. S2) is otherwise inhibited by sub-optimal coupling due to the natural hydrodynamic resistance of the hinged faces (Supplemental Material [41]). The range of φ u values where no coupling is advantageous broadens asymmetrically towards positive α values as φ l increases towards π. In these situations, stronger coupling tends to push the upper level's faces against the substrate, thus slowing down the convergence to the target structure. As already noted for two exemplary target structures (Fig. 2), the value of α that optimizes folding rate also roughly optimizes yield (Fig. S4b-d). In fact, for most structures, the ratio between the yield achieved at α opt and the maximum yield at any α (max(yield), Fig. S4b) is close to one (Fig. S4c), and the distance between these two values of α is often close to zero (|α max(yield) − α opt | ≈ 0, Fig. S4d). Larger separations between these two α values are possible but often coincide with regions where the yield is relatively insensitive to the exact α value (Fig. S4c) or where advantages over no coupling are negligible (Fig. 3c).
To rationalize these observations on the emergence of an optimal coupling parameter α and to provide a design rule for its choice, we can map our results into a first passage problem to a target (Supplemental Material [41]) [46]. This choice is justified by recent work that has demonstrated how, for single-level pyramids with less than five faces, the folding time is dominated by the closing of the first pair of faces -event that can be described as a first passage problem in a two dimensional random walk [37]. Thus, the total folding time for single-level structures depends on the initial location of its faces [37]. Similarly, for two-level structures, if we assume that the lower level folds much faster than the upper one, the closing of the first pair of upper level's faces after the lower level has closed completely (t 3 in Fig. 1) is the leading event that dominates the folding of most structures. The optimal value of α should then be the one that leads to the distribution of angles in the upper level when the lower level has closed (t 2 in Fig. 1) that minimizes t 3 (i.e. the time for the first two faces of the upper level to reach their target angle φ u ). Starting from a location drawn from a one-dimensional Gaussian distribution N (µ, σ) with varying mean µ(α) and fixed variance σ and reflected at the boundaries to ensure its proper normalization (Supplemental Material [41]), the mean first passage time to a trap at location φ φ φ = (φ u , φ u ) of a two-dimensional random walk on a lattice with reflective boundaries can  (Fig. 1a). As φ l increases from acute to obtuse, the target structure of the lower level (gray faces) transitions from an inverted pyramid (I, II, IV) as in Fig. 2b to a regular pyramid (VI) as in Fig. 2a via truncated pyramids (III, V). For φu going from negative to positive values, the upper level goes from an umbrella shape, either convex (IV) or concave (VI), to a pyramid (I, III, V), truncated or not, through a flat plane (II). (b) Phase diagram of the optimal coupling parameter αopt for the different structures. The colour indicates whether a negative (blue) or positive (red) α is optimal (Fig. S4a). (c) Highest folding rate max(kα) at any α for different target structures. Folding rates are normalized to k0. In (b-c), the black dashed isolines show the optimal values of α that lead to the highest folding rates and the black solid lines highlight structures whose folding is optimal in the absence of coupling (α = 0). Each data point is obtained as an average of 5000 folding events. be calculated as [46] where T θ0→φ φ φ is the first passage time to the trap starting from position θ 0 = (θ 01 , θ 02 ) on the lattice with probability P (θ 0 ) = N (µ, σ) |θ01 N (µ, σ) |θ02 (Supplemental Material [41]). Figure 4a shows that τ f (µ) should be minimized when α is chosen so that µ(α) ≈ φ u . Our numerical results indeed confirm that, for a given value of α, the means of the distributions of the angles in the upper level after to lower has closed completely (exemplary distributions shown in Fig. 4b) as a function of the lower level's target angle φ l closely resemble the upper level's target angles φ u of the structures with the highest folding rate at that α (Figs. 4c and S5), thus effectively providing a design rule for the choice of the optimal coupling parameter α between levels based on the target geometry.
In conclusion, we have shown how the folding of microscopic hierarchical two-level kirigami structures in a fluid develops a natural degree of coupling between levels due to hydrodynamic resistance. This coupling is an essential factor to account for when it comes to understanding the folding dynamics of similar structures as both folding rate and yield are impacted in its absence. Moreover, we show how rational design rules (e.g. based on solving a first passage time problem) can be used to explain the emergence of an optimal coupling parameter for maximizing folding rates and yield as a function of the geometry of the target structure. In the future, we envisage that these design rules could be adopted to engineer hydrodynamic resistances and active hinges for designing fast self-folding hierarchical multi-level kirigami structures for applications in, e.g., soft robotics [16][17][18][19][20] and mechanical actuators and metamaterials [9,10,[13][14][15].
Acknowledgements. MPB, QXP and GV are grateful to the studentship funded by the A*STAR-UCL Research Attachment Programme through the EPSRC M3S CDT (EP/L015862/1). NAMA acknowledges financial support from the Portuguese Foundation for Science and Technology (FCT) under the contracts no. UIDB/00618/2020 and UIDP/00618/2020. BJW is supported by the Royal Commission for the Exhibition of 1851. NAMA and GV acknowledge support from the UCL MAPS Faculty Visiting Fellowship programme.  Fig. 1 for α = −1, α = 0 and α = 1. The vertical dashed lines represent the means µ(α) = θ i u (t2) α of the same-color distributions. Each distribution is obtained from 50000 structures. (c) For a given value of α and lower level's target angle φ l (see also Fig. S5), the means of these distributions (solid lines) are well captured by the range of possible structures with the highest folding rate at that α value (Fig. 3c) once the numerical uncertainty in determining max(kα) is accounted for (shaded areas). For a given structure, this uncertainty around the position of max(kα) is calculated taking the largest between the discretization step in α and the range of data points around the peak falling within the standard error.