{'fraction_non_alphanumeric': 0.03904109589041096, 'fraction_numerical': 0.01155449672424062, 'mean_word_length': 4.360095770151636, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of rippling in incommensurate one-dimensional coupled chains

Graphene and other recently developed 2D materials exhibit exceptionally strong in-plane stiffness. Relaxation of few-layer structures, either free-standing or on slightly mismatched substrates occurs mostly through out-of-plane bending and the creation of large-scale ripples. In this work, we present a novel double chain model, where we allow relaxation to occur by bending of the incommensurate coupled system of chains. As we will see, this model can be seen as a new application of the well-known Frenkel-Kontorova model for a one-dimensional atomic chain lying in a periodic potential. We focus in particular on modeling and analyzing ripples occurring in ground state configurations, as well as their numerical simulation.

Introduction. Two-dimensional layered crystals such as graphene, nature's thinnest elastic material, exhibit exceptionally strong in-plane stiffness but are also highly flexible with extremely low flexural rigidity compared to conventional membranes. Mono-and few-layer graphene structures display spontaneous creation of long-range ripples which have been experimentally observed either in suspended conditions [14] or supported on slightly mismatched substrates [5,21,24]. Ripples in graphene is a multiscale phenomenon that is expected to influence its mechanical and electronic properties [17], and inducing periodic ripples in a controlled fashion may open new perspectives, e.g., for graphene-based devices [3,15].
Various mechanisms have been described as responsible for the appearance of ripples in two-dimensional materials. Spontaneous ripples in suspended graphene monolayers have been attributed to thermal fluctuations [8] or stress [3]. Graphene monolayers supported by an almost-commensurate substrate relax to form a periodic vertical corrugation pattern [5,18,24,23]. Sharp out-of-plane folds called ripplocations have been predicted in van der Waals (vdW) homostructures such as MoS 2 multilayers, similar to line defects [13]. In a ripplocation, an additional local line of atoms is inserted in one layer leading to the formation of an out-of-plane wrinkle, while the other layer remains flat.
In this paper, we will study a different mechanism yet, where the spontaneous atomic-scale relaxation of free-standing systems of incommensurate vdW bilayers leads to a simultaneous long-range rippling of the bilayer system [19]. Note that vdW multi-layer structures tend to form naturally incommensurate stackings, either due to a relative rotation of the crystalline orientations or to a natural mismatch between the respective lattice constants. To model mathematically this multiscale phenomenon, we present a new double chain model, where we allow relaxation to occur by bending of the incommensurate coupled system of chains, as shown in Figure 1. We focus in particular on modeling and analyzing rigorously ripples occurring in ground state configurations, as well as their numerical simulation. Our model was  [19].
The model (2.6)-(2.8) developed in this work can be seen as a new application of the well-known Frenkel-Kontorova model for a one-dimensional atomic chain lying in an incommensurate periodic potential [1,2,6,4]. The role of the periodic potential is here assumed by the local interlayer contribution to the interaction energy, which is influenced by the local atom stacking configuration of the bilayer. Indeed, staggered atoms (maximizing the number of neighbors) minimize the energy compared to locally aligned configurations as seen on Figure 2 (see Section 1 for details). We will show how, by taking on some bending strain, the bilayer system is able to significantly reduce its energy by maximizing the area of staggered configurations.
Let us note that it is not possible for perfectly inextensible 2D sheets to be curved in more than one direction. As a result, a reasonable 2-dimensional model for the rippling of incommensurate layers should account for the slight extensibility of the individual layers. Thus, many simplifications used in this work, such as Eq. (2.3), cannot be used outside of our one-dimensional chain toy model. This paper is organized as follows. First, we motivate our approach by presenting a simple numerical atomistic simulation to illustrate the phenomenon and motivate our model (see also the related example in [19]). In Section 2, we introduce a mathematical model for the rippling mechanism, and we discuss the various assumptions that allow our explicit construction. In Section 3, we recall the classical analytical results on the Frenkel-Kontorova model due to Aubry and Le Daeron [2] that allow us to identify the ground state configurations. In Section 4, we show how the popular supercell method of approximating the incommensurate system by periodic approximations can be rigorously justified. Finally, we present in Section 5 numerical computations that illustrate and support our analytical results.
1. Motivating numerical example. We present first a simple numerical atomistic simulation to illustrate the phenomenon. Two coupled one-dimensional chains are fitted to represent a single zigzag row of dimers in graphene [19]. Within each chain, harmonic bonds model the interaction between nearest neighbors, and a harmonic angular spring is centered on each atom involving its two nearest neighbors. Between the two chains, atoms interact with a long-range Lennard-Jones potential 1 .Throughout  this section, we work in Lennard-Jones units: = 2.39 meV is the unit of energy and σ = 3.41Å is the unit of length for graphene-graphene interaction, see [9].
To simulate an incommensurate system, we create a periodic supercell by choosing the equilibrium length to be slightly different in each chain: we take N 1 = 232 and N 2 = 233 atoms respectively in each chain for a total length of 116 σ. Periodic boundary conditions are applied at the ends. Numerical parameters are fitted to match a zigzag row of dimers in graphene [12,19], leading to the parameterization given in Table 1. Note that the bond and angular spring constants k and k θ for each chain are related to the Young modulus Y and bending modulus D of graphene as where d is the width of the dimer row and l the equilibrium length of the bonds in the atomistic chain. In the initial configuration, each chain is in its isolated equilibrium state, and their relative distance is chosen as h = 1.063 σ. The total energy of the system is a function of all atomic positions R 1 i , i = 1, . . . , N 1 (bottom chain) and R 2 j , i = 1, . . . , N 2 (top chain) and also on the periodic cell length L, and can be divided into three components: however, it is also too smooth to correctly describe energy variations between different stacking configurations. Truly realistic computations should thus use a more accurate potential, such as the graphitic disregistry-dependent potential developed by Kolmogorov and Crespi [11]. For simplicity, and without loss of generality, we will use in this paper the Lennard-Jones potential, acknowledging that this actually leads to underestimating the formation of ripples.

3
The intralayer energies U 1 intra , U 2 intra take the form where it is understood that periodic boundary conditions are applied: R 1 0 , R 1 N1+1 , etc. refer to the appropriate atomic position in the adjacent periodic images. The interlayer energy U inter is where the inner sum runs over all the neighbors of R 1 i in the top layer that fall initially within the cutoff radius of 29 σ, whether in the simulation cell or the periodic images.
The numerical relaxation of atomic positions and total length of this coupled system leads to the rippled shape presented in Figure 1 and the partition of energy presented in Table 2. The total length of the system is reduced by 0.188%. A first observation from this computation is that the energy reduction achieved by bending appears to be about twice the potential energy in the angle springs in the relaxed configuration. Moreover, the two chains relax nearly inextensibly by producing a coherent ripple in the vertical direction by bending together. The inplane strain energy, modeled here by the potential energy in the bond springs, is an order of magnitude smaller than the other energies.
The driving mechanism for this relaxation process can be infered by measuring the disregistry of the atoms in layer 2 compared to the atoms in layer 1. This is a measure of the local atomic configuration, obtained by projecting the position of each atom in layer 2 on the curve describing layer 1. By measuring the distance of the projection to the nearest atom to the left on this curve, we can distinguish between favorable (staggered) configurations, ∆ = .25 σ, and unfavorable (aligned) configurations, ∆ = 0.
The disregistry obtained in the previous numerical simulation is plotted on Figure 3. Clearly, the relaxation tends to increase the number of atoms close to the most favorable configuration. On the other hand, the distance between the layers does not vary much.
Our conclusion is that this relaxation process can be correctly modeled by considering only the bending degrees of freedom and neglecting any in-plane strain in each layer. The driving mechanism behind the formation of coherent ripples is the competition between two main contributions: 1. the bending energy in the angular bonds due to chain curvature; 2. the interlayer potential energy variations between local configurations arising from the disregistry across the layers. 2. Derivation of the model. We will now develop a mathematical model for the determination of the ground state of a coupled system of parallel periodic atomic chains which are infinite and truly incommensurate: unlike the previous numerical example, there exists no periodic supercell for the coupled system. This is achieved by choosing an irrational number to be the ratio l 2 /l 1 of equilibrium length for the intralayer bonds.
2.1. Geometry. First, we describe the model geometry of the coupled system. Our main assumption is that the chains act inextensibly due to the strong in-plane modulus, as observed numerically in Section 1. Consecutive atoms in each chain will be separated by a fixed distance, respectively 1 for the chain C 1 and α for chain C α where α is an irrational real number. The chains are allowed to bend in the twodimensional plane, while they remain separated by a fixed distance h. Consequently, C 1 and C α form two infinite smooth parametric curves in the plane separated by a fixed distance, as seen on Figure 4(a).
To derive the equations of the model, we assume that the curvature of the system remains small and varies slowly. Consequently, we use the arc length along either curve to measure the distance between atoms in the respective chain. Furthermore, we choose to describe the whole system by following the curve C 1 . This leads us to introduce the following quantities.
Abscissa. Let s ∈ R denote the arc length parameter along the curve C 1 , and γ 1 (s) : R → R 2 be a natural parametrization of C 1 : γ 1 (s) = 1. Atoms on the chain Let the atoms of chain C α be positioned at points {R α j } j∈Z . Then, their projection on C 1 can be parameterized by a strictly increasing sequence of abscissas as seen on Figure 4(a), which we denote where n 1 (s) is the normal vector to the curve C 1 at abscissa s. Curvature. While the arc length between consecutive atoms R α j along C α is fixed as α, this is not the case for their projections on C 1 . As seen on Figure 4(b), the arc Curvature κ α j (a) Abscissa and curvature at projected positions. length between projections of consecutive atoms, i.e. s α j+1 − s α j for j ∈ Z, depends on the curvature of the curve C 1 , which we denote κ. Note that we choose the convention κ < 0 for the configurations presented in Figure 4. We obtain then, approximating locally the curvature by the constant value κ α j = κ(s α j , s α j+1 ) : Assuming that |hκ α j | 1, we will make the first-order approximation This leads to the formula for the curvature 2.2. Potential energy. To formulate our model, we consider only two contributions to the potential energy of the system of coupled chains: an interaction energy and a bending energy. For simplicity, the bending energy is modeled as a quadratic term in the curvature variable. It remains to construct a suitable approximation for the inter-chain potential energy as a function of the variables {s α j } j∈Z . Let V be an interatomic potential contributed by two atoms each in a different chain depending only on their relative distance, such as the Lennard-Jones potential. The total interaction energy is formally: In this expression, we have isolated the contribution of each atom in chain C α as an inter-chain potential energy created collectively by all atoms of chain C 1 . Due to the long-range nature of V , this inter-chain potential depends not only on the abscissa s α j of the projection of the atom on the curve, but also on the shape of the whole curve as we see on Figure 4(a). However, due to the fast decay of the inter-atomic potential, the main contribution to the inter-chain potential is local and can be deduced by approximating positions of the neighboring atoms belonging to the curved chain C 1 .
Let us construct such an approximation explicitely. For each j ∈ Z, a parabola is used as a locally second-order accurate approximation C 1 to the exact curve C 1 , in the neighborhood of the atom positioned at R α j . This parabola, parametrized as s → γ 1 (s), is uniquely determined (see Figure 5) by the distance h between the curves and the local curvature κ α j of curve C 1 , choosing R α j as origin of a local system of coordinates. The approximate atomic positions R 1 i = γ 1 (i) along this parabola are then deduced from the additional input of the projected abscissa s α j , using the arc length ∆ j = s α j (mod 1). This geometrical construction allows us to introduce now a disregistry-and curvature-dependent approximate inter-chain potential. For given values of s, κ ∈ R and h > 0, we obtain the sequence of positions { R 1 i } i∈Z as above, relative to the origin R of the coordinates. Then, we propose as an approximation to V 1→α (R) introduced in (2.4) the quantity: It is clear that the potential V per (s; κ, h) is a periodic function of the variable s with period 1. Note that in the following, we will drop the dependency of V per on the fixed variable h.
Remark 2.1. This construction of V per gives an explicit formula for the approximate potential. Hence, it can be used to precompute numerically and accurately an approximate inter-chain potential, for calculations aimed at finding ground state configurations.
The total potential energy of the full system is then modeled as the formal sum: where β is the angle spring constant modeling the resistance to bending of the coupled system of chains. Thanks to Eq. (2.3), this energy is formally where v(s α j , s α j+1 ), the local potential energy at site j, is given by 3. Ground state analysis. The model described by Eq. (2.6) or (2.7) is a generalized Frenkel-Kontorova model. There exists a vast literature on the Frenkel-Kontorova model, which we will not attempt to summarize here: see e.g. [4]. We only present the basic theory necessary for our analysis.
Following [2], we will make the following assumptions on the function v: 1. v(s, t) has a lower bound: there exists B such that (1, 1) and satisfies a symmetry condition:

v(s, t) is invariant by translation by
3. v(s, t) is a twice continuously differentiable function, and satisfies for all s, t ∈ R, where C is a finite positive constant. Note that for the function v defined by (2.8), conditions (3.1) and (3.2) are satisfied whenever V per is bounded from below. For the twist condition (3.3) to hold, stronger conditions are necessary, namely V per must be twice continuously differentiable with ∂ 2 v ∂s∂t bounded from below and β large enough.
Under these conditions, Aubry-Mather theory leads to the following description and results, for which details are to be found in the seminal paper [2]: Definition 1. A minimum energy configuration is a sequence {s j } such that any finite change δ j of a finite set of atoms necessarily increases the energy, that is: for any J < J, and any choice of δ j with δ j = 0 for j < J and j ≥ J.
The set of minimum energy configurations, which includes all possible boundary conditions, can be further divided between defect configurations and ground state configurations. To make this distinction, one derives the equilibrium equations (3.5) ∂v ∂s j (s j , s j+1 ) + ∂v ∂s j (s j−1 , s j ) = 0.
These equations can be interpreted as the motion equation of a dynamical system with discrete time j by introducing the conjugate variable for s j , Then (p j+1 , s j+1 ) are recursively defined from (p j , s j ) by the implicit system derived from (3.5): Note that under the assumptions above, ∂v ∂s (s, t) is a monotonic function of t for fixed s, and hence system (3.7) has a unique solution for a given (s j , p j ). One then introduces the bounded variable in the unit torus T, and then defines a non-linear operator T on the cylinder R × T onto the cylinder, By construction, T is an area-preserving one-to-one map. We then define the orbit of a point (p 0 , θ 0 ) as an infinite sequence (p j , θ j ) = T j (p 0 , θ 0 ) for −∞ < j < ∞.
Definition 2. An orbit (p j , θ j ) is recurrent if there exist an integer sequence j k → ∞ such that lim k→∞ (p j k , θ j k ) = (p 0 , θ 0 ) . Definition 3. 1. A defect configuration {s j } is a stationary configuration (satisfying (3.5)) such that the associated orbit { T j (p 0 , θ 0 )} is non-recurrent. 2. An elementary defect is a defect configuration which is also a minimum energy configuration.

A ground-state configuration is a minimum energy configuration which is not a defect configuration (the associated orbit is recurrent).
Using the previous definitions, the following two theorems, which were first stated and proved by Aubry and Le Daeron in [2], largely identify the minimum energy configurations and ground states of the Frenkel-Kontorova model: For any minimum energy configuration, there exists two real numbers l (the atomic mean distance) and ω (the phase), such that for any j, s j and jl + ω belong to the same interval [k j /2, (k j + 1)/2] where k j = Int(2s j ).
Note that the number l is also called the rotation number of the orbit corresponding to a minimum energy configuration {s j }. While it can in principle be arbitrarily chosen, in the case of the present model we are interested mainly in the choice l = α, which models the case of a bilayer system which is flat on average. Other choices could be interpreted as systems that have a global curvature.
1. The set of minimum energy configurations with atomic mean distance l is totally ordered, i.e., for two such configurations {s j } and {t j }, one of the assertions holds true for all j: s j < t j , s j = t j , or s j > t j . 2. The whole set of ground state configurations with irrational atomic mean distance l can be parameterized with one or two hull functions f which are strictly increasing: (a) When f is continuous, a unique function f parameterizes all ground states; (b) When f is discontinuous, two determinations f + and f − are necessary to parameterize all ground states. They are both strictly increasing and correspond to the right continuous or left continuous determinations of the same discontinuous function. (In this case, the set of discontinuity (d) Any ground-state configuration is determined by a phase ω and a determination of f , either f + or f − such that In particular, incommensurate ground-states are represented by trajectories which are rotating either on a smooth Kolmogorov-Arnold-Moser torus (when f is continuous) or on a discontinuous Cantor set or Cantorus (when f is discontinuous). A striking consequence is that ground-state trajectories always belong to the possibly zero-measure set of non-chaotic trajectories, in particular the ground-state ensemble is described by a Cantorus typically embedded in a chaotic region of the map. The transition between these two regimes has been called breaking of analyticity [2,1,6]

Approximation by periodic configurations: a partial justification.
Incommensurate ground-states are typically constructed as the limit of commensurate ground-states [1,2,6] with rotation number α n = p n /q n converging to the irrational α. However, in general only a subsequence of such configuration sequences converge.
In fact, we will see that when there exists a smooth invariant torus for T with rotation number α, then it is the limit (e.g. in the Hausdorff metric) of such periodic orbits; however, when the set of ground states with rotation number α forms a Cantor set, one may only prove that it is the limit of a particular sequence of commensurate ground states. For example, a given sequence may in principle converge to an incommensurate defect configuration, i.e., an orbit which is asymptotic to the Cantorus in both directions.
A validation of this approach, at least in the case of a smooth hull function, can be obtained by the use of perturbation theory. The starting point is given by the following results, see e.g. [7,22] for the proof: Proposition 6. Assume that T : T × R is a 2-dimensional symplectic map of class C r with r > 1. Suppose that T admits a C r (resp. analytical) invariant circle Γ, homotopic to T × {0}, on which the motion is C r -conjugate (resp. analytical) to a rigid rotation with rotation number α Diophantine of type (K, τ ), i.e. such that Then, given any nonnegative integer k < (r − 1)/τ , there exists a symplectic C r−1−kτ (resp. analytical) mapping Ψ defined in a neighborhood of Γ, mapping Γ to T × {0}, and having a C r−1−kτ (resp. analytical) inverse in a neighborhood of Γ, such that where ∆ : R → R and R : T × R → T × R are C r−1−kτ (resp. analytical) functions, ∆(0) = 0, and the remainder R satisfies in a neighborhood of Γ, Using the mapping ψ, we have thus a new set of coordinates (Φ,H) defined in a neighborhood of Γ, and such that Γ is the curveH = 0. In particular, the fact that T is approximately integrable in the vicinity of Γ leads to the following result: Lemma 7. If |α−p/q| is small enough, all periodic orbits of type p/q are contained in the strip where K depends only on the system and on the circle.
Proof. The idea of the proof, see [7] for details, is that since T is a small perturbation of an integrable mapping, there exists invariant K.A.M. tori whose rotation number α is such that |α − p/q| ≤ |α − p/q| and p/q belongs to the interval (α, α ). Then, because the map which to a rotation number associates the corresponding invariant circle is Lipschitz, the circle of rotation number α is contained in the strip Finally, as a consequence of the twist property, the periodic orbit of type p/q has to be contained between the circles of rotation number α and α.
We can now easily show that the computation of periodic ground states allows us to approximate the incommensurate ground states, i.e., the associated hull function.
To do so, let us consider a periodic minimum energy configuration, characterized equivalently by a finite sequence s 0 , . . . , s q−1 and s q = s 0 + p, or the initial point (p 0 , θ 0 ). Without loss of generality, we choose s 0 ∈ [0, 1). We define a function F q : [0, 1] → [s 0 , s 0 + 1] by associating: • to each point x j = {j p q } the value {s j − s 0 } + s 0 for j = 0, . . . , q − 1, where {·} denotes the fractional part, • to x q = 1 the value s 0 + 1, • finally, the periodic function x → F q (x)−x is written as the linear interpolant of the previously defined values: This function can then easily be extended to R by the rule F q (x + 1) = F q (x) + 1.
Proposition 8. Let the mapping T and a smooth invariant circle Γ be as in Proposition 6 with k ≥ 1. Then, for |α − p/q| small enough, there exists a constant C depending only on the system and on the circle, such that if f is a hull function for the set of ground states associated with Γ, there exists a phase ω ∈ R such that Remark 4.2. When α is of type (K, τ ) with τ large, there are values of q for which the dominant error term is the interpolation error, O 1/q 2 . The error bound can then be improved by using e.g. spline or trigonometric interpolation.
where the last estimate follows from the construction of F q and (4.4). Since f belongs to C r and by linearity of the interpolation, we obtain that the rate of convergence of the linear interpolant satisfies This result provides us with an a priori estimate for the convergence rate of the supercell approximations. It justifies in particular the use of the so-called artificial strain as a measure of the error introduced the system by stretching one or both of the chains to obtain a commensurate system.
Let us stress that the estimate is valid only in the case where the mapping T has a smooth invariant circle Γ with rotation number α, or equivalently, when the hull function is smooth (see Theorem 5). In general, when the bending modulus β decreases enough, the system will undergo a commensurate-incommensurate transition at a value β c which depends on α. The hull function then becomes discontinuous and the estimate breaks down.

Numerical examples.
We propose now to illustrate the previous analysis by some numerical examples, both in the smooth and the discontinuous case. For the purpose of computations in this section, we use the following set of parameters: where we recall that the distance between atoms is 1 for chain C 1 and α for chain C α , the distance between chains is h and the inter-atomic potential between atoms of different chains is chosen as the Lennard-Jones potential The ratio α is chosen here arbitrarily to satisfy the following conditions: first, the two chains have almost the same period, α ≈ 1, so we expect a large Moiré pattern and the formation of ripples, and second, α is of Diophantine type, see Proposition 6. Both conditions are satisfied here since the golden ratio (1 + √ 5)/2 is of Diophantine type ( √ 5, 2), see e.g. [10], and 13/8 is one of its rational approximant.

Realistic bending modulus: smooth ripples.
We present first the numerical results when choosing the realistic value β = 764 for the bending modulus β (see Section 1 and [12,19]). The conditions (3.1) and (3.3) are satisfied in this case. The numerical procedure is as follows. For any choice of rational approximant α ≈ p/q, we set up a periodic system with initial values s p,q j,0 = j p q , j = 1 . . . q, and we minimize numerically the energy per periodic cell given by (2.6). From the relaxed values obtained for {s p,q j }, we can then obtain approximations F q for the hull function by Formula (4.2). We can also reconstruct approximately the parametric curves modeling the chains using the curvature formula (2.3) and plot the ripple geometry.
We compute first a reference solution by using the periodic approximant p m = 2555, q m = 2566 for which the error (also called artificial strain) is approximately |α − p m /q m | ≤ 6 · 10 −8 .
Results are presented in Figures 6(a) and 6(b). Note that due to the scaling difference of the two axis, the oscillations of the chains are greatly exaggerated in Figure 6(a), giving the false impression of a greatly varying distance between the chains (see the orthonormal view on top). As expected, the relaxation process increases the number  of atoms in a favorable staggered configuration, s ≈ .5 (mod 1). We observe that both the shape of the ripples and the hull function show a smooth behavior. The overall shape of the twin curves seem almost perfectly periodic (even if the atoms are not distributed periodically on the curve) with the Moiré pattern period L M ≈ 106σ. This indicates that a purely continuous approach is likely to be succesfuly in capturing the overall behavior, since atomic details (here, the atomic disregistry) integrate out at the mesoscopic scale.
Next, we evaluate numerically the convergence behavior of the numerical approach. We consider all periodic approximants p/q to α with q ≤ 1000, where p is the closest integer to αq, and we compute for each the relaxed configuration {s p,q j } j=1...q . We evaluate then the error with respect to the reference configuration computed above, using the approximate error norm The result of this calculation is presented in Figure 6(c). The computed error falls rapidly on a single straight line as a function of decreasing α − p/q, showing good agreement with the estimate (4.3). This supports the intuitive idea that the artificial strain |α − p/q| controls the error, rather than the size q of the periodic supercell.
Remark 5.1. Note that the distance between layers h above is determined by minimizing the energy with respect to h only, e.g. for two flat layers, where the energy can be computed numerically.
Remark 5.2. It is also possible to deal with a varying distance h, if we make the assumption that this distance is locally determined as a function of the position of the atom S α j : indeed, this amounts to using a new potential given by the formula (see also (2.5) and Figure 5):

Weak bending modulus: nonsmooth behavior.
Finally, we present the results for a weakened system of chains with bending modulus β = 10. Note that the twist condition (3.3) is still valid in this case. The new results are presented in Figure 7. The relaxation pattern in this case is quite different from the previous case: the ripples present a sharp peak separating regions of almost constant curvature, following the Moiré pattern. This allows the two chains to form commensurate regions where the atoms are all in the favorable staggered position, separated by solitons or domain walls which are classical in the study of the Frenkel-Kontorova model [2,4].
The hull function is discontinuous in this case, with a large plateau at half-integer values corresponding to favorable staggered configurations and well-defined transition values forming smaller plateaus separated by forbidden gaps. This indicates a different mechanical behavior from the smooth case, as the sharp peaks are exponentially localized solitons, mechanically pinned by the Peierls-Nabarro potential [20,16,4].
Another difference is a small, but nonzero secondary oscillation of the system of ripples, which can be observed by comparing the vertical position of the maxima of the reconstructed chains on Figure 7(a). This indicates the onset of the discrete interaction between solitons and substrate (Peierls-Nabarro potential) as well as between solitons at the mesoscopic scale. This is in contrast to the smooth case where atomic details could be integrated out in a continuum description at the mesoscale.
Numerically, there is also a marked difference in the convergence behavior. We observe on Figure 7(c) that the computed errors do not fall onto a single line as in the smooth case. On average, we observe that the errors fall exponentially fast for |α − p/q| > 10 −5 , and then stagnate. Note that in this case, there is no smooth invariant circle with rotation number α for the map T (see section 4), so the estimate (4.3) does not hold. Hence, in principle some subsequences of the periodic approximations may converge to metastable defect configurations.  Conclusion. In summary, we have presented a new model for the spontaneous relaxation for free-standing systems of incommensurate van der Waals bilayers. This relaxation phenomenon is driven by the trade-off between bending energy and the weak interlayer interactions and the associated disregistry-dependent local configuration potential. Note that this is a different mechanism from the widely studied intrinsic monolayer rippling observed, e.g., in suspended graphene [14], which can be explained by phononic thermal fluctuations [8] and stress [3]. It differs also from ripples caused by interactions with a periodic substrate [18].
To study mathematically the relaxation of such systems, we have introduced a new one-dimensional double chain model in which consecutive atoms in each chain are constrained to remain at a fixed distance. Relaxation occurs through spontaneous bending of the chains to minimize the interlayer interaction energy. This model can be seen as a new application of the generalized Frenkel-Kontorova model [4].
As a result, we have identified the ground state configurations of the system, and we have given rigorous error estimates for the application to our model of the popular supercell approach used to approximate incommensurate systems. We have also presented numerical results supporting the analysis, both in the case of a realistic bending modulus where smooth quasi-periodic ripples can be numerically observed, and for a weaker bending modulus where commensurated domains separated by sharply peaked solitons appear. The transition between the two regimes is the famous commensurateincommensurate transition [2].
The next step will be to extend this model to the two-dimensional case. A particular difficulty is that perfectly inextensible two-dimensional sheets cannot be curved in more than one direction, and thus the model should be extended to account for small, but nonzero, in-plain strains.
6. Acknowledgments. This work was supported in part by ARO MURI Award W911NF-14-1-0247. Mitchell Luskin was also supported in part by the Radcliffe Institute for Advanced Study at Harvard University.