A Discrete-to-Continuum Model of Weakly Interacting Incommensurate Chains

TL;DR

This paper develops a continuum variational model from a discrete atomic chain model to describe incommensurate layered materials like graphene on a substrate, capturing key behaviors observed in the discrete system.

Contribution

It introduces a formal discrete-to-continuum derivation for incommensurate atomic chains, bridging microscopic and macroscopic descriptions.

Findings

Continuum model accurately reproduces discrete chain behavior.

Large commensurate regions separated by incommensurate domain walls.

Numerical solutions confirm qualitative and quantitative agreement.

Abstract

In this paper we use a formal discrete-to-continuum procedure to derive a continuum variational model for two chains of atoms with slightly incommensurate lattices. The chains represent a cross-section of a three-dimensional system consisting of a graphene sheet suspended over a substrate. The continuum model recovers both qualitatively and quantitatively the behavior observed in the corresponding discrete model. The numerical solutions for both models demonstrate the presence of large commensurate regions separated by localized incommensurate domain walls.

Figures

Equations

$$E({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2})=E_s({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2})+E_b({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2})+E_w({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2}),$$

$$E_s({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2})=\sum _{i=1}^{N_2-1}\frac{k_s}{2}{\left(\frac{\Vert {\mathbf{b}}_i\Vert -h_2}{h_2}\right)}^2,$$

$$E_b({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2})=\sum _{i=1}^{N_2-2}\frac{k_b}{2}{\left({\theta }_{i+1}-{\theta }_i\right)}^2,$$

$${\mathbf{q}}_{i+1}-{\mathbf{q}}_i=\Vert {\mathbf{b}}_i\Vert (\cos ({\theta }_i),\sin ({\theta }_i)),$$

$$E_b({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2})=k_b\sum _{i=2}^{N_2-1}\frac{\Vert {\mathbf{b}}_i\Vert \Vert {\mathbf{b}}_{i-1}\Vert -{\mathbf{b}}_i\cdot {\mathbf{b}}_{i-1}}{\Vert {\mathbf{b}}_i\Vert \Vert {\mathbf{b}}_{i-1}\Vert },$$

$$E_w({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N_2})=\omega \sum _{i=1}^{N_2}\sum _{j=1}^{N_1}g\left(\frac{\Vert {\mathbf{p}}_j-{\mathbf{q}}_i\Vert }{\sigma }\right),$$

$$g(r)=r^{-12}-2r^{-6}.$$

$${\mathcal{C}}_1^0=\{(s,0):s\in [0,L]\}\text{mbox}and{\mathcal{C}}_2^0=\left\{(s,\sigma ):s\in [0,L]\right\}.$$

$$\left\{(s+u(s),\sigma +v(s)):s\in [0,L]\right\}.$$

$$t=\frac{s}{L},\bar{\xi }=\frac{u}{L},\bar{\eta }=\frac{v}{L},\bar{E}=\frac{E}{\omega }$$

$$\epsilon =\frac{\sigma }{L},{\delta }_1=\frac{h_1}{\sigma },{\delta }_2=\frac{h_2}{\sigma },{\gamma }_e=\frac{k_s}{\omega {\delta }_2},{\gamma }_b=\frac{k_b{\delta }_2}{\omega }.$$

$${\mathcal{C}}_1=\left\{(t,0):t\in [0,1]\right\}\text{mbox}and{\mathcal{C}}_2=\left\{(t+\bar{\xi }(t),\epsilon +\bar{\eta }(t)):t\in [0,1]\right\}.$$

$${\delta }_1-{\delta }_2=\mathcal{O}(\epsilon ),$$

$$\Gamma (t)=(t+\epsilon \xi (t),\epsilon +\epsilon \eta (t))\text{mbox}witht\in [0,1].$$

$$\begin{gathered} {\mathbf{b}}_i=\epsilon ({\delta }_2+\xi (t_{i+1})-\xi (t_i),\eta (t_{i+1})-\eta (t_i)) \\ =\epsilon {\delta }_2\left(1+\epsilon {\xi }^{\prime }(t_i)+\frac{{\delta }_2{\epsilon }^2}{2}{\xi }^{\prime \prime }(t_i)+\frac{{\delta }_2^2{\epsilon }^3}{6}{\xi }^{\prime \prime \prime }(t_i),\epsilon {\eta }^{\prime }(t_i)+\frac{{\delta }_2{\epsilon }^2}{2}{\eta }^{\prime \prime }(t_i)+\frac{{\delta }_2^2{\epsilon }^3}{6}{\eta }^{\prime \prime \prime }(t_i)\right)+o\left({\epsilon }^3\right). \end{gathered}$$

$$\begin{gathered} {\bar{E}}_s[\xi ,\eta ]=\frac{{\gamma }_e{\delta }_2{\epsilon }^2}{2}\sum _{i=1}^{N_2-1}\{{\left[{\xi }^{\prime }(t_i)+\frac{\epsilon }{2}\left({\left({\eta }^{\prime }(t_i)\right)}^2+{\delta }_2{\xi }^{\prime \prime }(t_i)\right)\right]}^2 \\ -\frac{{\epsilon }^2}{3}\left[3{\left({\xi }^{\prime }(t_i){\eta }^{\prime }(t_i)\right)}^2-3{\delta }_2{\xi }^{\prime }(t_i){\eta }^{\prime }(t_i){\eta }^{\prime \prime }(t_i)-{\delta }_2^2{\xi }^{\prime }(t_i){\xi }^{\prime \prime \prime }(t_i)\right]\}+o\left(N_2{\epsilon }^4\right), \end{gathered}$$

$${\bar{E}}_b[\xi ,\eta ]=\frac{{\gamma }_b{\delta }_2{\epsilon }^4}{2}\sum _{i=2}^{N_2-1}{\left[{\eta }^{\prime \prime }(t_i)\right]}^2+o\left(N_2{\epsilon }^4\right).$$

$$\begin{gathered} {\bar{E}}_s[\xi ,\eta ]\sim \frac{{\gamma }_e\epsilon }{2}{\int }_0^1\{{\left[{\xi }^{\prime }(t)+\frac{\epsilon }{2}\left({\left({\eta }^{\prime }(t)\right)}^2+{\delta }_2{\xi }^{\prime \prime }(t)\right)\right]}^2 \\ -\frac{{\epsilon }^2}{3}\left[3{\left({\xi }^{\prime }(t){\eta }^{\prime }(t)\right)}^2-3{\delta }_2{\xi }^{\prime }(t){\eta }^{\prime }(t){\eta }^{\prime \prime }(t)-{\delta }_2^2{\xi }^{\prime }(t){\xi }^{\prime \prime \prime }(t)\right]\}dt=:{\mathcal{E}}_s^{\epsilon }[\xi ,\eta ] \end{gathered}$$

$${\bar{E}}_b[\xi ,\eta ]\sim \frac{{\gamma }_b{\epsilon }^3}{2}{\int }_0^1{\left[{\eta }^{\prime \prime }(t)\right]}^{2}dt=:{\mathcal{E}}_b^{\epsilon }[\xi ,\eta ].$$

$${\mathcal{E}}_s^{\epsilon }[\xi ,\eta ]=\frac{{\gamma }_e\epsilon }{2}{\int }_0^1\left\{{\left({\xi }^{\prime }(t)+\frac{\epsilon }{2}{\left({\eta }^{\prime }(t)\right)}^2\right)}^2-{\epsilon }^2{\left({\xi }^{\prime }(t){\eta }^{\prime }(t)\right)}^2-\frac{{\epsilon }^2{\delta }_2^2}{12}{\left({\xi }^{\prime \prime }(t)\right)}^2\right\}dt.$$

$${\mathcal{E}}_w^{\epsilon }[\xi ,\eta ]=\frac{1}{\epsilon }{\int }_0^{1}G\left(t,\xi ,\eta \right)dt.$$

$$\sum _{j=-\infty }^{\infty }g(d_{ij}/\epsilon ),$$

$$d_{ij}=\sqrt{{\left(j\epsilon {\delta }_1+k_i+\epsilon \xi \left(\epsilon {\delta }_{2}i\right)\right)}^2+{\left(\epsilon +\epsilon \eta \left(\epsilon {\delta }_{2}i\right)\right)}^2}.$$

$$k_i=\text{mod}(i\epsilon {\delta }_2,\epsilon {\delta }_1)=i\epsilon ({\delta }_2-{\delta }_1)-i_1\epsilon {\delta }_1=i\epsilon {\delta }_2({\delta }_1-{\delta }_2)/{\delta }_2-i_1\epsilon {\delta }_1,$$

$$k(t):=\epsilon \alpha t,$$

$$G(t,\xi (t),\eta (t))$$

$$=\sum _{j=-\infty }^{\infty }g\left(\frac{\sqrt{((j-i_1)\epsilon {\delta }_1+k(t)+\epsilon \xi (t){)}^2+(\epsilon +\epsilon \eta (t){)}^2}}{\epsilon }\right)$$

$$=\sum _{j=-\infty }^{\infty }g\left(\sqrt{(j{\delta }_1+\alpha t+\xi (t){)}^2+(1+\eta (t){)}^2}\right),$$

$$\begin{gathered} {\mathcal{E}}^{\epsilon }[\xi ,\eta ]:=\frac{{\gamma }_s\epsilon }{2}{\int }_0^1\left\{{\left({\xi }^{\prime }(t)+\frac{\epsilon }{2}{\left({\eta }^{\prime }(t)\right)}^2\right)}^2-{\epsilon }^2{\left({\xi }^{\prime }(t){\eta }^{\prime }(t)\right)}^2-\frac{{\epsilon }^2{\delta }_2^2}{12}{\left({\xi }^{\prime \prime }(t)\right)}^2\right\}dt \\ +\frac{{\gamma }_b{\epsilon }^3}{2}{\int }_0^1{\left({\eta }^{\prime \prime }\right)}^{2}dt+\frac{1}{\epsilon }{\int }_0^{1}G(t,\xi ,\eta )dt. \end{gathered}$$

$${\mathcal{E}}^{\epsilon }[\xi$$

$$-{\gamma }_s\epsilon ({\xi }^{\prime \prime }+\epsilon {\eta }^{\prime }{\eta }^{\prime \prime })+{\epsilon }^{-1}{G}_{\xi }(t,\xi ,\eta )$$

Full text

1 Introduction

A graphene sheet is a single-atom thick sheet of carbon atoms arranged in a hexagonal lattice. Graphene has exceptional physical properties and yields insights into the fundamental physics of two-dimensional materials. These facts have motivated an extensive effort to model and simulate graphene and other related carbon nanostructures.

The study in this paper is motivated by discrete-to-continuum modeling that describes the deformation of a graphene sheet suspended over a substrate. For a suspended graphene sheet, the positions of the atoms on the sheet are determined by the strong, bonded interactions between nearest neighbors on the sheet and by the weak, non-bonded interactions with nearby atoms on the substrate. A mismatch between the geometries of the hexagonal lattice of graphene and the substrate lattice can induce strain in the sheet. This strain can be relaxed by both in-plane and out-of-plane, atomic-scale displacements of the atoms on the sheet.

Insight into the response of graphene to this lattice mismatch can be obtained within the framework of the classical Frenkel-Kontorova theory [1]. For slightly mismatched one-dimensional lattices, this theory predicts relatively large commensurate regions separated by localized incommensurate regions. Hence we expect for suspended graphene that the local adjustments of atoms on the sheet may create large domains where the two lattices are commensurate and the interlayer energy is minimized. At the same time, there should be localized incommensurate regions, or domain walls, where strain may be relaxed by out-of-plane displacement.

Commensurate regions separated by domain walls have been observed in simulations of relaxed moiré patterns in graphene sheets [2, 3]. When two lattices with different lattice geometries or the same geometry but different orientations are stacked, a larger periodic pattern, called a moiré pattern, emerges. These moiré patterns are a strictly visual effect. However, if the atoms on one or both of the lattices are then relaxed to accommodate the mismatch between the lattices, additional patterns can occur. In [2], the authors study these relaxed moiré patterns by simulating interacting, identical graphene lattices where one lattice is slightly rotated with respect to the other. In [3], the authors report on similar simulations for a slightly rotated graphene lattice interacting with a hexagonal boron nitride substrate, which also has the structure of a hexagonal lattice with a slightly larger lattice constant than that of graphene. In both papers, simulations predict a two-dimensional pattern of wrinkles intersecting at small, misfit regions, called hot spots, exhibiting large out-of-plane displacements. The wrinkles separate large, flat domains of commensurate regions.

For the model we develop in this paper, the starting point is a pair of parallel curves that ‘carry’ atoms. The lower curve models a rigid substrate, and its atoms are fixed at a prescribed spacing. The upper curve models a graphene sheet, and the atoms on the upper curve can displace. The equilibrium spacing between the atoms on the upper curve is also prescribed. The two lattices are mismatched if the interatomic spacing between the atoms on the lower curve does not equal the equilibrium spacing between atoms on the upper curve.

We start with a discrete description of the stretching and bending energies of the atoms on the upper curve and the weak interaction energy between the atoms on the upper and lower curves, respectively. From this we attain a continuum description when a small parameter, representing the ratio of the typical spacing between the curves to the length of the curves, is small but not equal to zero. The minimizers of this continuum energy represent equilibrium configurations of the upper curve.

The continuum energy that we obtain has a Ginzburg-Landau-type structure with the elastic contribution that corresponds to the classical Föppl-von Kármán energy [4]. A number of recent studies have considered related problems of wrinkling of thin elastic sheets that are, e.g., bonded to a compliant substrate with a large compressive misfit [5, 6, 7], pulled down by the force of gravity [8], or are floating on a fluid [9, 10, 11].

The novel aspect of the continuum energy derived in the present work is the potential describing the weak interactions between the curves. Although a continuum description, the weak energy retains information about the local mismatch between the original discrete lattices. The presence of lattice mismatch constitutes the principal difference between the adhesion forces that are associated with the weak interactions between the lattices and the forces that are usually considered in the standard problem for a thin film bonded to an adhesive substrate. The wrinkling patterns in a typical adhesion problem develop primarily due to the misfit strains. On the other hand, the wrinkles observed in incommensurate lattices appear not only because of the mismatch in strains, but also in order to accommodate the difference between the number of atoms on the two chains while maximizing the regions of registry. For example, it does not make a difference whether the deformable chain has one more or one fewer atom than the rigid chain per period—in both cases, a single wrinkle will appear wherever there is an extra atom or a ”vacancy” on the deformable chain. The out-of-plane deflection at this position occurs because the equilibrium distances between the chains are different if the chains are in registry or out of registry.

From our total continuum energy we derive the Euler-Lagrange equations, which are then solved numerically. We present some basic comparisons between discrete simulations and the predictions of our continuum model. For curves corresponding to slightly mismatched chains, our model predicts large commensurate regions separated by domain walls, or wrinkles, formed by localized out-of-plane ridges. The spacing of these domain walls can be determined by the need to accommodate a certain number of extra atoms on the upper curve. Qualitatively, our solutions exhibit a pattern of wrinkles similar to the predictions of the atomistic simulations in [2].

The discrete-to-continuum modeling in this paper is analogous to the approach taken in [12], in which the authors derive a continuum theory of multi-walled carbon nanotubes by upscaling an atomistic model. The atomistic model includes a bending energy related to strong covalent bonds between atoms in the same wall and an interaction energy between the atoms in adjacent walls. Numerical solutions using the continuum model that results from upscaling show that, for sufficiently large radii, the cross-section of a double-walled nanotube polygonizes. In this setting, the straight sections of the polygonal cross-section are commensurate regions. The corners of the polygon are the domain walls, and are analogous to the localized ridges predicted by the model developed in this paper.

A discrete model, similar to the one considered in this paper, was recently used in [13] to demonstrate numerically that the spontaneous atomic-scale relaxation of free-standing systems of incommensurate van der Waals bilayers leads to a simultaneous long-range rippling of the bilayer system.

This paper is organized as follows. In Section 2, we formulate a discrete energy of the system of a graphene sheet over a substrate. In Section 3, we derive a continuum energy that keeps track of the mismatch of the spacing between the atoms on each curve. The next section includes numerical results that compare the atomistic model with the continuum model. Furthermore, in this section we show how the different parameters give rise to different material deformations.

2 Atomistic Model

For simplicity, here we model the formation of isolated wrinkles in a graphene layer supported by a substrate within the framework of a one-dimensional model. Note that a similar description can be readily developed for graphene bilayers (cf. [12]).

Suppose that we have a discrete system that consists of two chains of atoms $\hat{\mathcal{C}}_{1}$ and $\hat{\mathcal{C}}_{2}$ that are $L$ units long. The atoms on the bottom chain $\hat{\mathcal{C}}_{1}$ are spaced $h_{1}$ apart and cannot move. This chain describes a rigid substrate. The atoms on the top chain $\hat{\mathcal{C}}_{2}$ can move. This chain represents a deformable graphene layer that is nearly inextensible and has a finite resistance to bending. Each atom on the top chain interacts with its two nearest neighbors via a strong bond potential, represented here by a stiff linear spring such that the equilibrium spacing between the atoms on $\hat{\mathcal{C}}_{2}$ is $h_{2}$. The resistance to bending is modeled by torsional springs between adjacent bonds. Further, all atoms on the first chain are assumed to interact with all atoms on the second chain via interatomic van der Waals potential. In what follows we will refer to $\hat{\mathcal{C}}_{1}$ as the rigid chain and to $\hat{\mathcal{C}}_{2}$ as the deformable chain.

We assume that, in the reference configuration (Figure 1(a)), the chains are parallel and separated by a distance $\sigma$. Here $\sigma$ is equal to the equilibrium distance between two atoms interacting via the van der Waals forces. Note that this reference configuration is not in equilibrium. If the ratio between $\sigma$ and the equilibrium bond length $h_{1}$ is large enough, from the point of view of an atom on $\hat{\mathcal{C}}_{2}$, its van der Waals interaction with all atoms on $\hat{\mathcal{C}}_{1}$ can be represented by an interaction with the curve representing $\hat{\mathcal{C}}_{1}$ with a uniform atomic density [14]. In equilibrium, the curves $\hat{\mathcal{C}}_{1}$ and $\hat{\mathcal{C}}_{2}$ are then given by the two straight parallel lines. The distance between these lines should be slightly smaller than $\sigma$ in order to accommodate attractive forces from more distant atoms. This approximation, however, ignores possible registry effects that are significant in determining the shape of the deformable chain $\hat{\mathcal{C}}_{2}$. In fact, the only situation in which the two straight parallel chains would correspond to an equilibrium configuration is when $h_{1}=h_{2}$. Indeed, in this case, all atoms on $\hat{\mathcal{C}}_{2}$ would occupy the positions above the midpoints between the atoms on $\hat{\mathcal{C}}_{1}$ and the system would be in global registry.

Here we are concerned with the situation when $h_{1}\neq h_{2}$, but $|h_{1}-h_{2}|/h_{1}\ll 1$. Under these assumptions, global registry cannot be achieved in an undeformed configuration because the distance between midpoints of neighboring intervals in $\hat{\mathcal{C}}_{1}$ is not equal to $h_{2}$. It follows that in order to achieve equilibrium, the deformable chain would have to adjust by some combination of bending and stretching.

Let the current positions of the $N_{2}$ atoms on the chain be given by the vectors $\left\{{\bf q}_{1},\dots,{\bf q}_{N_{2}}\right\}\subset\mathbb{R}^{2}$. For every $i=1,\ldots,N_{2}-1$, we represent the bond between the atom $i$ and the atom $i+1$ by the vector ${\bf b}_{i}={\bf q}_{i+1}-{\bf q}_{i}$. Then, the total energy of the system is given by

[TABLE]

where the stretching energy is defined by a harmonic potential

[TABLE]

with $k_{s}$ being the spring constant. The bending between the adjacent links of the chain is penalized by introducing torsional springs connecting these links, and therefore, the associated bending energy is given by

[TABLE]

where $k_{b}$ is the torsional constant and $\theta_{i}$ is the angle between the $i$-th link and the $x$-axis defined by

[TABLE]

for every $i=1,\ldots,N_{2}-1$. Assuming that $\left|\theta_{i+1}-\theta_{i}\right|\ll 1$ for all $i=1,\ldots,N_{2}-1$, in the sequel we will consider the expression

[TABLE]

for the bending energy that is equivalent to (2) to leading order.

The energy of the weak van der Waals interaction between $\hat{\mathcal{C}}_{1}$ and $\hat{\mathcal{C}}_{2}$ is defined by

[TABLE]

where $g$ is a given weak pairwise interaction potential and ${\bf p}_{j}=(p_{j},0)\in\mathbb{R}^{2},\ j=1,\ldots,N_{1}$ are the positions of the atoms on the rigid chain. The parameters $\sigma$ and $\omega$ define the equilibrium interatomic distance and the strength of the potential energy (4), respectively. In what follows, we assume that $g$ is the classical Lennard-Jones 12-6 potential given by

[TABLE]

3 Continuum Model

As a first step in formally deriving a continuum model, we assume that the chains of atoms are embedded in two sufficiently smooth curves $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. See Figure 1(a). Hence, the lower curve $\mathcal{C}_{1}$ is straight and rigid, while the upper curve $\mathcal{C}_{2}$ can deform. We denote these two curves in the reference configuration by $\mathcal{C}^{0}_{1}$ and $\mathcal{C}^{0}_{2}$. We assume

[TABLE]

Letting $h_{i},\ i=1,2,$ denote the equilibrium spacing for the atoms on $\mathcal{C}^{0}_{i}$, we assume that the atoms on $\mathcal{C}^{0}_{i}$ are a distance $h_{i}$ apart in the reference configuration. Note, that we do not assume that the system of two chains is stress-free in the reference configuration.

We let $(u(s),v(s))$ be the displacement of the point $(s,\sigma)$ on $\mathcal{C}^{0}_{2}$. Hence the deformed curve $\mathcal{C}_{2}$ is given by

[TABLE]

In particular, an atom at the point $(s_{i},\sigma)$ on $\mathcal{C}^{0}_{2}$ is displaced to the point $(s_{i}+u(s_{i}),\sigma+v(s_{i}))$. See Figure 1(b).

We assume $\sigma<<L$, i.e., that the spacing between the chains is much less than the length of the chains. To exploit this, we introduce the rescalings

[TABLE]

and the nondimensional parameters

[TABLE]

We obtain with a slight abuse of notation that

[TABLE]

We assume that $\delta_{i}=\mathcal{O}(1),\ i=1,2$, i.e., the distance between the atoms within each chain is comparable to the distance between the chains (and hence both are much smaller than the length of the chains). Furthermore, in order to observe the registry effects on a macroscale, we assume that

[TABLE]

so that the mismatch between the equilibrium spacings of the two chains is small.

We assume $\bar{\xi}=\varepsilon\xi$, $\bar{\eta}=\varepsilon\eta$, so that $\mathcal{C}_{2}$ is parameterized by

[TABLE]

See Figures 2(a) and 2(b). These scalings for the displacements are appropriate for small deformations considered here and eventually will lead to expressions for the strains similar to those for Föppl-von Kármán theory.

In the rescaled coordinates, the atoms in $\mathcal{C}_{2}^{0}$ are located at the points ${\bf q}_{i}^{0}=(t_{i},\varepsilon)$, were $t_{i}=\varepsilon\delta_{2}i$ for $i=1,\ldots,N_{2}$. The $i$-th atom is then displaced to the point ${\bf q}_{i}=(t_{i}+\varepsilon\xi(t_{i}),\varepsilon+\varepsilon\eta(t_{i}))$ for every $i=1,\ldots,N_{2}$.

3.1 Elastic Energy Contribution

We now take advantage of the fact that $\varepsilon\ll 1$. By using the Taylor expansions of $\xi(t_{i})$ and $\eta(t_{i})$ in $\varepsilon$, the bond ${\bf b}_{i}$ between the atom $i$ and the atom $i+1$ can be expressed in the form

[TABLE]

Substituting these expansions into the expressions (1) and (3) for the stretching and bending energies, respectively, we can redefine both energies in terms of values of $\xi$ and $\eta$ at $t_{i}$, where $i=1,\ldots,N_{2}-1$. Preserving the notation for the energies, we have that

[TABLE]

and

[TABLE]

Here the number of terms in the expansion of $\bar{E}_{s}$ is chosen so as to match the power of $\varepsilon$ in the lowest order term in the expansion of $\bar{E}_{b}$.

We now recall that $\mathcal{C}_{2}^{0}$ has a unit length in nondimensional coordinates and that the spacing between the atoms is equal to $\varepsilon\delta_{2}\ll 1$. Hence, the number of atoms in $\mathcal{C}_{2}^{0}$ is approximately $\frac{1}{\varepsilon}$ and therefore

[TABLE]

and

[TABLE]

From now on, we assume that the admissible functions $\xi$ and $\eta$ satisfy periodic boundary conditions, i.e., $\xi(0)=\xi(1)$, $\eta(0)=\eta(1)$, $\xi^{\prime}(0)=\xi^{\prime}(1)$, $\eta^{\prime}(0)=\eta^{\prime}(1)$, etc. With this assumption, after integrating by parts, the stretching energy in (15) can be written as

[TABLE]

Here the term ${\left(\xi^{\prime}(t)+\frac{\varepsilon}{2}{\left(\eta^{\prime}(t)\right)}^{2}\right)}^{2}$ is a one-dimensional version of the Föppl-von Kármán energy describing large deflections of thin flat plates. Note that the remaining higher-order nonlinear elastic terms coupling the vertical and horizontal components of deformation are nonpositive and thus problematic from the point of view of establishing existence of a solution to the corresponding variational problem. In Section 3.3 we argue that neglecting these terms should not affect the global behavior of a minimizer of the full continuum energy.

3.2 Van der Waals Energy Contribution

We now discuss the contribution to the energy from the van der Waals interactions, that is, the continuum version of (4), which has the form of

[TABLE]

The novelty of our model is in defining a function $G$ that gives a continuum description of the mismatch of the spacing between the atoms on each chain. We shall first focus in the inner sum of (4) and try to estimate the interaction of a given atom on the deformable chain with all the atoms on the rigid chain. We accomplish this by tracking the offset between the two lattices embedded in $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ as a function of $t$.

Our starting point is to pick an atom $i$ on $\mathcal{C}_{2}$. A discrete description of the total interaction energy between this atom and the atoms on $\mathcal{C}_{1}$ is given by

[TABLE]

where $d_{ij}$ is the distance between the fixed atom $i$ on $\mathcal{C}_{2}$ and the atom $j$ on $\mathcal{C}_{1}$. Here we replace the finite rigid chain with a larger rigid chain that contains infinitely many atoms. Due to the rapid decay of the potential function $g$ as it argument tends to $\infty$, this should have no effect on the total interaction energy unless the atom $i$ lies very close to an endpoint of the deformable chain. Without loss of generality, we will assume that, in the reference configuration, the leftmost atom on the deformable chain $\mathcal{C}_{2}^{0}$ lies directly above an atom on $\mathcal{C}_{1}$.

To write down an expression for $d_{ij}$, we let $k_{i}$ denote the offset between the atomic lattices of the two chains in the reference configuration as measured at the atom $i$ on $\mathcal{C}_{2}^{0}$. Then $k_{i}$ is the distance between the projection of the atom $i$ onto $\mathcal{C}_{1}$ and the first atom on $\mathcal{C}_{1}$ to the left of this projection. We will label this atom on the rigid chain as the atom with the index $j=0$. See Figure 3. It is now clear that

[TABLE]

Furthermore, since the atom at the left end of $\mathcal{C}_{2}^{0}$ is directly above the [math]-th atom on $\mathcal{C}_{1}$, if the horizontal component of the position of the atom $i$ on $\mathcal{C}_{2}$ is $t_{i}=\varepsilon\delta_{2}i$, then

[TABLE]

for some $i_{1}\in\mathbb{N}$. See Figure 4. Hence for any $t$ we define

[TABLE]

where $\alpha:=(\delta_{2}-\delta_{1})/\varepsilon\delta_{2}=\mathcal{O}(1)$. And therefore, our function $G$ that gives a continuum description of the mismatch of the spacing between the atoms on each chain is defined by

[TABLE]

where we translated the index $j$ by $i_{1}$ on the last step.

3.3 Continuum Energy

Putting together all the contributions, one choice for the continuum energy of the system is

[TABLE]

Note that the variational problem associated with the energy (24) is not well-posed since the negative sign in front of ${\left(\xi^{\prime\prime}\right)}^{2}$ leads to $\min{\mathcal{E}^{\varepsilon}}=-\infty$. This issue is due to our choice to terminate the expansion (12) of the bond vector at $O\left(\varepsilon^{4}\right)$; once the additional terms are included, the coefficient in front of the term containing the third derivative of $\xi$ squared will be positive. The downside of incorporating higher derivatives, however, is that the model quickly becomes more complicated.

From the Ginzburg-Landau structure of the energy in (24), if we expect the minimizers of (24) to develop wrinkles of characteristic width $\varepsilon$, then we should expect all derivatives of the minimizers to appear as powers of $\varepsilon^{-1}$ inside the wrinkled regions. Accordingly, all elastic terms in (24) will then contribute roughly the same amount to the overall energy. In fact, a quick glance at (12), indicates that all terms in that expansion can also have the same magnitude inside the wrinkles, possibly invalidating the asymptotic procedure that led to (24).

On the other hand, in the regions between the wrinkles, all minimizers should have bounded derivatives and the terms with higher powers of $\varepsilon$ in (24) should simply provide small corrections to the lower order contributions. The situation is not unlike that arising in continuum modeling of crystalline solids, where the structural defects—such as dislocations—can only be described in terms of their influence on the global strain field, without properly resolving the cores of the defects. In order to fully resolve these cores, an alternative approach would be to use a quasicontinuum method that combines continuum description in the bulk with the discrete description in a vicinity of each defect [15].

With these observations in mind, we will adopt the following strategy to formulate a purely continuum theory. Since the negative terms in (24) should only contribute inside the wrinkles and since the asymptotic procedure that led to (24) is likely not to be valid inside the wrinkled regions, we neglect these terms from now on. We conjecture that the overall effect of this change on the structure of minimizers will be restricted to a perturbation in the cross-section of the wrinkles. Numerical simulations in the next section appear to support this statement and, in fact, show that the predictions of the proposed continuum model are close to the results of discrete simulations, even inside the wrinkles. We set

[TABLE]

The energy functional in (25) can be viewed as a generalization of similar energies for one-dimensional Frenkel-Kontorova chains with the elastic contribution similar to that in Föppl-von Kármán theory.

Note that both the discrete and continuum nondimensional models contain the small parameter $\varepsilon$ and, in particular, the continuum model cannot be thought of as a limit of the discrete model as $\varepsilon\to 0$. Instead, we conjecture that both models converge to the same asymptotic limit as $\varepsilon\to 0$ in the appropriate sense. The limit has to be understood within the framework of $\Gamma$-convergence [16] so that both energies are $\Gamma$-equivalent [17]. Hence the number of terms retained in the expansion of the discrete problem in order to obtain the continuum problem should be sufficient to reproduce the behavior of the discrete system for a small $\varepsilon$. The proof of $\Gamma$-equivalence is a subject of a future work.

The Euler-Lagrange equations for (25) are

[TABLE]

Here (26) corresponds to the horizontal force balance and the function $G_{\xi}$ describes how the tensile/compressive force in the chain changes with axial position. Similarly, (27) is the vertical force balance, including the vertical force $G_{\eta}$ that arises from coupling between the chains mismatch and the out-of-plane displacement of the deformable chain. Note that the vertical deformation in the von Kármán contribution to (25) is multiplied by the small parameter $\varepsilon$, in contrast with a typical form of this expression.

The natural boundary conditions are as follows

[TABLE]

These boundary conditions have natural interpretations in elasticity theory: (28) corresponds to the assumption of no applied compression, while (29) and (30) indicate that there is no applied moment or shear force. The boundary conditions and the force balances (26)-(27) can be modified in the standard way to include applied loads and body forces.

4 Numerical Results

In this section, we study the predictions of our continuum model by numerically solving the 2-point boundary-value problem defined by the Euler-Lagrange equations (26) and (27) with the boundary conditions (28)–(30).

We make several basic comparisons between the predictions of our continuum model and those of the discrete model. The discrete simulations are conducted by using dissipation-dominated (gradient-flow) dynamics based on the discrete energy described in Section 2. Recall that this energy has 3 terms that correspond to the stretching (1) and bending (2) energies of the deformable chain and the van der Waals interaction between the deformable and the rigid chains (4). Assuming gradient-flow dynamics gives a system of ordinary differential equations. These are solved numerically until the solution equilibrates, which yields solutions that are minimizers of the discrete energy. As an initial condition for both the discrete and continuum simulations, we assume that the rigid and deformable chains of atoms are parallel and at a distance $\sigma$ apart.

Figure 5 shows the result of an atomistic simulation for which $\sigma=1.0$, $\omega=1.0$, $h_{1}=1.0$, $h_{2}=.99$, $k_{b}=100.0$, and $k_{s}=10.0$. There are two relatively large regions on which the atoms on the deformable chain are uniformly spaced and the chain is parallel to the rigid chain. These two flat regions are separated by a single narrow region with a relatively large vertical displacement, or wrinkle. The atoms on the rigid chain are at positions $(j,0)$ with $j$ an integer. Hence the van der Waals interaction from the rigid chain creates potential wells located above every point $j+1/2$ on the $x$-axis.

As seen in inset (a) in Figure 5, in the regions where the deformable chain is parallel to the rigid chain, the atoms fall into the potential wells of the van der Waals interaction. Hence in these regions the atoms are 1 unit apart. However, because $h_{2}=.99$, over the length of the system there is one ‘extra’ atom on the chain that has no potential well. As seen in inset (b) in Figure 5, this extra atom is accommodated by a wrinkle, which minimizes the cost of placing $n$ atoms over $n-1$ potential wells. This wrinkle forms the domain wall between the large commensurate regions to the left and the right.

One way to understand the configuration of wrinkles that forms in a system of two incommensurate chains is to consider the corresponding moiré pattern. As noted in the introduction, when two lattices with different geometries or the same geometry but different orientations are stacked, a larger periodic pattern, called a moiré pattern, emerges [2, 3]. A similar, one-dimensional pattern arises when two incommensurate periodic chains of atoms are overlaid. In Figure 5 this pattern is shown below the plot of the discrete solution. The atoms of the rigid and deformable chains in the reference configuration are both represented by circles. At the left endpoint, the atoms of the deformable chain lie exactly above the midpoint between the atoms of the rigid chain, so that the corresponding circles show the minimum overlap. As one moves to the right along the chains, because of the difference between the lattice parameters $h_{1}$ and $h_{2}$, the atoms of the deformable chain progressively shift in relation to the atoms of the rigid chain until the circles lie exactly on the top of each other. This creates a pattern of lighter and darker regions, where the darker regions correspond to optimal atomic registry. On the other hand, the atoms of the deformable chain are out of registry in the lighter regions, and this is where the wrinkles form. Note that each light region corresponds to a single extra atom on the deformable chain, hence the number and locations of light regions in the moiré pattern should predict the number and the locations of wrinkles. For instance, in Figure 5, there is exactly one extra atom on the deformable chain and, as a consequence, exactly one light region and one wrinkle.

Figures 6 and 7 show the numerical solutions of both discrete and continuum models for parameter values

[TABLE]

and for two different values of $h_{2}$. In both cases, we see that the continuum solution has the same structure as the discrete solution, with a single wrinkle separating two commensurate regions. The continuum solution is close to the discrete solution except that the amplitude of the wrinkle in the continuum solution is significantly smaller than that in the discrete solution.

To further understand the details within the wrinkle, we note that when $h_{1}<h_{2}$, there are fewer atoms per period on the deformable chain than there are atoms on the rigid chain. Consequently, the deformable chain contracts between the wrinkles so that interatomic distance in the deformable chain matches that in the rigid chain. The contraction in the bulk causes the expansion of distances between atoms of the deformable chain inside the wrinkles thus making the wrinkles relatively wide (the left plot in Figure 7). On the other hand, when $h_{1}>h_{2}$, there are more atoms per period on the deformable chain, leading to bulk expansion and wrinkles that are relatively narrow (the right plot in Figure 7).

Note that only three atoms form the body of the wrinkle in the discrete solution in the two cases shown in Figures 6 and 7. This is to be expected since the nondimensional width of the wrinkle and the nondimensional interatomic distance are of order $\varepsilon$ (these are of order $1$ in dimensional units). With so few atoms, there is no reason to expect that the continuum approximation should remain valid inside the wrinkles. The problem is further compounded by large gradients of the continuum solution in these singular regions. Indeed, as Figure 7 demonstrates, the match between the discrete and continuum solutions is reasonable away from the center atom inside the wrinkle, while there is a significant discrepancy near the center atom. A further problem, associated with contraction inside the wrinkle, can be seen in the plot on the right in Figure 7. The continuum solution forms a narrow loop near the maximum, indicating self-penetration of the deformable chain. Hence, the continuum solution for the combination of parameters (31) is unphysical near the top of the wrinkle. Note, however, that the nonphysicality is restricted to the region where there is only a single atom of the discrete solution. Indeed, as Figure 8 demonstrates, the continuum displacements predict the discrete displacements well at the atomic positions for the values (31) of the parameters used to produce the plot on the right in Figure 7. We emphasize that the accuracy of the continuum approximation is the same for both plots in Figure 7. There is no self-penetration in the plot on the left simply because near the wrinkle the chain experiences local expansion rather than contraction.

The elastic constants $k_{s}$ and $k_{b}$ are the same in (31), so that the cost of stretching and bending is comparable for the solutions depicted in Figures 6 and 7. That bending is relatively expensive tends to spread bending deformation, causing the dips that form on both sides of the wrinkles as well as depressing the amplitude of the wrinkles themselves. In particular, large bending stiffness causes the deformable sheet to “undershoot” the equilibrium distance between the chains in registry at the edge of a wrinkle. A subsequent “bounce-off” away from the substrate is the cause of the dips adjacent to the wrinkles. A macroscopic version of this effect is well known for elastic beams on a liquid [9]. An analysis similar to the one presented in [9] can be performed here to study the shape of the dips. Indeed, the ability to analyze the system of ordinary differential equations to establish the properties of equilibrium solutions is one of the clear advantages of working with the continuum rather than the discrete model.

We now fix the strength of the Lennard-Jones interaction $\omega=1$ and investigate the influence of the elastic constants on the discrete and continuum solutions. Generically speaking, we expect that with larger elastic constants the solution should incur a larger penalty for elastic deformation and, consequently, should exhibit smaller gradients. This in turn should lead to wider wrinkles that involve more atoms per wrinkle. Both effects should have a positive influence on the accuracy and applicability of the continuum approximation near the wrinkles.

In Figure 9, we set $k_{s}=100.0$ and $k_{b}=1.0$ so that stretching deformation is significantly more expensive than bending. As a result, the bending stiffness is small and the dips next to the wrinkles disappear, both in the case of local stretching (Figure 9, left) and local compression (Figure 9, right). As expected, there are more atoms involved in forming a wrinkle, and the continuum approximation matches the discrete solution very well. Further, the amplitude of the wrinkles is larger than that observed in Figure 6. Note that, even though the number of atoms that form the wrinkle is now larger, the wrinkle itself still serves to accommodate only one extra atom on the deformable chain.

In Figure 10, we consider the opposite case, $k_{s}=1.0$ and $k_{b}=100.0$, so that bending is significantly more expensive than stretching. Because the bending stiffness is now large, the dips next to the wrinkles are enhanced and the wrinkles have a smaller amplitude than in Figure 6. The combination of large bending stiffness and compression causes the self-penetration near the top of the wrinkle in the plot on the right in Figure 10. Note that, typically, bending in two-dimensional materials is much cheaper than in-plane compression [18, 19], so that in practice one is unlikely to encounter the extreme situation depicted in the plot on the right in Figure 10.

The solutions for smaller value of the parameter $\varepsilon$—when the size of the system is assumed to be larger given the same interatomic distance—are shown in Figures 11–13. For the relatively large stretching constant, $k_{s}=10.0$, it appears that all continuum solutions match well their discrete counterparts.

If the deformable chain has three extra atoms over the length of the system, then three equispaced wrinkles form as shown in Figure 14. The locations and the number of wrinkles are predicted by the moiré pattern at the bottom of Figure 14.

The amplitude of the wrinkles is determined by how close an atom on the deformable lattice wants to be to the rigid chain when it is in registry (i.e, it is positioned above the midpoint between two atoms of the rigid chain) versus when it is out of registry (i.e., it is positioned exactly above one of the atoms of the rigid chain). As we can see in Figure 15, the equilibrium distance for an atom in registry and out of registry are, respectively, $\sim 0.86$ and $\sim 0.98$. The lower of these values corresponds to the distance between the chains when they are in registry in Figures 6–14, while the higher value closely approximates the height of the wrinkles observed in the same figures, unless the amplitude is depressed due to a relatively large bending stiffness.

We conclude this section by emphasizing that the wrinkling predicted by the continuum model occurs because the interaction energy $G$ retains information about the discrete character of the substrate. Recall that, in the nondimensional variables, $h_{1}$ is the spacing between atoms on the rigid chain and $\sigma$ is the equilibrium separation of the Lennard-Jones atom-to-atom potential. Because of the fast decay of this potential, the equilibrium separation essentially determines the range of the van der Waals interaction between atoms. If $h_{1}$ is close to $\sigma$, an atom on the deformable chain interacts with only the one or two nearby atoms, both to the left and to the right on the rigid chain. The strength of the interaction varies significantly in the regions above the ‘gaps’ between atoms on the rigid chain, which creates deep potential wells. As explained above, these potential wells along with the mismatch drive the formation of isolated wrinkles. Conversely, when $h_{1}$ is much smaller then $\sigma$, an atom on the deformable chain interacts with a relatively large number of nearby atoms on the rigid chain and the gaps between these nearby atoms are relatively small. Hence the atom essentially feels the average of this interaction, which is like a continuum approximation in which the strength of the interaction depends only on the vertical distance from the line of fixed atoms. One would not expect the formation of isolated wrinkles in this case.

To support these observations, we note that for the solutions shown in Figures 6—14, $h_{1}=\sigma=1.0$. We contrast the isolated wrinkles observed in those cases with the discrete and continuum solutions shown in Figure 16.

For the discrete solution in Figure 16, $h_{1}=1.01$ and $\sigma=3.0$ and the corresponding configuration is approximately sinusoidal. Furthermore, the amplitude of the solution is on the order of $10^{-4}$, while the height of the wrinkles in Figures 6—14 is on the order of $10^{-2}$. The continuum solution for parameters corresponding to those used in the discrete simulation shows an excellent match with the discrete solution.

5 Conclusions

We have applied an upscaling procedure to develop a mesoscopic continuum model of the cross-section of a graphene sheet interacting with a rigid crystalline substrate. Without making any a priori assumptions on the structure of the continuum model, we established that the elastic part of the model resembles the Föppl–von Kármán theory of thin rods but without thickness as that notion is meaningless for the 2D materials (or their 1D reduction). Our approach allows us to identify the potential function that keeps track of microscopic registry effects on the mesoscale and also determines the correct orders of magnitude of various contributions to the continuum energy in terms of the small parameter $\varepsilon$. The main novelty of this work is in the development of the continuum contribution for the weak interaction between the deformable and rigid chains. Although a continuum description, this energy retains discrete information about any mismatch between the chains.

Numerical simulations demonstrate that the predictions of the mesoscopic model are generally in close correspondence with the configurations obtained via a discrete molecular dynamics approach. In both cases, relaxation of slightly mismatched chains resulted in formation of wrinkles, where flat, approximately uniform domains, are separated by wrinkles characterized by relatively large out-of-plane deformations. Fixing the strength of the Lennard-Jones potential, and using a flat, undeformed configuration as the initial condition, we considered a number of different parameter regimes. When the elastic constants corresponding to stretching and bending are of order $1$, the wrinkles contain a small number of atoms and the quality of the continuum approximation near the top of the wrinkle is poor, while the approximation is accurate both near the base of the wrinkle and in the flat regions of registry. In the case when the stretching constant is of order $>10$ and is larger than the bending constant, which is of practical importance, a solution of the continuum model approximates well the corresponding solution of the discrete model, even near the top of the wrinkle. This improved approximation occurs because there are more atoms per wrinkle and because the solution exhibits relatively smaller gradients for a given $\varepsilon$.

Compared to the number of atoms on the rigid chain, there can be more or fewer atoms on the deformable chain. The number of wrinkles is then equal to the difference between the number of atoms on the chains, with each wrinkle compensating for exactly one extra atom or vacancy. The locations of the wrinkles correspond to the regions of the largest misfit in the moiré pattern for a given system of incommensurate deformable and rigid chains.

Unlike what is seen in analogous thin film/substrate systems with misfit strains, here the wrinkling is observed not only when the deformable chain is compressed but also when it is stretched. This is true because the equilibrium distance between an atom of the deformable chain and the rigid chain is strictly larger in the region of misregistry than in the region of registry. The difference between these equilibrium distances determines the amplitude of a wrinkle.

Periodic boundary conditions were imposed in this work primarily for simplicity, and other boundary conditions can be considered. The model can also be extended in a standard way to include applied body forces. In this more general setting it would be possible to consider a problem in which the registry effects can be combined with a standard compressive misfit between the deformable and rigid chains. It would be of interest to understand how the wrinkling pattern would respond to the applied loads and whether different types of singular regions—e.g., wrinkles and delamination blisters—may coexist in the same system.

Continuum modeling that retains discrete registry effects is important both for solving computational problems more efficiently and for allowing theoretical insight into mesoscopic pattern formation in graphene and other 2D materials. Here one of the principal advantages is the ability to fully utilize the powerful tools provided by the theory of differential equations and variational calculus. The similarity between the structure of the continuum model that was derived in this paper and the existing models of a thin film on an adhesive substrate should allow for an extension of the known results for these macroscopic systems to 2D materials. For example, the formation of dips in the vicinity of the wrinkles, which was discussed in the previous section, may likely be explained using the techniques developed in [9]. Continuum modeling may also facilitate the study of the influence of atomic relaxation of slightly mismatched graphene/substrate lattices on electronic properties of the system [2].

In addition to the analysis of the one-dimensional model, future work will include a rigorous verification of our conjecture that both the discrete and continuum models converge to the same asymptotic limit in the appropriate sense as $\varepsilon\to 0$. Here the limit has to be understood within the framework of $\Gamma$-convergence [16]. This analysis should justify the number of terms that was retained in the expansion of the discrete problem in order to obtain the continuum problem. Further, we plan to develop an extension of the derivation presented in this paper to a full model for a 2D material consisting of two slightly mismatched, interacting lattices.

Acknowledgment: This work was supported by the National Science Foundation grant DMS-1615952.