Folding pathways to crumpling in thermalized elastic frames

TL;DR

This paper investigates how thermal fluctuations induce folding and crumpling in elastic frames with holes, revealing a sequence of origami-like folds as temperature increases, using molecular dynamics simulations.

Contribution

It introduces a detailed simulation study of crumpling pathways in elastic frames with holes, highlighting the role of origami-like folding at different scales.

Findings

Folding occurs in a sequence of origami-like steps with increasing temperature.

Normal-normal correlation functions quantify the number of folds at various temperatures.

The system approaches a crumpled phase through a hierarchical folding process.

Abstract

The mechanical properties of thermally excited two-dimensional crystalline membranes can depend dramatically on their geometry and topology. A particularly relevant example is the effect on the crumpling transition of holes in the membrane. Here we use molecular dynamics simulations to study the case of elastic frames (sheets with a single large hole in the center) and find that the system approaches the crumpled phase through a sequence of origami-like folds at decreasing length scales when temperature is increased. We use normal-normal correlation functions to quantify the temperature-dependent number of folds.

Equations

$$\mathcal{H}={\mathcal{H}}_{\text{stretch}}+{\mathcal{H}}_{\text{bend}},$$

$${\mathcal{H}}_{\text{stretch}}$$

$${\mathcal{H}}_{\text{bend}}$$

$${\mathcal{R}}_{\text{g}}^2$$

$$C=\frac{1}{N}\bigl{(}\langle\mathcal{H}^{2}\rangle-\langle\mathcal{H}\rangle^{2}\bigr{)}\ .$$

$$⟨{\hat{\mathbf{n}}}_O\cdot {\hat{\mathbf{n}}}_X{⟩}_L=1-\frac{kT}{2\pi {\kappa }_0}\left[{\eta }^{-1}+\log \left(\frac{{\ell }_{\text{th}}}{a}\right)+\mathcal{C}\frac{kT}{{\kappa }_0}{\left(\frac{{\ell }_{\text{th}}}{L}\right)}^{\eta }\right].$$

$$C(r)=\sum _s⟨\hat{\mathbf{n}}(s)\cdot \hat{\mathbf{n}}(s+r)⟩.$$

$${\ell }_{\text{p}}$$

$$F(k)=\biggl{\langle}\biggl{|}\sum_{r}C(r)\mathrm{e}^{\mathrm{i}kr}\biggr{|}^{2}\biggr{\rangle}$$

Full text

I Introduction

Studies of two-dimensional materials such as graphene Katsnelson (2012) have stimulated renewed exploration of the statistical mechanics of elastic membranes. Thermalized elastic membranes Nelson (2002); Nelson et al. (2004) have a low-temperature extended (“flat”) phase with long-range order in the surface normals and a scale-dependent bending rigidity ($\kappa$) and Young’s modulus ($Y$) Nelson and Peliti (1987); Aronovitz and Lubensky (1988); Guitter et al. (1989); Le Doussal and Radzihovsky (1992); Zhang et al. (1993); Bowick et al. (1996). Without strong distant self-avoidance these materials are also believed to undergo a crumpling transition when the microscopic bending rigidity is comparable to the scale of thermal fluctuations Bowick and Travesset (2001); Nelson et al. (2004); Kantor and Nelson (1987); Bowick et al. (1996); Cuerno et al. (2016). The crumpling transition, however, has never been convincingly observed in a physical system. The challenge is to find a thin material that exhibits both the flat phase and the crumpled phase. One might expect that soft flexible systems, such as pure amphiphilic bilayers Shum et al. (2008), would then be the natural setting to observe the transition from a flat phase to a crumpled phase. But here length scales become important. The flat phase itself is stabilized by the strong thermal stiffening of the bending rigidity resulting from soft flexural phonons that give rise to isotropic thermal corrugations. These fluctuations set in only above the thermal length scale, $\ell_{\text{th}}$, beyond which the renormalization of the bending rigidity becomes comparable to the microscopic (bare) bending rigidity $\kappa_{0}$. This scale, which follows from the renormalization group flow of the bending rigidity, is given by $\ell_{\text{th}}\sim\kappa_{0}/\sqrt{Y_{0}k_{B}T}$, where $Y_{0}$ is the bare Young’s modulus. Even for the diblock copolymers of Ref. Shum et al. (2008), however, this length is relatively large, $\ell_{\text{th}}\sim 300$ nm at room temperature. It is therefore difficult to fabricate soft systems large enough to harness these thermal effects. Graphene, on the other hand, has a thermal length scale of the order of nanometers at room temperature, making it possible to observe a crinkled flat phase Nicholl et al. (2015); Blees et al. (2015). The difficulty now is accessing the crumpling transition which, considering that for graphene $\kappa_{0}\approx 1.25$ eV Fasolino et al. (2007), would occur at temperatures of order $50\,000$ K!

In Yllanes et al. (2017) it was shown that the crumpling transition can be dramatically lowered by either perforating the membranes with a periodic array of holes or by excising a single large hole to form a thin frame. In both cases there are two effects which lower the crumpling temperature. There is less material to bend and there are fewer multiples of the thermal length scale contributing to the stiffening of the renormalized bending rigidity. For perforated membranes the growth is cut off by the mean spacing between holes and for thin frames by the narrow frame width. Here we explore crumpling for a thin hexagonal frame so as to preserve as many of the symmetries of the triangular lattice as possible. We find that crumpling is indeed enhanced, but also uncover a striking pathway to the crumpled phase. A well-known elementary bending mode is an origami-like fold along a line Kantor and Jaric (1990); Di Francesco and Guitter (1994); Santangelo (2017). The energy cost for such a fold is proportional to the length of the fold line rather than the typical area cost for introducing a complex network of crinkles designed to crush a macroscopic piece of membrane Vliegenthart and Gompper (2006); Lahini et al. (2017).

With our molecular dynamics simulations we find that crumpling occurs via a sequence of these soft origami folds. The successive stages of crumpling may be characterized by the distribution of the fold lines. The patterns observed correspond to the fold lines for the exactly solvable model of planar folding first studied in Di Francesco and Guitter (1994).

The rest of the paper is organized as follows: Section II presents our coarse-grained crystalline model and describes our simulations. In Section III we define an order parameter sensitive to the degree of crumpling based on normal-normal correlations and describe the folding process qualitatively. This description is made more quantitative in Section IV with the use of the full normal-normal correlation function. Finally, Section V examines the effect of changing the size of the hole and Section VI summarizes our conclusions. We include an appendix that considers the normal-normal correlation function in Fourier space.

II Model

We consider a two-dimensional sheet of elastic material with a large hole in the center (a frame), which we model as a crystalline membrane Bowick and Travesset (2001). In this representation, the sheet is discretized with a tiling of equilateral triangles of side $a=1$, which defines a lattice of unbreakable but elastic bonds. In order both to facilitate triangulation and to approach a circular symmetry we use a hexagonal geometry with inner radius $R_{1}$ and outer radius $R_{2}=R_{1}+W$ (Figure 1). This system can be regarded as six ribbons of length $R=(R_{1}+R_{2})/2$ stitched together to form a hexagon. On this triangular lattice we define a standard coarse-grained elastic Hamiltonian Seung and Nelson (1988):

[TABLE]

where

[TABLE]

The sum in the stretching term is over all pairs of lattice nearest neighbors $\langle i,j\rangle$. The bending energy, on the other hand, depends on the angles between the normals $\boldsymbol{n}_{\alpha}$, $\boldsymbol{n}_{\beta}$ of all pairs of triangles $(\alpha,\beta)$ that share a side. Each of these pairs of triangles defines a dihedral, which is how we will refer to each term in the bending Hamiltonian in the rest of the paper.

The coupling constants $\epsilon$ and $\tilde{\kappa}$ are directly related to the Young’s modulus and bare bending rigidity of continuum elastic theory ($Y_{0}=2\epsilon/\sqrt{3}$, $\kappa_{0}=\sqrt{3}\tilde{\kappa}/2$ Seung and Nelson (1988)). As temperature increases, the frame transitions from a thermally excited flat phase, with long-range order in the normals, to a crumpled phase (see Figure 1). This crumpling transition has already been thoroughly studied for unperforated sheets Bowick and Travesset (2001); Nelson et al. (2004); Kantor and Nelson (1987); Le Doussal and Radzihovsky (1992); Bowick et al. (1996); Cuerno et al. (2016) as well as for sheets with perforations Yllanes et al. (2017). In the following, therefore, we do not focus on the crumpling transition itself, but instead study an intermediate regime, where we show that crumpling is realized through a series of folding pathways.

Aside from an overall constant factor that defines the energy units, what differentiates one elastic material from another is the ratio $\epsilon/\tilde{\kappa}$. Following Bowick et al. (2017); Yllanes et al. (2017) we set $\epsilon=1440\tilde{\kappa}/a^{2}$, which approximates graphene values without making the system too rigid and hard to thermalize. In practice, this choice makes very little difference in the final results, since the crumpling point depends only logarithmically on the Young’s modulus (see, e.g., Košmrlj and Nelson (2016) and compare the crumpling temperatures in Cuerno et al. (2016) and Yllanes et al. (2017)).

We have carried out Molecular Dynamics simulations of model (1) for different system sizes and a wide temperature range (working in an $NVT$ ensemble with a standard Nosé-Hoover thermostat Nosé (1984); Hoover (1985)). Our simulations were implemented in the HOOMD-blue package Glaser et al. (2015); Anderson et al. (2008) and run on Tesla GPUs. Following Yllanes et al. (2017), we use a timestep of $\Delta t=0.0025\tau$, where $\tau=\sqrt{ma^{2}/kT}$ is the Lennard-Jones unit of time (expressed in terms of the particle mass $m$) and we use natural units with $a=m=1$. All energies are measured in units of $kT$. For each system size and temperature we run for $10^{9}$ MD steps and use a jackknife procedure Amit and Martin-Mayor (2005) to estimate statistical errors. For our larger frames with $R_{1}=70,R_{2}=82$ these runs take from 60 to 70 hours on a Tesla K40m GPU.

III The crumpling order parameter

The standard method to study the crumpling transition itself is to consider the radius of gyration of the frame

[TABLE]

In this equation and in the following, $\langle\cdot\rangle$ will denote a thermal average. At low temperatures, $\mathcal{R}_{\text{g}}^{2}\sim R_{2}^{2}$, while in the crumpled phase $\mathcal{R}_{\text{g}}^{2}\sim\log(R_{2}/a)$. Panels (a) and (b) of Figure 2 show $\mathcal{R}_{\text{g}}$ as a function of temperature for two frame sizes.

A very clear signal of the phase transition can also be obtained by plotting the specific heat $C$ of our elastic frames a function of temperature, which can be computed with a fluctuation-response formula Yllanes et al. (2017):

[TABLE]

This quantity has a very sharp peak at the transition point of $kT_{\text{c}}\approx 4.0\tilde{\kappa}$ for both values of $R_{1}$ studied (see Figures 2c and d).

Finite-size scaling studies of the crumpling transition have already been done in considerable detail both for pristine Cuerno et al. (2016) and perforated Yllanes et al. (2017) sheets. Here, we are interested in following the detailed geometry of the frame as the transition is approached. To this end, we have found that a system size of $R_{1}=30$ and $R_{2}=42$ provides a good compromise: The frame is large enough for finite-size effects to be negligible (compare the position of the peaks for $R_{1}=30$ and $R_{1}=70$ in Figure 2) and to provide good statistics, yet small enough for us to follow individual folds and creases in snapshots. In what follows we shall always consider this particular system size.

Even though $\mathcal{R}_{\text{g}}$ is a useful proxy for a crumpling order parameter, it is illuminating to study normal-normal correlations directly. Because our frames are made up of six ribbons of length $R=(R_{1}+R_{2})/2$ stitched together to form a hexagon (see Fig. 1) we first consider a long elastic ribbon at a low temperature. Then the correlation between the normal at the beginning ($O$) and end ($X$) of the ribbon is Košmrlj and Nelson (2016)

[TABLE]

In this equation, $\eta\approx 0.8$ is a critical exponent, $\mathcal{C}$ is a positive constant of order unity and $\ell_{\text{th}}^{2}=2\sqrt{3}\pi^{3}\tilde{\kappa}^{2}/3kT\epsilon$ is the thermal length scale. For sufficiently small $kT$, this correlation is positive even for $r_{OX}=L\to\infty$. Conversely, we could in principle approximate the crumpling transition as the point where $\langle\boldsymbol{n}_{O}\cdot\boldsymbol{n}_{X}\rangle\approx 0$ for large $r_{OX}=L$.

This test is easily done by considering our frame to be a ribbon of width $R_{2}-R_{1}=W$, with periodic boundary conditions in the azimuthal direction. If we then consider a central hexagon of radius $R=(R_{1}+R_{2})/2$, we can identify $\langle\boldsymbol{n}_{O}\cdot\boldsymbol{n}_{X}\rangle$ with $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{D}\rangle$, the correlation between diagonally opposed vertices (see Figure 3) 111We actually consider the average of $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{D}\rangle$, $\langle\boldsymbol{n}_{B}\cdot\boldsymbol{n}_{E}\rangle$ and $\langle\boldsymbol{n}_{C}\cdot\boldsymbol{n}_{F}\rangle$..

If we plot $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{D}\rangle$ as a function of temperature, however, we find that it goes to zero well before the crumpling temperature of $kT_{\text{c}}\approx 4.0\tilde{\kappa}$, although the data is rather noisy. We can obtain a better plot by studying instead the correlation between opposing vertices along the same side of the hexagon ($\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{B}\rangle$ and five symmetric equivalents). This quantity goes to zero at a higher temperature than $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{D}\rangle$, but still well before the actual crumpling transition. Furthermore, for a finite temperature range, $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{B}\rangle$ is negative, indicating anticorrelation between normals.

Let us consider the meaning of this negative correlation. The easiest way to achieve $\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{B}=-1$ is to create a single fold between the two points. Of course, our frame is fluctuating thermally. A measurable negative ensemble averaged $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{B}\rangle$, however, suggests that typical instantaneous configurations of the frame might have well defined origami folds, separated by a distance comparable to $|\boldsymbol{r}_{A}-\boldsymbol{r}_{B}|$, but whose location varies with time. Upon recalling that the correlation for the even more distant points such as $A$ and $D$ goes to zero at an even lower temperature [Figure 3-(c)], we can hypothesize that there exists an intermediate regime before crumpling where the system concentrates its bending energy on a small number of well-defined folds. We expect that the distance between single folds is a function of temperature, which becomes quite small as $T$ approaches the crumpling transition.

We can test this idea of a temperature-dependent folding scale by plotting a heat map of bending energy on the frame (Figure 4). For low temperatures, there is hardly any bending energy in the system. As $T$ is increased, we find that almost all the bending energy is concentrated in one or two radial segments, indicating well-defined creases. The number of such folds increases with temperature until, at the crumpling transition $kT_{\text{c}}\approx 4.0\tilde{\kappa}$, the whole system is crumpled (i.e., folded at a microscopic scale).

IV Normal-normal correlation function

We can make the analysis more quantitative by studying $\hat{\boldsymbol{n}}(r)$, the normal to the sheet along the central hexagon of radius $R=(R_{1}+R_{2})/2$, as a function of the arc length. Then we can compute

[TABLE]

This function gives us more complete information than the correlation between particular pairs of points of the previous section and is also more precise, since we can average over the whole curve of annular midpoints around the hexagon.

The result for $R_{1}=30$, $R_{2}=42$ is plotted in Figure 5, where several regimes can be identified. At the lowest temperatures, the system is deep in the flat phase and $C(r)$ has a positive asymptote as $r\to\infty$. At the highest temperatures, in the crumpled phase, $C(r)$ is trivial: it is negative for $r=2$ (opposite vertices of one dihedral) and zero beyond that since the system is folded at a microscopic scale. Qualitatively similar behavior was observed, e.g. Kantor and Nelson (1987), for tethered surfaces with the geometry of a parallelogram. Our interest here is the existence of an intermediate regime where $C(r)$ has negative values for a finite $r$ window, which shifts and shrinks as $T$ approaches $T_{\text{c}}$. Just prior to this anomalous behavior, we find that $C(r)$ is well described by an exponential decay $C(r)\simeq\exp(-r/\ell_{\text{p}})$, where $\ell_{\text{p}}$ is a persistence length (see, e.g. Hoover (1985); Košmrlj and Nelson (2016); Yllanes et al. (2017)). See the curves with $kT/\tilde{\kappa}\leq 2.00$ in Figure 5.

Since the $C(r)$ are very smooth functions and can be measured to high precision, we can easily find a minimum for each one in the folding regime. This minimum then indicates the average distance between folds (or, alternatively, their average number) for each temperature (Figure 6).

V Role of the system geometry

As we have mentioned, the $C(r)$ function has been extensively studied for unperforated membranes. In these cases, there is no folding regime and the system suddenly switches between a flat phase where $C(r)>0$ for all $r$ and an extremely compact crumpled phase where $C(r=2)<0$ and $C(r>2)=0$. This behavior arises because in unperforated membranes the system presumably cannot find a privileged direction along which to fold in order to minimize its bending energy.

We have tested these ideas by redoing our simulations for thicker frames (smaller holes). In particular, keeping $R_{2}=42$, we have run simulations for $R_{1}=10,20$ ($W=22,32$, respectively, compared with the $W=12$ studied above). In this analysis, we have found that for $W=22$ a folding regime can still be identified (though it spans a much narrower temperature interval than for $W=12$), while for $W=32$ the behavior is already qualitatively the same as for sheets without a hole 222And quantitatively very similar, with $kT_{\text{c}}$ approaching the value without a hole, although, unlike for our thin frames, in this case finite-size effects in $kT_{\text{c}}$ are noticeable..

This study is summarized in Figures 7 to 9. First, Figure 7 plots the radius of gyration and $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{B}\rangle$ for the three frame widths. As discussed in the previous sections, the interval between the temperature at which $\langle\boldsymbol{n}_{A}\cdot\boldsymbol{n}_{B}\rangle=0$ (appearance of very large folds) and $T_{\text{c}}$ corresponds to the folding regime. This is very wide for $W=12$, noticeably narrower for $W=22$ and has essentially disappeared for $W=32$.

Moreover, if we consider (Figure 8) the same bending energy heat map we used for our original frames, we see that even at a temperature just $5\%$ below $T_{\text{c}}$ no folds can be identified for our $W=32$ frame. The would-be folding regime has already shrunk beyond our ability to detect it. Finally, the picture is confirmed by considering $C(r)$ again in Figure 9, where the correlation functions suddenly change from having a positive asymptotic value (flat phase) to collapsing completely. This is the same pattern observed for membranes without holes in Ref. Kantor and Nelson (1987).

VI Conclusions

We have studied the equilibrium configurations of thermalized elastic frames as temperature is increased and found that, as the crumpling transition is approached, most of the bending energy is concentrated on a growing number of origami-like folds. The scale of these folds is a well-defined function of temperature.

The folding pathway explored here neglects distant self-avoidance. In Ref. Yllanes et al. (2017) we argued that the crumpling transition studied here for phantom frames with large holes would persist even in the presence of self-avoidance in an appropriate thermodynamic limit. The non-trivial width-dependent scaling of the thermally generated renormalized persistence length associated with our hexagonal frames

[TABLE]

suggests that the appropriate limit is $R,W\to\infty$ with fixed $W(W/\ell_{\text{th}})^{\eta}/R$ for frames with width $W$ and edge length $R$. Thus, we expect that a sharp transition survives for hexagonal frames, where both a crumpled and a flat phase would survive in a polymer-like large-size limit. In addition, we expect that the origami-like folding pathway uncovered here will not be significantly affected by distant self-avoidance for most of the run-up to the crumpling transition from low temperatures. Indeed, few self-intersections are evident in the partially crumpled images shown in Figs. 1 and 4. How self-avoidance affects the folding at the transition itself when $kT\approx 4.0\tilde{\kappa}$ is an interesting subject for future investigation.

Some insight into the effect of self-avoidance on folding patterns without holes follows from the study of Vliegenthart and Gompper of forced crumpling of self-avoiding elastic sheets at zero temperature Vliegenthart and Gompper (2006). Here, a comparison of gradual forced crumpling with and without self-avoidance eventually reveals a radial bias, just as we find here with a single large hole. The patterns of wrinkles in this computer experiment are qualitatively similar for both phantom and self-avoiding membranes.

Acknowledgements.

This research was supported by the NSF through the DMREF grant DMR-1435794 and DGE-1068780, as well as by the Syracuse University Soft Matter Program. Work by DRN was also supported through the NSF DMREF program, via NSF grant DMR-1435999 and via the Harvard Materials Science Research and Engineering Center, through NSF grant DMR-1420570. The research of MJB was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. DY thanks the KITP for hospitality during part of this project and acknowledges funding by Ministerio de Economía, Industria y Competitividad (MINECO) (Spain) through grants no. FIS2015-65078-C2 and PGC2018-094684-B-C21 (also partly funded by the EU through the FEDER program) and the resources and assistance provided by BIFI-ZCAM (Universidad de Zaragoza), where we carried out most of our simulations on the Cierzo supercomputer.

Appendix A The normal-normal correlation in Fourier space

In the main part of this paper we have worked with the correlation function in real (position) space. However, its Fourier transform

[TABLE]

provides an interesting complementary picture. Now the folding regime leads to a maximum in $F(k)$ that shifts to smaller wavenumbers as $T$ increases (Figure 10). In addition, considering $F(k)$ as a function of temperature offers an alternative way of finding the crumpling point (Figure 11).

Appendix B The crumpling process

In the main text of the paper, we considered the folding pathways as the temperature was increased up to the crumpling transition, which we studied through the analysis of configurations extracted from the equilibrium evolution of the system. We can obtain a complementary picture by studying the non-equilibrium evolution of a system at $T_{\text{c}}$ which started from a completely flat configuration. We have done this in Figures 13 and 13 for our systems with $R_{2}=42$ and $R_{1}=10,30$ (each at their corresponding $T_{\text{c}}$). For the latter, individual folds quickly appear at several points along the frame and serve to nucleate crumpled regions. It takes longer for the thicker frame to crumple, since the system cannot create linear folds that easily. In both cases, once the system has crumpled completely, the microscopic fold patterns correspond to those of the planar folding problem, as in Di Francesco and Guitter (1994).