Thermal buckling and symmetry breaking in thin ribbons under compression
Paul Z. Hanakata, Sourav S. Bhabesh, Mark J. Bowick, David R. Nelson, David Yllanes
TThermal buckling and symmetry breaking in thin ribbons under compression
Paul Z. Hanakata a, ∗ , Sourav S. Bhabesh b,1 , Mark J. Bowick c , David R. Nelson d , David Yllanes e,f,c a Department of Physics, Harvard University, Cambridge, MA 02138, USA b Amazon Web Services (AWS), Washington DC Metro Area, USA c Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA d Department of Physics, Lyman Laboratory of Physics and School of Engineering and Applied Sciences, Harvard University, Cambridge,MA 02138, USA e Chan Zuckerberg Biohub, San Francisco, CA 94158, USA f Instituto de Biocomputaci´on y F´ısica de Sistemas Complejos (BIFI), 50009 Zaragoza, Spain
Abstract
Understanding thin sheets, ranging from the macro to the nanoscale, can allow control of mechanical properties such asdeformability. Out-of-plane buckling due to in-plane compression can be a key feature in designing new materials. Whilethin-plate theory can predict critical buckling thresholds for thin frames and nanoribbons at very low temperatures, aunifying framework to describe the effects of thermal fluctuations on buckling at more elevated temperatures presentssubtle difficulties. We develop and test a theoretical approach that includes both an in-plane compression and an out-of-plane perturbing field to describe the mechanics of thermalised ribbons above and below the buckling transition. Weshow that, once the elastic constants are renormalised to take into account the ribbon’s width (in units of the thermallength scale), we can map the physics onto a mean-field treatment of buckling, provided the length is short comparedto a ribbon persistence length. Our theoretical predictions are checked by extensive molecular dynamics simulations ofthin thermalised ribbons under axial compression.
Keywords:
1. Introduction
Thin sheets, possibly with embedded kirigami cuts,have been the object of intense recent study [1]. A carefuldesign allows membranes with cuts to stretch far beyondtheir pristine limits [2–7], to have non-linear post-bucklingbehaviours [8, 9], and even to exhibit complex motionssuch as roll, pitch, yaw, and lift [10]. Many of these noveleffects arise due to out-of-plane deflections, i.e., escape intothe third dimension. With such mechanical versatility andstraightforward actuation, kirigami sheets have been usedas building blocks for soft robots, flexible biosensors andartificial muscles [11, 12]. A full theoretical framework forthis rich phenomenology must rest on a thorough under-standing of the fundamental mechanical effects. In partic-ular, out-of-plane motion in simple kirigami systems (e.g.,a sheet with a single slit) have been described as an Eulerbuckling problem [10]. The buckling of pillars and plateshas been studied for centuries, but a unifying theory tounderstand buckling in nanosystems when thermal fluctu-ations become important, as in the case of molecularly thinmaterials such as MoS and graphene [13], is still lacking.In the classical description, the dimensionless F¨oppl-von K´arm´an number vK = Y W L /κ , where Y is the 2D ∗ Corresponding author.
Email address: [email protected] (Paul Z. Hanakata) Work completed prior to joining AWS.
Young’s modulus, κ is the bending rigidity, W and L are respectively the T = 0 width and length of the rib-bon, can be used to quantify the ease of buckling a thinsheet out of plane at zero temperature. The picture ismore complicated for thermalised membranes [14], where Y and κ become scale dependent and, in particular, thebending rigidity is dramatically enhanced [15–20]. Thislongstanding theoretical prediction is consistent with animportant study of graphene ribbons by Blees et al. [2].Using a cantilever setup, the effective bending rigidity ofmicron-size graphene at room temperature was found toincrease by a factor of roughly 4000 relative to the zero-temperature microscopic value. Although it is possiblethat some of this increase may be due to quenched randomdisorder in the graphene ribbons [21], these room temper-ature experiments nevertheless demonstrate a striking en-hancement over the T = 0 density functional theory pre-dictions [22]. When thermal fluctuations are important,classical Euler buckling predictions break down. In fact,in such an entropy-dominated high-temperature setting,some aspects of nanoribbon behaviours have more in com-mon with linear polymers with long persistence length [23].In this letter, we investigate (i) to what temperatureclassical Euler buckling still holds, (ii) how we can lo-cate buckling transitions in fluctuating ribbons under com-pression, and (iii) how these buckling transitions changewith temperature and with the ribbon dimensions. To Preprint submitted to Extreme Mechanics Letters December 14, 2020 a r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec his end, we develop a mean-field theory (MFT) approachto the buckling of thermalised ribbon under longitudinalcompression and use molecular dynamics simulations tocheck our predictions. The applicability of our MFT isdetermined by two crucial length scales: First, the ther-mal length (cid:96) th ∼ κ/ √ Y k B T , where k B is the Boltzmannconstant, T is the temperature, Y and κ are the micro-scopic 2D Young’s modulus and bending rigidity respec-tively. And second, the one-dimensional persistence length (cid:96) p = 2 κW /k B T . We are interested in the regime (cid:96) th 2. Model and methods We consider a rectangular sheet of size L × W , with L > W , which is discretised by a triangular lattice of un-breakable bonds, in the crystalline membrane paradigm [26].The triangular lattice used here can be considered as aconvenient dual representation to the honeycomb latticeusually employed to model graphene. We use the notation L to distinguish the T = 0 rest length from the projectedlength after thermal shrinking or compression. Neighbour-ing nodes are connected by harmonic springs and there isan energy cost when the normals ( nnn α ) of neighbouring pla-quettes are not aligned. The total energy is given by H = k (cid:88) (cid:104) i,j (cid:105) ( | rrr i − rrr j | − a ) + ˆ κ (cid:88) (cid:104) α,β (cid:105) (1 − nnn α · nnn β ) (1)where k is the harmonic spring constant, ˆ κ is the micro-scopic bending rigidity and a is the preferred length be-tween two neighbouring nodes which also sets our unitof length. The first sum is over neighbouring nodes andthe second over neighbouring triangular plaquettes. Aschematic is shown in Fig. 1(a). Our discretised bare elas-tic constants are related to the bare continuum ones by κ = √ κ/ Y = 2 k/ √ n α n β a (a) L relax W ε =0, T=0 ε > ε c , T>0 L ε (b) (c) ε =0, T>0 L (d) Figure 1: (a) Schematic of nodes on a triangular lattice. (b) Ribbonin its flat rest configuration at T = 0, with length L and width W .(c) Relaxed ribbon at T > L down to athermal projected length L relax . (d) Ribbon at a non-zero tempera-ture and compressed beyond the critical buckling strain (cid:15) c . In panels(b)–(d) each end of a ribbon is clamped to have width W and thecolour map shows the height of each node relative to z = 0 at thetwo ends. We used OVITO software to visualise the ribbons [25]. Since we are interested in relatively narrow ribbons,we use L ∼ a and W ∼ a (2500 nodes). Followingprevious work [28–30], we set k = 1440ˆ κ/a , which givesus a F¨oppl-von K´arm´an number of vK ∼ , comparableto micron-size 2D materials such as graphene and MoS .As we change the temperature, keeping k/ ˆ κ fixed, twocrucial length scales, the thermal and persistence lengths,will vary [23, 29]: (cid:96) th = (cid:115) π κ k B T Y , (2) (cid:96) p = 2 κW k B T . (3)We want here to adapt the zero-temperature theory totemperatures high enough for thermal renormalisation tobecome significant. The temperature should not, however,be so high that (cid:96) p becomes small compared to L (i.e., westay far away from the ribbon crumpling regime). In sim-ulations we fixed L , W and k/ ˆ κ while varying ˆ κ and T .We simulated over a temperature range 10 − ≤ k B T / ˆ κ ≤ − (cid:46) W /(cid:96) th (cid:46) . In the following weshall use W /(cid:96) th as the natural variable for the tempera-2ure scaling of the system, and focus on the regime where W > (cid:96) th . We use the HOOMD package [31] to simulate model (1)in the NVT ensemble with a Nos´e-Hoover thermostat. Inorder to study the buckling dynamics, we clamp the ribbonby fixing the nodes on the first two rows at both ends. Wevary the distance between the clamped edges to induce thedesired strain. Importantly, we thus operate in a constant-strain ensemble.Because of thermal fluctuations, the ribbon shrinksfrom from its T = 0 rest length L . We define L relax asthe projected natural length at which all stress componentsare zero and define the incremental compressive strain as (cid:15) = 1 − L (cid:15) /L relax , where L (cid:15) is the projected length at agiven compressive strain (cid:15) . At finite T we have thereforethe inequalities L (cid:15)> < L relax < L , illustrated in Fig. 1.Following [29], we use a timestep of ∆ t = 0 . τ where τ is the Lennard-Jones time τ = (cid:112) ma /k B T andwe use natural units of mass and energy m = a = 1. Ourclamped systems are simulated in the NVT ensemble for10 steps, saving a snapshot every 10 steps. For eachchoice of parameters, we simulate either 5 or (more com-monly) 10 independent runs. We use a jackknife method(see, e.g., [32]) to estimate statistical errors. 3. Theoretical expectations The most dramatic signature of the buckling transi-tion occurs in stress-strain curves. Fig. 2 shows the stressas a function of the strain as measured from our simula-tions when T = 0, then at a low T such that W /(cid:96) th =0 . 3, and finally at a more elevated temperature such that W /(cid:96) th = 8 . 5, where thermal fluctuations have a strongereffect. The computed Young’s modulus, critical stress, andcritical strain for T = 0 are within 10% of the theoreticalpredictions [ Y simulation Y theory = 0 . , σ simulationc σ theoryc = 0 . , (cid:15) simulationc (cid:15) theoryc =0 . W , with displacements uniform along the y direc-tion, as E/W = 12 Y (cid:90) L/ − L/ (cid:18) d u ( x )d x (cid:19) d x − σ xx d + 12 κ (cid:90) L/ − L/ (cid:18) d h ( x )d x (cid:19) , (4)where u ( x ) is the displacement field along the x axis, h ( x )the displacement perpendicular to the ribbon, and − σ xx d represents the work done by a force F = W σ xx to com-press the ribbon an amount d along ˆ x relative to its natu-ral length L . Here, Y and κ are the 2D Young’s modulus -0.5 0 0.5 1 1.5 2 2.5 ε / ε ctheory σ / σ c t heo r y Simulations at T =0-5 0 5 10 15 20 ε / ε ctheory -3-2.5-2-1.5-1-0.500.511.522.533.54 σ / σ c t heo r y W / l th ~0.3 0 1000 2000 3000 ε / ε ctheory -4-2024681012 W / l th ~8.5 (a) (b) (c) d d -0.6 -0.4 -0.2 0 0.2 0.4 0.6-101234567 -0.4 -0.2 0 0.2 0.400.511.52 Figure 2: (a) Stress as a function of compressive strain scaled by theirzero-temperature critical values ( σ theoryc and (cid:15) theoryc ) for a ribbon at(a) T = 0, (b) W /(cid:96) th = 0 . W /(cid:96) th = 0 . 3. The dottedlines are linear fits in the pre-buckling (small- (cid:15) ) regime and in thepost-buckling regime (cid:15) > (cid:15) c . The scaled critical stress, which isproportional to the renormalised bending rigidity κ R , increases withincreasing W /(cid:96) th (or temperature). In contrast, the slope ( Y R )becomes smaller with increasing W /(cid:96) th . Note that the very differenthorizontal scales in (b) and (c). and bending rigidity which measure the compressional andbending energies respectively. In the compressed, but un-buckled, state the strain is (cid:15) = d u ( x )d x = d/L and fromEq. (4), the compressional energy is E comp = W Y d /L .In this regime, we minimise over d to find Hooke’s Law σ = F/W = Y (cid:15) , which accounts for the first, linear partof the stress-strain curve. Beyond the critical strain (cid:15) c ,however, the system prefers to trade compressional energyfor bending energy. As we shall discuss below, for tangen-tial boundary conditions at two ribbon ends, as is the casefor our simulations, we have the usual buckling instabilitywhen (cid:15) > (cid:15) c , (cid:15) c = 4 π κ/ ( Y L ) [33].What is the incremental stress δσ associated with anadditional strain δ(cid:15) = δd/L when (cid:15) > (cid:15) c ? To this end, weassume the compressional energy vanishes. We can nowregard x as a coordinate embedded in the ribbon. Notethat the tipping angle θ ( x ) of the normal away from the z -axis is given by θ ( x ) ≈ d h d x , so that the additional energyassociated with the buckled state can be rewritten as∆ E/W ≈ κ (cid:90) L/ − L/ (cid:18) d θ ( x )d x (cid:19) d x − δσd, (5)3here − δσW d is the extra work done beyond the bucklingtransition by the stress increment δσ . Once buckling leadsto a ribbon with a well-developed looping arch, i.e., (cid:15) (cid:29) (cid:15) c we expect that d θ ( x )d x ∼ πdL so that the normal turns anangle δθ ∼ π/ d ∼ L/ 2. The energy associatedwith Eq. (5) is then ∆ E ∼ W κd /L . Upon minimisingthis expression over d , we obtain δσ xx = const. κdL ≈ const. κL δ(cid:15). (6)Thus, the slope of the stress-strain curve beyond (cid:15) c , oncethe buckling transition becomes well developed, should beof the order κ/L , as might have been guessed from di-mensional analysis.We conclude that the ratio of the pre- and post-bucklingslopes is ∼ Y L /κ , i.e., it is of the order the F¨oppl-vonK´arm´an number ∼ in our simulations! Hence, it is notsurprising that the zero-temperature stress-strain curvelooks nearly flat in Fig. 2(a). There is, however, a hintof a non-zero slope at finite temperatures in our simula-tions when W /(cid:96) th = 0 . 3, which becomes more pronouncedwhen W /(cid:96) th =8.5. As discussed below, we attribute thisenhanced post-buckling slope to a strong W -dependentupward renormalization bending rigidity κ → κ R , due tothermal fluctuations. Moreover, by rescaling the stressand the strain with their respective zero-temperature crit-ical buckling compression and strain, we can see that thecritical strain and critical buckling compression increasewith increasing T , or equivalently increasing W /(cid:96) th , asshown in Fig. 2(b) and (c).The argument above cannot tell us the details of whathappens close to (cid:15) c , where one must account for deli-cate balance between compression and bending energies.To understand this regime, we now construct a simpleLandau-like theory of the buckling transition, appropriateto the constant-strain ensemble enforced by our constantNVT simulations. As the ribbon is compressed along the longitudinal x direction it can both compress and deflect out of plane inthe z direction. We work in the Monge representation anddenote the vertical displacement by h ( x, y ). In this deriva-tion we denote the instantaneous projected length aftera compression d (to produce a dimensionless compressivestrain (cid:15) ) by L (cid:15) . To control the buckling order parameter,we also impose an out-of-plane electric field E coupled tothe height of a charged ribbon, generating a potential en-ergy V ⊥ = − (cid:82) L (cid:15) / − L (cid:15) / ρ E h ( (cid:126)x ) d x , where ρ = Q/ ( L W )is the charge density. To describe a ribbon in a gravita-tional field we simply need to substitute ρ = m/ ( L W )and E = g . We assume a large F¨oppl-von K´arm´an num-ber Y L W /κ (easily achieved for graphene and Mo S ), Both the critical strain (cid:15) c ∼ κ/ ( Y L ) and the post-bucklingslope ∼ κ/L vanish in the thermodynamic limit L → ∞ . in which the stretching along the ribbon will be compara-tively small. The total free-energy cost is given by G = 12 (cid:90) d x (cid:104)(cid:0) ∇ h (cid:1) + 2 µu ij + λu kk (cid:105) − ρ E (cid:90) d xh − σ (cid:90) d x ( ∂ x u x ) , (7)where u ij = ( ∂ i u j + ∂ j u i ) / ∂ i h )( ∂ j h ) / σ xx = σ de-notes a uniaxial stress at the clamped edges. Notice that,since the centre-of-mass height is h C M = W L (cid:82) d xh ( (cid:126)x ),we can write G = G −E Qh CM and the thermally averagedcentre-of-mass height h CM in the full fluctuating which weare only approximating here (cid:104) h CM (cid:105) = 1 Z (cid:90) D [ h, u i ] h CM e − ( G −E Qh C M ) /k B T , (8)where Z = (cid:82) D [ h, u i ] e − E/k B T is the partition function.Since we are interested in the buckling response due toan external field, we also study the height susceptibilitydefined as χ = d (cid:104) h C M (cid:105) / d E . Upon using Eq. (8) we obtain χ ∝ (cid:104) h (cid:105) − (cid:104) h CM (cid:105) . (9)We can further simply the physics into a 1D buck-ling problem. We approximate h ( x, y ) ≈ h ( x ) and de-fine charge density ρ = Q/L , an effective 1D bendingrigidity and Young’s modulus given by κ = κW and Y = Y W , respectively, where κ and Y denote T = 0 val-ues of the elastic constants. Within a Monge representa-tion, we can approximate L (cid:15) + d (cid:39) L (cid:15) + (cid:82) L (cid:15) / − L (cid:15) / (cid:0) d h d x (cid:1) d x ,where the strain (cid:15) is given by (cid:15) = d/L relax . The total en-ergy then consists of bending, stretching and work doneby the external compressive force F and an out-of-planefield, G [ h, E ] = κ (cid:90) L (cid:15) / − L (cid:15) / d x (cid:18) d h d x (cid:19) + Y L (cid:15) (cid:34)(cid:90) L (cid:15) / − L (cid:15) / d x (cid:18) d h d x (cid:19) (cid:35) − F (cid:90) L (cid:15) / − L (cid:15) / d x (cid:18) d h d x (cid:19) − ρ E (cid:90) L (cid:15) / − L (cid:15) / h d x. (10)Note that we have eliminated, or “integrated out”, thein-plane phonons. See Appendix C for a detailed deriva-tion of Eq. (10), which incorporates our constant-strainboundary conditions. Note also the non-local character ofthe second, stretching term. Lifshitz and Cross [34] havedescribed equations of motion for micro-electromechanicaldevices with a similar non-local term. The ansatz of thefirst buckling mode h ( x ) = h M (cid:104) (cid:16) πxL (cid:15) (cid:17)(cid:105) , whichallows for a height h M midway between the clamps andsatisfies the boundary conditions d h d x (cid:12)(cid:12) x = ± L (cid:15) / = 0, then4eads to an expansion in the buckling amplitude h M G = π L (cid:15) (cid:18) κ π L (cid:15) − F (cid:19) h + π Y L (cid:15) h − ρL (cid:15) E h M . (11)Note that, although Eq. (11) resembles a Landau theorynear a critical point, the expansion parameter depends ina non-trivial way on the system dimension L (cid:15) . Note alsothat the single mode approximation only makes sense closeto the transition; many more Fourier modes would be re-quired to describe the fully developed post-buckling loop-ing arch that develops for large strains, as in Fig. 2(a). T = 0For E = 0, we can minimise Eq. (11) over h M to obtaina critical 2D compressive stress of σ c = 4 π κ/L (cid:15) c , and acorresponding critical buckling strain (cid:15) c = σ c /Y = π κY L (cid:15) c ,where L (cid:15) c is the projected length at the critical bucklingstrain. These are the critical load and critical strain ofclassical Euler buckling with tangential boundary condi-tions [33]. The buckling amplitude is then h M = (cid:40) , (cid:15) < (cid:15) c (or σ < σ c ) , ± L (cid:15) c π (cid:113) (cid:15) − κπ Y L (cid:15) c , (cid:15) ≥ (cid:15) c (or σ ≥ σ c ) . (12)To test the above approach, we compared simulations at T = 0 with the analytical predictions. These simulationsreproduced the square-root scaling predicted by the the-ory and yielded consistent values for the Young’s modulus,critical stress and critical strain (see Appendix A for de-tails and plots). At the critical point the system becomes sensitive toexternal perturbation. In analogy with the magnetic sus-ceptibility of an Ising system, within the MFT we candefine a height susceptibility as the linear response to auniform out-of-plane external field, χ ≡ ∂h M ∂ E (cid:12)(cid:12)(cid:12)(cid:12) E =0 = (cid:40) QY π L (cid:15) c W ( (cid:15) c − (cid:15) ) − if (cid:15) < (cid:15) c Q Y π L (cid:15) c W ( (cid:15) − (cid:15) c ) − if (cid:15) > (cid:15) c . (13)This response function diverges at the buckling transition.Hence, the system becomes infinitely sensitive to the out-of-plane field E as the buckling transition is approached.Note also that χ is larger for small aspect ratio W /L (cid:15) .Eq. (11) predicts a non-linear dependence of the bucklingamplitude h M on E when (cid:15) = (cid:15) c h M = L (cid:15) c (cid:18) Q E π Y W (cid:19) / . (14)The finite-temperature generalisation of this susceptibilityis given by Eq. (9). For a very large F¨oppl-von K´arm´an number vK we can approx-imate L (cid:15) as L As the temperature increases and the thermal length (cid:96) th (see Eq. (2)) becomes smaller than the membrane’s di-mensions, the elastic constants of the system are renor-malised. For ribbons with W < L this renormalisationis cut off by the width and leads to the following renor-malized elastic constants [23]: κ R ( W ) (cid:39) (cid:40) κ if W < (cid:96) th ,κ (cid:16) W (cid:96) th (cid:17) η if W > (cid:96) th , (15) Y R ( W ) (cid:39) Y if W < (cid:96) th Y (cid:16) W (cid:96) th (cid:17) − η u if W > (cid:96) th , (16)where η ≈ . − . 85 and η u ≈ . − . W and temperature-dependent stiffen-ing in the bending rigidity and softening in the Young’smodulus. Upon substituting the renormalised elastic con-stants in to the MFT, we obtain a scaling for the criti-cal 2D stress of σ c ∝ ( W /(cid:96) th ) η and a critical bucklingstrain (cid:15) c ∝ ( W /(cid:96) th ) η + η u . Because (cid:96) th ∝ T − / and using η ≈ . η + η u = 2, we seethat σ c and (cid:15) c are predicted to increase with increasingtemperature with non-trivial power laws, (cid:15) c ∼ T . and σ c ∼ T . . 4. Numerical results for finite temperature The MFT section explains how we can use the maxi-mum height h M of a compressed ribbon as an order param-eter for a buckling transition and estimate how the criticalstrain and critical stress will shift with increasing tem-perature. We now test this theoretical prediction againstmolecular dynamics simulations. We use the notation (cid:104) O (cid:105) for the average in the N V T ensemble of the observable O .It will be convenient to replace h M by the height of theribbon centre of mass h CM = N (cid:80) i h i as our order param-eter. There is, however, a subtle point to be considered. Inthe absence of an external field ( E = 0) our energy (10) isinvariant with respect to height-inversion symmetry. Thismeans that configurations with h CM = ± h have the sameprobability and would seem incompatible with the result (cid:104) h CM (cid:105) (cid:54) = 0 for (cid:15) > (cid:15) c . As with conventional magneticphase transitions, this apparent paradox is resolved by re-alising that, in the limit of large system sizes, the systemundergoes spontaneous symmetry breaking [38]. Formally,we could consider a small symmetry-breaking field E to es-tablish a preferred direction and take the double limit (cid:104) h CM (cid:105) = lim E→ lim A →∞ (cid:104) h CM (cid:105) A , (17)where A denotes the system size. Notice that if we reversedthe order of the limits (cid:104) h CM (cid:105) would always vanish. Thisis the situation in any computer simulation, where flips5 Applied field -2-1012 < h C M > ε−ε c =+0.0010 ε−ε c =+0.0007 ε−ε c =+0.0002 ε−ε c =-0.0003 ε−ε c =-0.0005 ε−ε c =-0.0008 ε−ε c =-0.0013 ε−ε c =-0.0029 -0.002 0 0.002 0.004 0.006 ε / δ ε c from stress-strain curve ε c from χ max (a) (b) Figure 3: (a) The height centre of mass (cid:104) h CM (cid:105) as a function of theout-of-plane field E for different strains relative to critical strain (cid:15) c obtained from the stress-strain curve. The slope (susceptibility) in-creases closer to the buckling transition. (b) The exponent 1 /δ asa function of (cid:15) . The critical strain obtained from the height sus-ceptibility χ (green vertical dashed line) does not coincide with the (cid:15) c obtained from the stress-strain curve (blue vertical dashed line).We see that h CM ∝ E far below the buckling transition and that itbecomes more sensitive (smaller 1 /δ ) as the system becomes closerto the transition. The 1 /δ exponent is close to 1/3 (red horizontalline) when (cid:15) is close to the value when χ is at maximum. between the up and down states are always possible aftera finite long time. The metastable dynamics for E = 0 andthe behaviour of the flipping time for a molecular dynamicssimulation will be considered in a future work [39].The previous discussion is in complete analogy to themagnetisation m of a magnetic system, where m plays therole of our height variable. In Appendix B we explore thebehaviour of the susceptibility χ = d h CM d E (cid:12)(cid:12) E =0 via simula-tions and find that these fluctuations become very large asthe buckling transition is approached from below. The definition of the broken-symmetry phase becomesdifficult for finite sizes, since the ribbon can always flip be-tween the up and down states. We can break this degen-eracy by applying an external field perpendicular to theplane. From Eq. (14) we expect steep curves of (cid:104) h CM (cid:105) E -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05 Compression ε < h C M > L / l p ~1.0, W / l th ~24 L / l p ~0.5, W / l th ~17 L / l p ~0.25, W / l th ~12 L / l p ~0.10, W / l th ~8.5 L / l p ~0.06, W / l th ~6.0 L / l p ~0.02, W / l th ~3.8 L / l p ~0.01, W / l th ~2.7 ε−ε c < h C M > π / L ε c Figure 4: Average of the squared centre-of-mass height (cid:104) h (cid:105) as afunction of compressive strain (cid:15) when E = 0. In the pre-bucklingregion (cid:104) h (cid:105) is close to zero, whereas in the post-buckling region (cid:104) h (cid:105) goes linearly with (cid:15) . The inset shows the dimensionless orderparameter (cid:104) h (cid:105) π /L (cid:15) c as a function of (cid:15) − (cid:15) c . The collapse of alldata with a slope of one, as in Eq. (19), agrees with our mean-fieldtheory. as E → 0, near the buckling transition, or equivalently d h CM d E ∝ E − / .We can test this prediction in MD simulations by chang-ing the perturbing field for compressions at a constant tem-perature. Specifically, we simulate ribbons with W /(cid:96) th ∼ . k B T / ˆ κ = 0 . 05) and apply an E up to 0.01. To savecomputing time we only simulated E > 0. In Fig. 3(a),we see that well below the buckling transition the aver-age centre-of-mass height (cid:104) h CM (cid:105) is weakly dependent onthe field. As we approach the critical buckling strain, d h CM d E (cid:12)(cid:12) E =0 becomes larger. Along the iso-strain (cid:15) = (cid:15) c where χ is at maximum, we expect (cid:104) h CM (cid:105) ∝ E /δ where δ = 3 (see Eq. 14). We can fit our data to calculate ex-ponent δ . In Fig. 3(b) we also plot the critical strain ob-tained from stress-strain curve and from the peak of theheight susceptibility χ . Interestingly, we find a propor-tionality between the critical strains obtained from stress-strain curves and critical strains obtained from the peaksof χ ; however, these two values do not coincide exactly(see Appendix C for more details). We find that 1 /δ is close to 1/3 as (cid:15) approaches (cid:15) c , where χ is maximum.We hope to investigate this proportionality in future work.Similar to magnetic-based memories, one could use the upand down buckling in a double-clamped ribbon to storeinformation, which can be controlled by compression (cid:15) ,temperature T , or perturbing out-of-plane field E . As we discussed earlier we can locate the buckling tran-sitions from stress-strain curves using data like those inFig. 2. We expect these curves will have a constant slope6lose to (cid:15) = 0, given by the Young’s modulus, and anotherslope ∼ κ R /L (cid:15) beyond the buckling point. The crossingpoint of the pre- and post-buckling curves gives the criticalbuckling load σ c and critical strain (cid:15) c (see Appendix Bfor more details).To provide a more quantitative test of the MFT, wecan use the relation h = (cid:16) L (cid:15) (cid:82) L (cid:15) / − L (cid:15) / h d x (cid:17) = h todefine a dimensionless buckling parameter at T = 0: h π L (cid:15) c ( T = 0) = (cid:15) − (cid:15) c ( T = 0) , (18)where (cid:15) c ( T = 0) = κπ Y L (cid:15) c ( T =0) . At finite temperature, weexpect the same relation to hold, with the corresponding (cid:15) c given by the renormalised constants: (cid:104) h (cid:105) π L (cid:15) c = (cid:15) − (cid:15) c . (19)Note that, at finite T , L (cid:15) c and L relax are temperature de-pendent. Fig. 4 shows (cid:104) h (cid:105) as a function of (cid:15) for differ-ent W /(cid:96) th . The linear dependence is clear. To test theMFT prediction, we subtract (cid:15) c , found previously fromthe stress-strain curve analyses, from (cid:15) . Remarkably, weindeed find a data collapse with a slope of one for (cid:15) > (cid:15) c , inaccordance with MFT and Eq. (19). At high temperatures,however, the transitions grow less sharp, presumably dueto finite-size effects. Next we examine how the elastic constants and criticalbuckling change with temperature. We plot κ R , Y R , and (cid:15) c , obtained from MD simulations, as a function of W /(cid:96) th in Fig. 5. At very low temperatures, when W /(cid:96) th (cid:28) L /(cid:96) p (cid:28) 1, these three parameters approach theirzero-temperature values. In this regime thermal fluctua-tions are weak, and thus our system behaves like a classicalribbon. In the W /(cid:96) th > κ R and softening in Y R . We test Eqs. (15)and (16) by fitting our data for W /(cid:96) th > L /(cid:96) p < Y R Y = A Y x − η u , κ R κ = A κ x η , (cid:15) c (cid:15) T =0c = A (cid:15) x η + η u , (20)where x = W /(cid:96) th .We first set the exponents to their expected values η = 0 . η u = 0 . A i to check forconsistency. The fits are excellent for the three quantities,with χ goodness-of-fit estimators per degree of freedomof χ Y / d.o.f. = 4 . / χ κ / d.o.f. = 3 . / χ (cid:15) / d.o.f. =5 . / W /(cid:96) th that can be accessed in thermalised simulations islimited. We have, however, obtained reasonable estimates -2 -1 W/l th -2 R eno r m a li z ed e l a s t i c c on s t an t s Y R / Y κ R / κ -2 -1 W / l th ε c / ε c ( ) L > l p , W > l th y~x -0.4 , χ /d.o.f=0.67 y~x , χ /d.o.f=0.50 L < l p , W > l th L < l p , W < l th y~x , χ /d.o.f=0.95 (a) (b) Figure 5: (a) Young’s modulus Y R , bending rigidity κ R and (b)critical strain (cid:15) c as a function of W /(cid:96) th . The expected theoreticalscaling of Eq. (20) is an excellent fit in the regime W > (cid:96) th , L < (cid:96) p with η = 0 . η u = 0 . of η u = 0 . η = 0 . Y R softens as the ribbon length L becomes comparableto the persistence length (cid:96) p . Very recent work by Mor-shedifard et al. also found an increase in buckling load ofsquare sheets with increasing temperature [24]. It has alsobeen shown in Ref. [40] that the critical buckling strainof Mo S sheets (described by a Stillinger-Weber poten-tial) increases with increasing temperature. To summarise,in the semi-flexible regime where L < (cid:96) p and (cid:96) th < W we find that the mechanics of thin ribbons becomes tem-perature dependent with Y R ∝ T − η u / , κ R ∝ T η/ , and (cid:15) c ∝ T ( η u + η ) / . 5. Conclusions In this letter we demonstrate that the buckling of ther-malized ribbons, when studied via molecular dynamicssimulations, can be described by a mean-field theory withrenormalized elastic constants when the ribbon length isshorter than the persistence length. We provide three in-dependent ways of locating the buckling transition. Inthe first approach we use the stress-strain curve to locatebuckling and indeed find that the buckling is delayed with7ncreasing temperature. The second approach is via heightfluctuations (Appendix B), in analogy with the study ofsusceptibility in magnetic systems. Such an increase inheight fluctuations close to the buckling transition was re-cently observed in the study of buckling of 1D colloidalsystems [41]. Lastly, we find that the height becomeshighly sensitive to an out-of-plane symmetry-breaking field E close to the transition.While the buckling transitions of thermalised nanorib-bons and phase transitions in magnetic systems seem toshare similar behaviours, the critical buckling strain issystem-size dependent ( (cid:15) c ∝ /L ), whereas the criticaltemperature of a magnetic system is typically indepen-dent of system size. Our simulations suggest regions inwhich the mean-field theory approximately holds. Theseregions are determined by the ratio between the systemsizes ( L , W ) and the relevant thermal lengths ( (cid:96) th , (cid:96) p ).In the low temperature regime ( L < (cid:96) p and (cid:96) th > W ),the classical (zero temperature) plate theory holds. Inthe intermediate (semi-flexible) regime where L < (cid:96) p and (cid:96) th < W we find that the mechanics of thin rib-bons with fixed width W can be described with a mean-field theory with temperature dependent elastic constants Y R ∝ T − η u / , κ R ∝ T η/ .Because of the softening in Y R and stiffening in κ R ,the buckling threshold increases with temperature, (cid:15) c ∝ T ( η u + η ) / . Normally η and η u are extracted from theFourier modes of height fluctuations and in-plane phonons.Here, we demonstrate that we can use an Euler buck-ling to measure these exponents directly. Current nano-fabrication techniques can create nanoribbons as thin as ∼ (cid:96) th ≈ 50 nm.It should therefore be possible to fabricate ribbons withwidth to thermal length ratio from roughly 0.01 to 100. Asimilar setup including an out-of-plane symmetry-breakingfield has been achieved experimentally [44]. The simu-lations and theory presented here provide predictions forbuckling of thermalised nanoribbons that can be tested ex-perimentally. The tunability of buckling via compression,temperature, and perturbing field could be the useful fordevelopment of mechanics-based non-volatile memories. Declaration of competing interests The authors declare that they have no known compet-ing financial interests or personal relationships that couldhave appeared to influence the work reported in this paper. Acknowledgements This research was supported in part by the NationalScience Foundation under grant no. NSF-PHY-1748958.Work by PZH and DRN was also supported through the NSF grant DMR-1608501 and via the Harvard MaterialsScience Research and Engineering Center, through NSFgrant DMR-2011754. The work of MJB was also par-tially supported by the NSF through the Materials Sci-ence and Engineering Center at UC Santa Barbara, DMR-1720256 (iSuperSeed). DY was supported by the ChanZuckerberg Biohub and received funding from the Minis-terio de Econom´ıa, Industria y Competitividad (MINECO,Spain); the Agencia Estatal de Investigaci´on (AEI, Spain)and Fondo Europeo de Desarrollo Regional (FEDER, EU)through grant no. PGC2018-094684-B-C21. DY and SBthank the KITP for hospitality during part of this project.PZH and DRN thank Abigail Plummer and Suraj Shankarfor helpful discussions. MB and DRN acknowledge help-ful conversations with Daniel Lopez. The computations inthis paper were run on the FASRC Cannon cluster sup-ported by the FAS Division of Science Research Comput-ing Group at Harvard University. Some simulations werecarried out on the Syracuse University HTC Campus Grid,which is supported by NSF-ACI-1341006. Appendix A. Numerical check of the T = 0 the-ory To check that our coarse-grained model is consistentwith the zero-temperature theory we simulated systemswith L = 100 a, W = 20 a, ˆ κ = 2 . , k/ ˆ κ = 1440 /a at T = 0. The energy is minimised using the FIRE algorithm.Recall that the connection between continuum elastic con-stants and those for a triangular lattice is κ = √ κ and Y = k √ . We plot the height amplitude h M and stress σ as a function of the compressive strain ( (cid:15) = 1 − L (cid:15) /L ) inFig. A.6. Our simulations produce a square-root scaling ofthe buckling amplitude, in agreement with the mean fieldtheory. The computed Young’s modulus, critical stress,and critical strain are within 10% of the theoretical pre-dictions [ Y simulation Y theory = 0 . , σ simulationc σ theoryc = 0 . , (cid:15) simulationc (cid:15) theoryc =0 . Appendix B. Stress-strain curve We fit data points close to (cid:15) = 0 to obtain Y R andfit data points beyond the buckling point to obtain thelinear asymptotic behaviour. We use the intersection ofthese two lines to estimate the critical buckling load σ c and critical buckling strain (cid:15) c . By plotting the scaledstress σ/σ theoryc as a function of (cid:15) , we can see that scaledcritical buckling load σ c /σ theoryc increases with increasing W /(cid:96) th (increasing T ), whereas the slope ( Y R ) decreaseswith increasing W /(cid:96) th , in accordance with the theoreticalexpectation (see Fig. B.7). Appendix B.1. The temperature-dependent critical strainfrom height susceptibility Since we are interested in the buckling response due toexternal field we study the height susceptibility defined as8 × -6 × -6 × -6 × -6 × -6 ε h M × -6 × -6 × -6 × -6 × -6 ε σ ε (x10 -6 )00.0050.01 h M (a) (b) d Figure A.6: The maximum height h M (a) and the longitudinal stress σ (b) at T = 0 as a function of the compressive strain (cid:15) = 1 − L (cid:15) /L .The inset in (a) shows the linear relationship between h and (cid:15) for (cid:15) > (cid:15) c . χ = d (cid:104) h C M (cid:105) / d E . We can directly obtain χ using heightfluctuations with Eq. (9).As discussed in the main text, the height of center ofmass h CM beyond buckling obtained from simulations offinite systems might flip after a long finite time. Thus h CM of independent runs average to zero. In simulationsof classical Ising spins it is common to take the absolutevalue of the order parameter [45], a strategy that can beadopted to our problem: χ [ | h CM | ] ∝ (cid:104) h (cid:105) − (cid:104)| h CM |(cid:105) . (B.1)Note that this quantity differs from the true susceptibil-ity (see eq 9). In MD simulations we can apply a smallsymmetry-breaking field to bias the system to buckle inone direction. Specifically, we simulated a system with W /(cid:96) th ∼ . χ [ | h CM | ] has a similar qualitative behavior andsimilar peak location to χ [ h CM ], as shown in fig. B.8. Tosave computing time we use χ [ | h CM | ] of eq. (B.1) to lo-cate the peaks. The susceptibilities ( χ [ | h CM | ]) for severaltemperatures as a function of the compressive strain (cid:15) areplotted in fig. B.9. Here and in following plots we indicatethe temperature through the ratio of the system’s width to Compressive strain ε -0.6-0.4-0.200.20.40.6 σ Compressive strain ε -20-1001020 σ / σ c t heo r y W / l th ~24 W / l th ~17 W / l th ~12 W / l th ~8.5 W / l th ~6 Figure B.7: Average stress σ as a function of compressive strain (cid:15) for a ribbon with at a temperature large enough so that W /(cid:96) th =8 . 5. The dashed lines are linear fits in the pre-buckling (small- (cid:15) )regime and in the post-buckling regime (cid:15) > (cid:15) c . The inset showsthe scaled stress σ/σ theoryc for different systems. The scaled criticalbuckling, which is proportional to the renormalised bending rigidity κ R , increases with increasing W /(cid:96) th (or temperature). In contrast,the slope ( Y R ) becomes smaller with increasing W /(cid:96) th . its thermal length, which is the appropriate scaling vari-able. We can clearly see that the buckling transition per-sists for finite T , while the position of the peaks increaseswith increasing W /(cid:96) th . This trend is consistent with ourtheoretical prediction that (cid:15) c should increase as the renor-malisation of the Young’s modulus and bending rigiditybecomes more and more important. We find a propor-tionality between the critical strains obtained from stress-strain curves and critical strains obtained from the peaksof χ ; however, these critical strain obtained from two dif-ferent approaches do not coincide exactly (see Fig. B.10). Appendix C. Variational approach Here we describe how eliminating in-plane displace-ment fields leads to a non-local stretching term in theGibbs energy. For a clamped 1D ribbon we write the effec-tive 1D Young’s modulus as Y = Y W and the bendingrigidity as κ = κW . The amount of work is − F d andthe compression distance d can be approximated as fol-lows, L (cid:15) + d (cid:39) (cid:90) L (cid:15) / − L (cid:15) / (cid:115) (cid:18) d h d x (cid:19) d x (C.1) d (cid:39) (cid:90) L (cid:15) / − L (cid:15) / (cid:18) d h d x (cid:19) d x. (C.2)9 ε < h C M > - < | h C M | > -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01Compressive strain ε < h C M > - < h C M > ε c from stress-strain curve ε c from stress-strain curve Figure B.8: (a) The susceptibility of the absolute value centre-of-mass height χ [ | h CM | ] as a function of the compressive strain (cid:15) with E = 0 and (b) the susceptibility of centre-of-mass height χ [ h CM ] asa function of the compressive strain (cid:15) with symmetry-breaking field E = 0 . W /(cid:96) th ∼ . 5. The locations of thepeaks are similar, and occur beyond the buckling strain (cid:15) c . We assume variations only in the x -direction. The Gibbsfree energy is given by G [ u x , h ] = κ (cid:90) L (cid:15) / − L (cid:15) / d x (cid:18) d h d x (cid:19) + Y (cid:90) L (cid:15) / − L (cid:15) / d x (cid:34) d u d x + 12 (cid:18) d h d x (cid:19) (cid:35) − F (cid:90) L(cid:15)/ − L(cid:15)/ d x (cid:18) d h d x (cid:19) , (C.3)where F = W σ xx . We will now focus on the middlestretching term G s controlled by Y . As is typically donein the 2D case, we focus on the vector-potential-like contri-bution, A ( x ) = (cid:0) d h d x (cid:1) , and we write the fields in Fourier -0.02 -0.01 0 0.01 0.02 Compressive strain ε χ L / l p ~1.0, W / l lth ~24 L / l p ~0.5, W / l th ~17 L / l p ~0.25, W / l th ~12 L / l p ~0.10, W / l th ~8.5 L / l p ~0.06, W / l th ~6.0 L / l p ~0.01, W / l th ~3.8 L / l p ~0.02, W / l th ~2.7 Figure B.9: The susceptibility of the absolute centre-of-mass height χ [ | h CM | ] as a function of the compressive strain (cid:15) . The position ofthe peak increases with the thermal length (cid:96) th . This is in accordancewith the theoretical expectation of a delayed buckling transition withincreasing T , due to thermal stiffening. The simulated systems coverthe temperature range of 0 . ≤ k B T/ ˆ κ ≤ . space as, d u d x = U + (cid:88) q (cid:54) =0 i q ˜ u ( q )e i qx (C.4) A = A + (cid:88) q (cid:54) =0 ˜ A ( q )e i qx , (C.5)where we have separated out the q = 0 modes. Thestretching energy G s is given by G s = Y (cid:90) L (cid:15) / − L (cid:15) / d x U + A + (cid:88) q (cid:54) =0 (cid:16) i q ˜ u ( q ) + ˜ A ( q ) (cid:17) e i qx × U + A + (cid:88) q (cid:54) =0 (cid:16) i q (cid:48) ˜ u ( q (cid:48) ) + ˜ A ( q (cid:48) ) (cid:17) e i q (cid:48) x = Y L (cid:15) ( U + A ) + Y (cid:88) q (cid:54) =0 (cid:88) q (cid:48) (cid:54) =0 (cid:90) L (cid:15) / − L (cid:15) / e i( q + q (cid:48) ) x (cid:16) i q ˜ u ( q ) + ˜ A ( q ) (cid:17) × (cid:16) i q (cid:48) ˜ u ( q (cid:48) ) + ˜ A ( q (cid:48) ) (cid:17) d x = Y L (cid:15) ( U + A ) + Y L (cid:15) (cid:88) q (cid:54) =0 | i q ˜ u ( q ) + ˜ A ( q ) | . (C.6)The stretching energy G s is clearly minimised when ˜ u ( q ) = − i ˜ A ( q ) /q . Upon imposing constant strain and the bound-10 .0005 0.001 0.0015 0.002 ε c from stress strain curve ε c f r o m χ m a x y = -4.2e-5 + 6.273*x Figure B.10: Critical strains obtained from the peaks of χ (eq. (B.1))as a function of critical strains obtained from stress-strain curves. .ary conditions we find, U = 1 L (cid:15) (cid:90) L (cid:15) / − L (cid:15) / d x d u d x = 1 L (cid:15) [ u ( L (cid:15) / − u ( − L (cid:15) / , (C.7)and similarly, A = 1 L (cid:15) (cid:90) L (cid:15) / − L (cid:15) / d x A = 1 L (cid:15) (cid:90) L (cid:15) / − L (cid:15) / (cid:18) d h d x (cid:19) d x. (C.8)Upon substituting Eq. (C.6) in the form G s = Y L (cid:15) A into Eq C.3, we obtain the free energy of Eq (10) of themain text, provided we include a contribution from thesymmetry-breaking field.We now discuss an out-of-plane field that couples tothe height. For instance, if we put uniform charges onthe ribbon and place it within a uniform electric field, thepotential energy is V ⊥ = − (cid:82) L (cid:15) / − L (cid:15) / ρ E h d x . After collectingterms, including an out-of-plane external electric field E ,we obtain the Gibbs energy G [ h, E ] = κ (cid:90) L (cid:15) / − L (cid:15) / d x (cid:18) d h d x (cid:19) + Y L (cid:15) (cid:34)(cid:90) L (cid:15) / − L (cid:15) / d x (cid:18) d h d x (cid:19) (cid:35) − F (cid:90) L (cid:15) / − L (cid:15) / d x (cid:18) d h d x (cid:19) − ρ E (cid:90) L (cid:15) / − L (cid:15) / h d x. (C.9)Close to the buckling transition we focus on the first buck-ling mode h = h M 12 (cid:104) (cid:16) πxL (cid:15) (cid:17)(cid:105) , where h M is the height amplitude, as an ansatz that insures tangential bound-ary conditions d h ( x )d x (cid:12)(cid:12)(cid:12) x ± = L (cid:15) / = 0. We then obtain Eq. (11)of the main text, G = π L (cid:15) (cid:18) κ π L (cid:15) − F (cid:19) h + π Y L (cid:15) h − ρL (cid:15) E h M . (C.10)It is helpful to write the above equation in terms of newparameters a, b, (cid:15) c G = a ( (cid:15) − (cid:15) c ) h + bh − ρL (cid:15) E h M , (C.11)where a = Y π L (cid:15) , b = π Y L (cid:15) , (cid:15) c = π κ Y L (cid:15) . Upon min-imising the Gibbs free energy by setting d G d h M (cid:12)(cid:12)(cid:12) E =0 = 0, wefind h M = (cid:40) , if (cid:15) < (cid:15) c ± (cid:112) a b ( (cid:15) − (cid:15) c ) = L (cid:15) c π (cid:113) (cid:15) − κπ Y L (cid:15) c , if (cid:15) > (cid:15) c (C.12) Appendix C.1. Susceptibility To obtain the susceptibility at zero external field wefirst solve d G d h M = 0, which leads to0 =2 ah M ( (cid:15) c − (cid:15) ) + 4 bh − ρL (cid:15) E E = 2 ρL (cid:15) [4 bh + 2 a ( (cid:15) c − (cid:15) ) h M ] . (C.14)We can now calculate the susceptibility and use ρ = Q/L and Y = Y W to obtain ∂h M ∂ E (cid:12)(cid:12)(cid:12)(cid:12) E =0 = (cid:40) QY π L (cid:15) c W ( (cid:15) c − (cid:15) ) − if (cid:15) < (cid:15) c Q Y π L (cid:15) c W ( (cid:15) − (cid:15) c ) − if (cid:15) > (cid:15) c (C.15)For a ribbon in a gravitational field simply replace E = g and Q = m , where m is the total mass. References [1] B. Grosso, E. Mele, Graphene gets bent, Physics Today 73(2020) 46. doi:10.1063/PT.3.4569 .[2] M. K. Blees, A. W. Barnard, P. A. Rose, S. P. Roberts, K. L.McGill, P. Y. Huang, A. R. Ruyack, J. W. Kevek, B. Kobrin,D. A. Muller, P. L. McEuen, Graphene kirigami, Nature 524(2015) 204–207. doi:10.1038/nature14588 .[3] T. C. Shyu, P. F. Damasceno, P. M. Dodd, A. Lamoureux,L. Xu, M. Shlian, M. Shtein, S. C. Glotzer, N. A. Kotov, Akirigami approach to engineering elasticity in nanocompositesthrough patterned defects, Nature materials 14 (8) (2015) 785. doi:10.1038/nmat4327 .[4] P. Z. Hanakata, Z. Qi, D. K. Campbell, H. S. Park, Highlystretchable MoS kirigami, Nanoscale 8 (1) (2016) 458–463. doi:10.1039/C5NR06431G . 5] Y. Tang, J. Yin, Design of cut unit geometry in hierarchicalkirigami-based auxetic metamaterials for high stretchability andcompressibility, Extreme Mechanics Letters 12 (2017) 77–85. doi:10.1016/j.eml.2016.07.005 .[6] A. Rafsanjani, K. Bertoldi, Buckling-induced kirigami, Phys.Rev. Lett. 118 (2017) 084301. doi:10.1103/PhysRevLett.118.084301 .[7] P. Z. Hanakata, E. D. Cubuk, D. K. Campbell, H. S. Park,Accelerated search and design of stretchable graphene kirigamiusing machine learning, Physical review letters 121 (25) (2018)255304. doi:10.1103/PhysRevLett.121.255304 .[8] M. Moshe, E. Esposito, S. Shankar, B. Bircan, I. Cohen, D. R.Nelson, M. J. Bowick, Kirigami mechanics as stress relief byelastic charges, Phys. Rev. Lett. 122 (2019) 048001. doi:10.1103/PhysRevLett.122.048001 .[9] Y. Yang, M. A. Dias, D. P. Holmes, Multistable kirigami fortunable architected materials, Physical Review Materials 2 (11)(2018) 110601. doi:10.1103/PhysRevMaterials.2.110601 .[10] M. A. Dias, M. P. McCarron, D. Rayneau-Kirkhope, P. Z.Hanakata, D. K. Campbell, H. S. Park, D. P. Holmes, Kirigamiactuators, Soft matter 13 (48) (2017) 9087–9092. doi:10.1039/C7SM01693J .[11] A. Rafsanjani, Y. Zhang, B. Liu, S. M. Rubinstein, K. Bertoldi,Kirigami skins make a simple soft actuator crawl, Sci-ence Robotics 3 (2018) eaar7555. doi:10.1126/scirobotics.aar7555 .[12] Y. Morikawa, S. Yamagiwa, H. Sawahata, R. Numano,K. Koida, T. Kawano, Donut-shaped stretchable kirigami: En-abling electronics to integrate with the deformable muscle,Advanced Healthcare Materials 8 (23) (2019) 1900939. doi:10.1002/adhm.201900939 .[13] M. Katsnelson, Graphene: carbon in two dimensions, Cam-bridge University Press, 2012.[14] D. Nelson, T. Piran, S. Weinberg, Statistical Mechanics of Mem-branes and Surfaces, 2nd Edition, World Scientific, Singapore,2004.[15] D. Nelson, L. Peliti, Fluctuations in membranes with crystallineand hexatic order, J. Phys. France 48 (1987) 1085–1092. doi:10.1051/jphys:019870048070108500 .[16] J. A. Aronovitz, T. C. Lubensky, Fluctuations of solid mem-branes, Phys. Rev. Lett. 60 (1988) 2634–2637. doi:10.1103/PhysRevLett.60.2634 .[17] E. Guitter, F. David, S. Leibler, L. Peliti, Thermodynamical be-havior of polymerized membranes, Journal de Physique 50(14)(1989) 1787–1819. doi:10.1051/jphys:0198900500140178700 .[18] P. Le Doussal, L. Radzihovsky, Self-consistent theory of poly-merized membranes, Phys. Rev. Lett. 69 (1992) 1209–1212. doi:10.1103/PhysRevLett.69.1209 .[19] Z. Zhang, H. T. Davis, D. M. Kroll, Scaling behavior of self-avoiding tethered vesicles, Phys. Rev. E 48 (1993) R651–R654. doi:10.1103/PhysRevE.48.R651 .[20] M. J. Bowick, S. M. Catterall, M. Falcioni, G. Thorleifsson,K. N. Anagnostopoulos, The flat phase of crystalline mem-branes, J. Phys. I France 6 (1996) 1321–1345. doi:10.1051/jp1:1996139 .[21] A. Koˇsmrlj, D. R. Nelson, Mechanical properties of warpedmembranes, Physical Review E 88 (1) (2013) 012136.[22] K. N. Kudin, G. E. Scuseria, B. I. Yakobson, C F, BN, andC nanoshell elasticity from ab initio computations, PhysicalReview B 64 (23) (2001) 235406. doi:10.1103/PhysRevB.64.235406 .[23] A. Koˇsmrlj, D. R. Nelson, Response of thermalized ribbons topulling and bending, Phys. Rev. B 93 (2016) 125431. doi:10.1103/PhysRevB.93.125431 .[24] A. Morshedifard, M. Ruiz-Garcia, M. J. A. Qomi, A. Kosm-rlj, Buckling of thermalized elastic sheets, arXiv preprint arXiv:2005.05949 .[25] A. Stukowski, Visualization and analysis of atomistic simulationdata with OVITO-the open visualization tool, Modelling andSimulation in Materials Science and Engineering 18. doi:10.1088/0965-0393/18/1/015012 . [26] M. J. Bowick, A. Travesset, The statistical mechanics ofmembranes, Phys. Rep. 344 (2001) 255–308. doi:10.1016/S0370-1573(00)00128-9 .[27] H. S. Seung, D. R. Nelson, Defects in flexible membranes withcrystalline order, Phys. Rev. A 38 (1988) 1005–1018. doi:10.1103/PhysRevA.38.1005 .[28] M. J. Bowick, A. Koˇsmrlj, D. R. Nelson, R. Sknepnek, Non-Hookean statistical mechanics of clamped graphene ribbons,Physical Review B 95 (10) (2017) 104109. doi:10.1103/PhysRevB.95.104109 .[29] D. Yllanes, S. S. Bhabesh, D. R. Nelson, M. J. Bowick, Thermalcrumpling of perforated two-dimensional sheets, Nat. Comm. 8(2017) 1381. doi:10.1038/s41467-017-01551-y .[30] D. Yllanes, D. R. Nelson, M. J. Bowick, Folding pathwaysto crumpling in thermalized elastic frames, Phys. Rev. E 100(2019) 042112. doi:10.1103/PhysRevE.100.042112 .[31] J. A. Anderson, J. Glaser, S. C. Glotzer, HOOMD-blue: APython package for high-performance molecular dynamics andhard particle monte carlo simulations, Computational MaterialsScience 173 (2020) 109363. doi:10.1016/j.commatsci.2019.109363 .[32] A. P. Young, Everything you wanted to know about Data Anal-ysis and Fitting but were afraid to ask, Springer, Berlin, 2015. arXiv:1210.3781 .[33] L. D. Landau, E. M. Lifshitz, Theory of Elasticity, 3rd Edition,Butterworth-Heinemann,, Singapore, 1999.[34] R. Lifshitz, M. Cross, Nonlinear dynamics of nanomechanicaland micromechanical resonators, Reviews of nonlinear dynamicsand complexity 1 (2008) 1–52. doi:10.1002/9783527626359 .[35] J.-P. Kownacki, D. Mouhanna, Crumpling transition and flatphase of polymerized phantom membranes, Phys. Rev. E 79(2009) 040101. doi:10.1103/PhysRevE.79.040101 .[36] J. Los, M. I. Katsnelson, O. Yazyev, K. Zakharchenko, A. Fa-solino, Scaling properties of flexible membranes from atomisticsimulations: application to graphene, Physical Review B 80 (12)(2009) 121405. doi:10.1103/PhysRevB.80.121405 .[37] R. Rold´an, A. Fasolino, K. V. Zakharchenko, M. I. Katsnel-son, Suppression of anharmonicities in crystalline membranesby external strain, Physical Review B 83 (17) (2011) 174104. doi:10.1103/PhysRevB.83.174104 .[38] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena,4th Edition, Clarendon Press, Oxford, 2005.[39] S. S. Bhabesh, P. Z. Hanakata, M. J. Bowick, D. R. Nelson,D. Yllanes, to be published.[40] J. W. Jiang, The buckling of single-layer MoS under uniaxialcompression, Nanotechnology 25 (35) (2014) 355402. doi:10.1088/0957-4484/25/35/355402 .[41] S. Stuij, J. M. van Doorn, T. Kodger, J. Sprakel, C. Coulais,P. Schall, Stochastic buckling of self-assembled colloidal struc-tures, Physical Review Research 1 (2019) 023033. doi:10.1103/PhysRevResearch.1.023033 .[42] P. Masih Das, G. Danda, A. Cupo, W. M. Parkin, L. Liang,N. Kharche, X. Ling, S. Huang, M. S. Dresselhaus, V. Meu-nier, et al., Controlled sculpture of black phosphorus nanorib-bons, ACS nano 10 (2016) 5687–5695. doi:10.1021/acsnano.6b02435 .[43] I. R. Storch, R. De Alba, V. P. Adiga, T. Abhilash, R. A. Bar-ton, H. G. Craighead, J. M. Parpia, P. L. McEuen, Young’smodulus and thermal expansion of tensioned graphene mem-branes, Physical Review B 98 (8) (2018) 085408. doi:10.1103/PhysRevB.98.085408 .[44] N. Lindahl, D. Midtvedt, J. Svensson, O. A. Nerushev, N. Lind-vall, A. Isacsson, E. E. Campbell, Determination of the bend-ing rigidity of graphene via electrostatic actuation of buckledmembranes, Nano letters 12 (2012) 3526–3531. doi:10.1021/nl301080v .[45] A. W. Sandvik, Computational studies of quantum spin sys-tems, in: AIP Conference Proceedings, Vol. 1297, AmericanInstitute of Physics, 2010, pp. 135–338..[45] A. W. Sandvik, Computational studies of quantum spin sys-tems, in: AIP Conference Proceedings, Vol. 1297, AmericanInstitute of Physics, 2010, pp. 135–338.