Elastic properties and mechanical stability of bilayer graphene: Molecular dynamics simulations

Graphene has become in last decades a paradigmatic example of two-dimensional and so-called van-der-Waals layered materials, showing large anisotropy in their physical properties. Here we study the elastic properties and mechanical stability of graphene bilayers in a wide temperature range by molecular dynamics simulations. We concentrate on in-plane elastic constants and compression modulus, as well as on the atomic motion in the out-of-plane direction. Special emphasis is placed upon the influence of anharmonicity of the vibrational modes on the physical properties of bilayer graphene. We consider the excess area appearing in the presence of ripples in graphene sheets at finite temperatures. The in-plane compression modulus of bilayer graphene is found to decrease for rising temperature, and results to be higher than for monolayer graphene. We analyze the mechanical instability of the bilayer caused by an in-plane compressive stress. This defines a spinodal pressure for the metastability limit of the material, which depends on the system size. Finite-size effects are described by power laws for the out-of-plane mean-square fluctuation, compression modulus, and spinodal pressure. Further insight into the significance of our results for bilayer graphene is gained from a comparison with data for monolayer graphene and graphite.


I. INTRODUCTION
Over the last few decades there has been a surge of interest in carbon-based materials with sp 2 orbital hybridization, such as fullerenes, carbon nanotubes, and graphene [1][2][3], continuously enlarging this research field beyond the long-known graphite.In particular, bilayer graphene displays peculiar electronic properties, which have been discovered and thoroughly studied in recent years [4,5].It presents unconventional superconductivity for stacking of the sheets twisted relative to each other by a precise small angle [6,7].Such rotated graphene bilayers show magnetic properties that may be controlled by an applied bias voltage [8,9].Also, localized electrons are present in the superlattice appearing in a moiré pattern, so that one may have a correlated insulator [10].Bilayer graphene displays ripples and out-of-plane deformations akin to suspended monolayers [11], thus giving rise to a lack of planarity which may be important for electron scattering [12].
A deep comprehension of thermodynamic properties of two-dimensional (2D) systems has been a challenge in statistical physics for many years [13,14].This question has been mainly discussed in the field of biological membranes and soft condensed matter [15,16], for which analyses based on models with realistic interatomic interactions are hardly accessible.In this context, graphene is a prototype crystalline membrane, appropriate to study the thermodynamic stability of 2D materials.This problem has been addressed in connection with anharmonic effects, in particular with the coupling between in-plane and out-of-plane vibrational modes [17,18].Bilayer graphene is a well-defined two-sheet crystalline membrane, where an atomic-level characterization is feasible, thereby permitting one to gain insight into the physical properties of this type of systems [17,[19][20][21][22].
Mechanical properties of graphene, including elastic constants, have been studied by using several theoretical [23][24][25][26] and experimental [27][28][29][30][31][32] techniques.These methods have been applied to analyze monolayer as well as multilayer graphene, including the bilayer [25,[33][34][35][36].In this context, a theory of the evolution of phonon spectra and elastic constants from graphene to graphite was presented by Michel and Verberck [23].In particular, for bilayer graphene on SiC, Gao et al. [32] have found a transverse stiffness and hardness comparable to diamond.More generally, mechanical properties of graphene and its derivatives have been reviewed by Cao et al. [37], and various effects of strain in this material were reported by Amorim et al. [38].
In this paper we extend earlier work on isolated graphene sheets to the bilayer, for which new aspects show up due to interlayer interactions and the concomitant coupling between atomic displacements in the outof-plane direction.We use molecular dynamics (MD) simulations to study structural and elastic properties of bilayer graphene at temperatures up to 1200 K. Especial emphasis is laid on the behavior of bilayer graphene under tensile in-plane stress and on its mechanical stability under compressive stress.MD simulations allow us to approach the spinodal line in the phase diagram of bilayer graphene, which defines its stability limit.We compare results found for the bilayer with data corresponding to monolayer graphene and graphite, which yields informa-tion on the evolution of physical properties from an individual sheet to the bulk.
The paper is organized as follows.In Sec.II we describe the method employed in the MD simulations.In Sec.III we present the phonon dispersion bands and the elastic constants at T = 0.In Sec.IV we present results for structural properties derived from the simulations: interatomic distances, interlayer spacing, and outof-plane atomic displacements.The in-plane and excess area are discussed in Sec.V, and in Sec.VI we analyze the elastic constants and compressibility at finite temperatures, along with the stability limit for compressive stress.Finite-size effects are studied in Sec.VII.The papers closes with a summary of the main results in Sec.VIII.

II. METHOD OF CALCULATION
In this paper we employ MD simulations to study structural and elastic properties of graphene bilayers as functions of temperature and in-plane stress.The interatomic interactions in graphene are described with a long-range carbon bond-order potential, the socalled LCBOPII [44], used earlier to perform simulations of carbon-based systems, such as graphite [44], diamond, [44] and liquid carbon [45].In more recent years, this interatomic potential has been utilized to study graphene [19,26,39,46], and in particular mechanical properties of this 2D material [47,48].The LCBOPII potential model was also used to conduct quantum pathintegral MD simulations of graphene monolayers [49] and bilayers [20], which allowed to assess nuclear quantum effects in various properties of this material.Here, as in earlier simulations [46,48,49], the original LCBOPII parameterization has been slightly modified in order to rise the bending constant κ of a graphene monolayer from 0.82 eV to a value of 1.49 eV, close to experimental results and ab-initio calculations [50].Values of the parameters employed here for the torsion term of the potential are given in Appendix A.1.
For the interlayer interaction we have considered the same parameterization as that previously used in simulations of graphene bilayers with this potential model [19,20], presented in Appendix A.2.For the minimumenergy configuration of bilayer graphene with AB stacking, we find an interlayer binding energy of 25 meV/atom (50 meV/atom for graphite) [19].
Our simulations were carried out in the isothermalisobaric ensemble, where one fixes the number of carbon atoms, 2N (i.e., N atoms per sheet), the in-plane stress tensor, {τ ij }, and the temperature, T .We have considered rectangular supercells with similar side lengths in the x and y directions in the layer plane, L x ≈ L y .These supercells included from N = 48 to 8400 carbon atoms per graphene sheet.Periodic boundary conditions were assumed for x and y coordinates, whereas C atoms were allowed to move freely in the out-of-plane z coordinate (free boundary conditions).
To keep a given temperature T , chains of four Nosé-Hoover thermostats were connected to each atomic degree of freedom [51].An additional chain including four thermostats was connected to the barostat which regulates the in-plane area of the simulation cell (xy plane), keeping the required stress {τ ij } [51,52].Integration of the equations of motion was performed by using the reversible reference system propagator algorithm (RESPA), which permits to consider different time steps for slow and fast degrees of freedom [53].For the atomic dynamics derived from the LCBOPII potential, we took a time step ∆t = 1 fs, which gave good accuracy for the temperatures and stresses discussed here.For fast dynamical variables such as the thermostats, we used δt = ∆t/4.
The configuration space has been sampled for T in the range from 50 to 1200 K. Given a temperature, a typical run consisted of 2×10 5 MD steps for system equilibration and 8 × 10 6 steps for calculation of ensemble averages.In Fig. 1 we show top and side views of an atomic configuration of bilayer graphene obtained in MD simulations at T = 800 K.In the top view, red and yellow balls stand for carbon atoms in the upper and lower graphene layers in AB stacking pattern.In the side view, one can see atomic displacements in the out-of-plane direction, clearly observable at this temperature.
To characterize the elastic properties of bilayer graphene we consider uniaxial stress along the x or y directions, i.e., τ xx = 0 or τ yy = 0, as well as 2D hydrostatic pressure P (biaxial stress) [54], which corresponds to τ xx = τ yy = −P , τ xy = 0. Note that P > 0 and P < 0 mean compressive and tensile stress, respectively.
For comparison with our results for graphene bilayers, we also performed some MD simulations of graphite using the same potential LCBOPII.In this case we considered cells containing 4N carbon atoms (four graphene sheets), and periodic boundary conditions were assumed in the three space directions.
Other interatomic potentials have been used in last years to analyze several properties of graphene, in particular the so-called AIREBO potential [55][56][57][58][59].Both LCBOPII and AIREBO models give very similar values for the equilibrium C-C interatomic distance and for the thermal expansion coefficient [49,56,59].For the Young's modulus of graphene, we find that the result obtained by employing the LCBOPII potential is closer to those yielded by ab initio calculations [56].

III. PHONON DISPERSION BANDS AND ELASTIC CONSTANTS AT T = 0
The elastic stiffness constants, c ij , of bilayer graphene calculated with the LCBOPII potential model in the limit T → 0 can be used as reference values for the finitetemperature analysis presented below.These elastic constants are calculated here from the harmonic dispersion relation of acoustic phonons.The interatomic force constants utilized to obtain the dynamical matrix were obtained by numerical differentiation of the forces using atom displacements of 1.5 × 10 −3 Å with respect to the equilibrium (minimum-energy) positions.
The phonon dispersion of bilayer graphene, calculated by diagonalization of the dynamical matrix is presented in Fig. 2 along high-symmetry directions of the 2D Brillouin zone.One finds 12 phonon bands, corresponding to four C atoms (2 per layer) in the crystallographic unit cell.Labels indicate the common names of the phonon bands: eight branches with in-plane atomic motion (LA, TA, LO, and TO, all of them two-fold degenerate), and four branches with displacements along the out-of-plane direction (ZA, ZO', and a two-fold degenerate ZO band).The phonon dispersion presented in Fig. 2 is analogous to those obtained for other empirical potentials and DFT calculations [60][61][62][63].We emphasize the presence of the flexural ZA band, which is parabolic close to the Γ point (ω ∼ k 2 ), and typically appears in 2D materials [64][65][66][67].Here k denotes the wavenumber, i.e., k = |k|, and k = (k x , k y ) is a wavevector in the 2D hexagonal Brillouin zone.
Note also the presence of the optical mode ZO', which does not appear in monolayer graphene, and in the case of the bilayer corresponds to the layer-breathing Ramanactive A 2g mode, for which a frequency of 89 cm −1 has been measured [68].The LCBOPII potential yields for this band at the Γ point (k = 0) a frequency of 92 cm −1 .This value is close to that found from ab initio calculations for graphene bilayers [60].
The interatomic potential LCBOPII was used before to calculate the phonon dispersion bands of graphene and graphite [61].However, the version of the potential employed in Ref. [61] was somewhat different from that considered here, which gives a description of the graphene bending closer to experimental results [46,50] (see Appendix A.1).
The sound velocities for the acoustic bands LA and TA along the direction ΓM, with wavevectors (k x , 0, 0), are given by the slope (∂ω/∂k x ) Γ in the limit k x → 0. The elastic stiffness constants can be obtained from these velocities by using the following expressions, valid for the hexagonal symmetry of graphene [69]: where ρ is the surface mass density of graphene.We find c 11 = 20.94eV Å−2 and c 12 = 4.54 eV Å−2 .Note that the dimensions of these elastic constants (force/length) coincide with those of the in-plane stress.These c ij can be converted into elastic constants C ij (units of force per square length), typical of three-dimensional (3D) materials, as the minimum-energy configuration of bilayer graphene.Taking d 0 = 3.3372 Å, we find for the bilayer C 11 = 1005 GPa and C 12 = 218 GPa, near the values found for graphite in the classical low-T limit, using the LCBOPII potential: 1007 and 216 GPa, respectively [67].

IV. STRUCTURAL PROPERTIES A. Interatomic distance
For bilayer graphene, the minimum-energy configuration for the LCBOPII potential corresponds to planar sheets and the interatomic distance between nearest neighbors in a layer amounts to 1.4193 Å.This distance turns out to be a little smaller than that found for monolayer graphene using the same interatomic potential (r 0 = 1.4199Å).This fact was noticed by Zakharchenko et al. [19] in their results of Monte Carlo simulations of graphene bilayers.
We have studied the change of the interatomic C-C distance r (actual distance in 3D space) as a function of 2D hydrostatic pressure, P , at several temperatures.In Fig. 3 we present the dependence of r on P at T = 300 and 1000 K.For T = 300 K, data derived from MD simulations are shown for three cell sizes : N = 240 (solid circles), 448 (open squares), and 3840 (open diamonds).We observe that the size effect on the interatomic distance is negligible, since differences between the results for different cell sizes are much smaller than the symbol size in Fig. 3.The data for T = 1000 K (open triangles) correspond to N = 240.Note that 2D hydrostatic pressure P > 0 corresponds to compressive stress.
Close to P = 0, this dependence can be fitted for T = 300 K to an expression r = r m + µP , where r m = 1.4230Å is the interatomic distance for the stress-free bilayer at this temperature and µ = −0.0289Å3 /eV.For T = 1000 K, we find r m = 1.4314Å and µ = −0.0302Å3 /eV.This slope is slightly larger than that found for T = 300 K.
In connection with the interatomic distance r, we note that for a strictly planar geometry, the area per atom for an ideal honeycomb pattern is given by S p = 3 √ 3 r 2 /4.At finite temperatures, however, the graphene layers are not totally planar and the actual in-plane area per atom is smaller than that given by the above expression, using the mean interatomic distance between nearest-neighbor C atoms.This is related with the so-called excess area and is discussed below in Sec.V.
In each hexagonal ring, two C-C bonds are aligned parallel to the y direction (vertical in Fig. 1, top image), and the four other bonds form an angle of 30 degrees with the x axis (horizontal direction).A compressive uniaxial stress along the y axis (τ yy < 0) causes a decrease in the length of the former bonds and an increase in the latter, as corresponds to a positive Poisson's ratio.The opposite happens for τ xx < 0.

B. Interlayer spacing
For the minimum-energy configuration we find an interlayer distance d 0 = 3.3372 Å, to be compared with a distance of 3.3371 Å obtained in Ref. [19] using an earlier version of the LCBOPII potential.At T = 300 K we obtain d = 3.374 Å, i.e., the interlayer distance increases somewhat due to bending of the graphene sheets caused by thermal motion.
The interlayer spacing is reduced in the presence of a tensile 2D hydrostatic pressure.Thus, for P = −0.5 eV/ Å2 and T = 300 K, we find d = 3.367 Å.This decrease is due to a reduction in the out-of-plane fluctuations under a tensile stress.The effect of this relatively high stress on the distance d is, however, smaller than the thermal expansion up to 300 K.In the presence of a compressive stress in the xy plane, one has an expansion of the interlayer spacing, but this kind of stress causes an instability of the bilayer configuration for relatively small values of P , as explained below.
The mean-square fluctuation (MSF) of the interlayer distance, (∆d) 2 = d 2 − d 2 , associated to thermal motion at finite temperatures, is related to the interaction between graphene layers.This MSF is expected to depend on the size of the considered simulation cell, and becomes negligible in the thermodynamic limit (N → ∞).In Fig. 4 we display (∆d) 2 derived from our MD simulations as a function of the inverse cell size for stressfree bilayer graphene at T = 100 K (circles) and 300 K (squares).One observes that (∆d) 2 → 0 for 1/N → 0, and grows linearly for increasing inverse cell size.
To connect these results with the energetics of bilayer graphene, we have calculated the interlayer interaction energy for several values of the spacing d near the distance d 0 corresponding to the minimum-energy configuration.The interaction energy per atom can be written as: where E 0 int is the energy for distance d 0 , and k is an effective interaction constant which is found to amount to 0.093 eV/ Å2 .Then, for the whole simulation cell (2N atoms in a bilayer), the energy corresponding to a distance d close to d 0 is Thus, thermal motion at temperature T , associated to the degree of freedom d, will cause a MSF of this variable, (∆d) 2 , given by the mean potential energy: where k B is Boltzmann's constant.This means that for a given temperature, (∆d) 2 scales as 1/N , as shown in Fig. 4.
As indicated in Sec.III, the phonon spectra of monolayer and bilayer graphene are similar.The main difference between them is the appearance in the latter of the ZO' vibrational band, which is almost flat in a region of k-space near the Γ point (see Fig. 2).As noted above, this vibrational mode of bilayer graphene is the layerbreathing Raman-active A 2g mode [68].The frequency of the ZO' band at Γ (which will be denoted here ω 0 ) can be related to the interlayer coupling constant k as ω 0 = (k N /M red ) 1/2 , with k N = 2N k and the reduced mass M red = N m/2 (m: atomic mass of carbon).We find ω 0 = 2(k/m) 1/2 .and putting for the coupling constant k = 0.093 eV Å−2 /atom, one obtains ω 0 = 92 cm −1 , which coincides with the frequency of the ZO' band derived from the dynamical matrix at the Γ point.Michel and Verberck [23] have studied the evolution of this frequency ω 0 with the number of sheets n in graphene multilayers.They found an increase of ω 0 for rising n, which saturates to a value of 127 cm −1 for large n (graphite).
The interlayer coupling in bilayer graphene was studied before by Zakharchenko et al. [19] by means of Monte Carlo simulations.The low-frequency part of the ZO' band was described by these authors using a parameter γ, which is related to the parameters employed here as γ = ρ ω 2 0 /4, where ρ is the surface mass density.From this expression we find γ = 0.035 eV Å−4 , which agrees with the low-temperature result derived from Monte Carlo simulations (Fig. 7 in Ref. [19]).
Fluctuations in the interlayer spacing of bilayer graphene at temperature T are related with the isothermal compressibility in the out-of-plane direction, χ z [20].In fact, χ z can be calculated from the MSF (∆d) 2 by using the expression [20] Using (∆d) 2 = 1.65 × 10 −4 Å2 for N = 960 at T = 300 K, we find χ z = 2.96 × 10 −2 GPa −1 .This value is a little larger than experimental results for graphite, of about 2.7 × 10 −2 GPa −1 [70,71].This is consistent with the fact that bilayer graphene is more compressible than graphite in the z direction, since in the latter each layer is surrounded by two other graphene layers, whereas in the former each sheet has a single neighbor.

C. Out-of-plane motion
The minimum-energy configuration for the graphene layers, i.e., the classical low-temperature limit, corresponds to planar sheets.At finite temperatures the graphene sheets are bent.This bending is directly related to the atomic motion in the out-of-plane z direction, whose largest vibrational amplitudes come from low-frequency ZA modes with long wavelength (small k).For stress-free graphene, the ZA phonon branch can be described close to the Γ point by a parabolic dispersion relation of the form ω(k) = κ/ρ k 2 , with the bending constant κ = 1.49eV (see Fig. 2).
It is known that the atomic MSF in the xy layer plane is relatively insensitive to the system size, but the out-ofplane MSF has important finite size effects.This dependence on N has been studied earlier for stress-free monolayer and bilayer graphene by means of Monte Carlo and molecular dynamics simulations, in particular using the LCBOPII interatomic potential [49].For system size N , one has an effective cut-off for the wavelength λ given by λ max ≈ L, where L = (N S p ) 1/2 , and S p is the inplane area per atom.Thus, the minimum wavenumber present for size N is k min = 2π/λ max , and the minimum frequency for ZA modes is so that ω min scales as N −1 .For a graphene sheet, we call r ≡ (x, y) the 2D position on the layer plane and h(r) is the distance to the mean plane of the sheet.In Fig. 5 we present the MSF of the atomic positions in the z direction for bilayer graphene, (∆h) 2 = h 2 − h 2 , as a function of 2D hydrostatic pressure P for various cell sizes at T = 300 K. Symbols represent results derived from our MD simulations, with cell size decreasing from left to right: N = 3840, 960, 448, 308, and 240.One observes first that (∆h) 2 appreciably increases for rising system size, as expected from earlier studies of 2D materials [72].For the largest size displayed in Fig. 5, N = 3840, we find at P = 0, (∆h) 2 = 0.22±0.01Å2 (not shown in the figure).The dependence of (∆h) 2 on N for stress-free bilayer graphene will be analyzed below in Sec.VII.Second, we also observe in Fig. 5 that the difference in atomic MSF between different system sizes is reduced for increasing tensile stress (P < 0).(∆h) 2 grows as the tensile stress is reduced (d(∆h) 2 /dP > 0), and eventually diverges at a size-dependent critical pressure P c (N ) > 0. Third, one sees that P c (N ) approaches zero for rising system size.
The dependence of the critical pressure P c on N will be discussed below in relation to fluctuations of the inplane area S p , which are also found to diverge in parallel with (∆h) 2 for each size N .The origin of this instability is related to the appearance of imaginary frequencies for vibrational modes in the ZA flexural band for pressure P c (N ).This will be discussed in Sec.VI in connec- tion with the 2D modulus of hydrostatic compression B p , which is found to vanish at P c .

V. IN-PLANE AND EXCESS AREA
The in-plane area, S p = L x L y /N , is the variable conjugate to the pressure P in the isothermal-isobaric ensemble considered here.Its temperature dependence, S p (T ), has been analyzed earlier in detail for monolayer and bilayer graphene from atomistic simulations [19,20,49].For the bilayer we find in the minimum-energy configuration (T = 0) an area S 0 = 2.6169 Å2 /atom, in agreement with earlier calculations [19,20].Here, we discuss the behavior of S p as a function of 2D hydrostatic stress, both tensile and compressive.
In Fig. 6 we display the dependence of S p on P for several cell sizes at T = 300 K. We present data for N = 308, 448, 960, and 8400.For tensile stress P < −0.05 eV/ Å2 , S p data for different cell sizes are indistinguishable at the scale of Fig. 6.In fact, we obtain a nearly linear dependence with a slope dS p /dP ≈ −0.12 Å4 /eV.However, differences appear close to P = 0, and even more for compressive stress (P > 0).For each size N , one observes a fast decrease in S p close to the corresponding stability limit of the planar phase.We obtain values of the in-plane area below 2.58 Å2 /atom, not shown in the figure for clarity.Changes in S p correspond to linear strain ǫ L as: S p = S 0 (1 + ǫ L ) 2 .This means that the vertical range in Fig. 6 corresponds to a strain range between ǫ L = −7.1 × 10 −3 (compression) and 4.4 × 10 −3 (tension).
In our MD simulations, carbon atoms are free to move in the out-of-plane direction (z coordinate), and the real surface of a graphene sheet is not strictly planar, hav-ing an actual area larger than the area of the simulation cell in the xy plane.Differences between the in-plane area S p and real area S r were considered earlier in the context of biological membranes [15,73,74] and in recent years for graphene, as a paradigmatic crystalline membrane [31,48].An explicit differentiation between both areas is relevant to understand certain properties of 2D materials [13].Some experimental techniques can be sensitive to properties connected to the area S r , whereas other methods may be suitable to analyze variables related to the area S p [31,75].
Here we calculate the real area S r of both sheets in bilayer graphene by a triangulation method based on the atomic positions along the simulation runs [20,48].In the sequel, S r will denote the real area per atom.The difference S r − S p has been called in the literature hidden area [31] or excess area [76,77].We consider the dimensionless excess area Φ for a graphene sheet, defined as [76,77] In the classical low-temperature limit, Φ vanishes, as the sheets become strictly planar for T → 0. We note that this is not the case in a quantum calculation, where one has Φ > 0 for T → 0 due to atomic zero-point motion [78].
The excess area is related to the amplitude of the vibrational modes in the z direction.This allows one to find analytical expressions for Φ in terms of the frequency of those modes.The instantaneous real area S * r may be expressed in a continuum approximation as [48,73,74] where h(r) represents the distance to the mean xy plane of the sheet, as in Sec.IV.C. Expanding h(r) as a Fourier series with wavevectors k = (k x , k y ) in the 2D hexagonal Brillouin zone, the real area S r = S * r may written as [13,15,48] where H(k) are the Fourier components of h(r) (see Appendix B).Taking into account that the MSF of a mode with frequency ω j (k) is given by k B T /mω j (k) 2 , m being the atomic mass of carbon, one finds for the excess area where the sum in j is extended to phonon branches with atomic motion in the z direction, i.e., ZA, ZO, and ZO'.For small k, the contribution of ZO and ZO' modes to the sum in Eq. ( 11) vanishes for k → 0, as in both cases ω j (k) converges to positive values.For the flexural ZA band with negligible effective stress, we have ω ZA (k) ∝ k 2 , and k 2 /ω ZA (k) 2 ∝ k −2 , so that the contribution of ZA modes with small k is dominant in the sum in Eq. (11).The minimum wavenumber k min available for cell size N scales as k min ∼ N −1/2 (see Sec. IV.C).Thus, its contribution to Φ grows linearly with N , and diverges for stress-free graphene in the thermodynamic limit.This divergence disappears in the presence of a tensile in-plane stress (even small), be it caused by internal tension or by an external pressure.
In Fig. 7 we display Φ for bilayer graphene as a function of 2D hydrostatic pressure P .Symbols represent data derived from our MD simulations at three temperatures: T = 300 K (circles), 600 K (squares), and 1000 K (diamonds).Dashed lines are guides to the eye.These data were obtained for system size N = 960.The excess area Φ increases as T is raised, in agreement with the growing amplitude of the out-of-plane vibrational modes (see Eq. ( 11)).In fact, a classical harmonic approximation (HA) for the vibrational modes predicts a linear increase of Φ with temperature.From the results shown in Fig. 7 for P = 0, we find the ratios Φ(1000 K)/Φ(300 K) = 3.14 and Φ(600 K)/Φ(300 K) = 1.93, a little less than the corresponding temperature ratios (3.33 and 2.0).For a pressure P = −0.5 eV/ Å2 , we find for those ratios at the same temperatures values of 3.28 and 1.99, respectively, closer to the harmonic expectancy.An in-plane tensile stress causes a decrease in the vibrational amplitudes in the z direction, and then the modes are better described by a harmonic approach.

VI. ELASTIC CONSTANTS AND COMPRESSIBILITY AT FINITE TEMPERATURES A. Temperature dependence
Using MD simulations one can gain insight into the elastic properties of materials under different kinds of applied stresses, e.g., hydrostatic or uniaxial.In particular, we consider elastic stiffness constants, c ij , and compliance constants, s ij , for 2D crystalline materials with hexagonal structure such as graphene.We call τ ij and e ij the components of the stress and strain tensors, respectively.τ ij is the force per unit length parallel to direction i, acting in the xy plane on a line perpendicular to the j direction.We use the standard notation for strain components, with e ij = ǫ ij for i = j, and e ij = 2ǫ ij for i = j [79,80].More details on elastic properties of 2D crystals can be found in Ref. [54].
In terms of the compliance constants, we have for applied stress {τ ij }: ( The matrix of stiffness constants {c ij } is the inverse of {s ij }, so that we have the relations In Fig. 8 we present the stiffness constants as a function of temperature, as derived from our MD simulations of bilayer graphene, using Eq. ( 12) (open circles).Panels (a) and (b) show results for c 11 and c 12 , respectively.Solid squares at T = 0, signaled by arrows, indicate results for c 11 and c 12 obtained from the phonon dispersion bands as indicated in Sec.II, using Eqs.( 1) and (2).We find that finite-temperature data for the stiffness constants converge at low T to the results of the HA for both c 11 and c 12 .For rising temperature, the stiffness constants decrease rather fast.This decrease is especially large for c 12 , which is found to be 1.18 eV/ Å2 at T = 1200 K vs the classical low-T limit of 4.54 eV/ Å2 .
Comparing the elastic constants c 11 and c 12 found here for bilayer graphene with those corresponding to monolayer graphene [81] and graphite [67] (normalized to one layer), we find that they increase for the sequence monolayer-bilayer-graphite.This agrees with the fact that interaction between layers reduces the amplitude of out-of-plane vibrational modes, thus favoring an increase in the "hardness" of the layers.This trend is similar to that discussed below for the 2D compression modulus B p .The Poisson's ratio ν can be obtained as the quotient c 12 /c 11 .This yields for T = 0 (HA) ν = 0.22.From the results of our simulations, we find ν = 0.15 and 0.09 for T = 300 and 1000 K, respectively, with an important reduction for rising T , as a consequence of the decrease in c 12 .Calculations based on the self-consistent screening approximation [82][83][84] (SCSA) predict for P = 0 in the thermodynamic limit (N → ∞) a Poisson's ratio ν = −1/3.A negative value for this ratio is also expected from the calculations presented by Burmistrov et al. [85] for N → ∞.From the results of our MD simulations, we do not find, however, any indication for a negative Poisson's ratio in the parameter region considered here.This is in line with earlier results of Monte Carlo simulations by Los et al. [26] for monolayer graphene in a region of system sizes larger than those considered here for bilayer graphene.
The 2D modulus of hydrostatic compression B p is defined for layered materials at temperature T as [54] Note the factor n (number of sheets) in the denominator, i.e., B p is the compression modulus per layer.P and S p appearing on the r.h.s. of Eq. ( 15) are variables associated to the layer plane, and in fact the pressure in the isothermal-isobaric ensemble used here is the conjugate variable to the area S p .One can also calculate the modulus B p on the basis of the fluctuation formula [48,86,87]  16).Open squares are data points obtained from the elastic constants s11 and s12.The dashed line is a guide to the eye.
where (∆S p ) 2 is the mean-square fluctuation of the area S p , which is calculated here from MD simulations at P = 0.This formula provides us with a practical procedure to obtain B p , vs calculating the derivative (∂S p /∂P ) T by numerical methods, which requires additional MD simulations at hydrostatic pressures close to P = 0.For some temperatures we have verified that results for B p found with both procedures coincide within statistical error bars, which is a consistency check for our results.
The modulus B p can be also obtained from the elastic constants of bilayer graphene.Taking into account Eq. ( 12), the change of the in-plane area, ∆S p due to a 2D hydrostatic pressure P , is given by ∆S p S p = e xx + e yy = −2(s 11 + s 12 )P .
Combining Eqs. ( 15) and ( 17), one finds which can be also written as B p = (c 11 +c 12 )/2.These expressions are valid for 2D materials with hexagonal symmetry.
In Fig. 9 we present the modulus B p of bilayer graphene as a function of T , as derived from our MD simulations.Solid circles represent results obtained from the in-plane area fluctuation (∆S p ) 2 by employing Eq. ( 16).Open squares are data points calculated from the elastic constants.Both sets of results coincide within error bars.At low temperature, B p converges to the value given by the expression where E is the energy.For bilayer graphene we have B 0 = 12.74 eV/ Å2 , which agrees with the extrapolation of finite-T results to T = 0.The modulus B p derived from MD simulations is found to decrease fast as the temperature is raised, and at T = 1200 K it amounts to about 60% of the low-T limit B 0 .monolayer graphene, the same interatomic potential yields in a HA: B 0 = 12.65 eV/ Å2 , somewhat less than the value found for the bilayer.This difference is larger at finite temperatures.Even though interlayer interactions are relatively weak, they give rise to a reduction in the vibrational amplitudes of out-of-plane modes, and as a result the graphene sheets become "harder" in the bilayer, so that the modulus B p increases with respect to an isolated sheet.Moreover, the difference between the modulus B p per sheet for bilayer and monolayer graphene grows for rising system size N (see Sec. VII).
The in-plane Young's modulus E p can be obtained from B p through the expression E p = 2B p (1 − ν).This yields for the bilayer at T = 0, E p = 19.87eV/ Å2 , which translated into units of force per square length gives E p /d = 0.954 TPa, similar to values appearing in the literature [37,56].

B. Mechanical instability under stress
The modulus B p is particularly interesting to study the critical behavior of bilayer graphene under 2D hydrostatic pressure.In Fig. 10 we present the dependence of B p on P , including tensile and compressive stresses.Symbols represent values derived from MD simulations for various cell sizes.From left to right: N = 8400, 960, 448, 308, and 240.For each size N , increasing in-plane compression causes a fast decrease in B p , which vanishes for a pressure P c (N ), where the bilayer graphene with planar sheets becomes mechanically unstable.This is typical of a spinodal point in the (P, T ) phase diagram [18,[88][89][90].For P > P c , the stable configuration corresponds to wrinkled graphene sheets, as observed earlier for monolayer graphene [48].
For given N and T , there is a pressure region (compressive stress) where bilayer graphene is metastable, i.e., for P < P c .The spinodal line, which delineates the metastable phase from the unstable phase, is the locus of points P c (N, T ) where B p = 0.This kind of spinodal lines have been studied earlier for water [91], as well as for ice, SiO 2 cristobalite [88], and noble-gas solids [89] near their stability limits.In recent years, this question has been investigated for 2D materials, and in particular for monolayer graphene [18,92].
According to Eq. ( 16), vanishing of B p for finite N corresponds to a divergence of the area fluctuation (∆S p ) 2 to infinity.Moreover, the MSF of the atomic z coordinate, (∆h) 2 , diverges at the corresponding spinodal pressure, as mentioned in Sec.IV.C.For graphene bilayers, this spinodal instability is associated to a soft vibrational mode in the ZA branch.In fact, for each N this instability appears for increasing P when the frequency of the ZA vibrational mode with minimum wavenumber, k min , reaches zero (ω min → 0).
Close to a spinodal point, the modulus B p behaves as a function of P as B p ∼ (P c − P ) 1/2 (see Appendix C).This pressure dependence agrees with the shape of the curves shown in Fig. 10 near the spinodal pressure P c for each size N .Note that P c moves to smaller compressive pressures as N rises.This size effect is analyzed below in Sec.VII.We also note that, for a given size N , the critical stress P c depends on temperature, as shown before for monolayer graphene [93].It was found that P c increases for rising temperature, as a consequence of a raise in vibrational amplitudes in the out-of-plane direction.We have checked that something similar happens for bilayer graphene, but a detailed study of this question requires additional MD simulations, which will be the subject of future work.

VII. SIZE EFFECTS
As noted above, some properties of 2D materials display important size effects.In this section, we concentrate on the size dependence of the MSF (∆h) 2 , the modulus B p , and the spinodal pressure P c for bilayer graphene, and study their asymptotic behavior for large N .
In Fig. 11 we present in a logarithmic plot the atomic MSF (∆h) 2 in the z direction as a function of system For comparison we also display data for monolayer graphene (circles) and graphite (diamonds), obtained from MD simulations with the LCBOPII interatomic potential.For the system sizes presented in Fig. 11, (∆h) 2 may be expressed in the three cases as a power of N : (∆h) 2 ∼ N α .We find for the exponent α values of 0.78, 0.69, and 0.56 for monolayer, bilayer graphene, and graphite, respectively.The MSF in the z direction can be written in a HA as: where the sum in j is extended to the phonon bands with atomic motion in the out-of-plane direction, i.e., ZA, ZO, and ZO' for bilayer graphene (as above in Eq. ( 11)).
The sum in Eq. ( 20) is dominated by ZA modes with wavevector close to the Γ point, i.e. small frequency ω.
The inputs of bands ZO and ZO' are almost independent of the system size, and they give a joint contribution of ≈ 8 × 10 −3 Å2 to (∆h) 2 in Eq. (20).
For the ZA band, putting a dispersion ρ ω 2 = κk 4 , one finds [72] (∆h where C = 6.03 is a constant.Thus, in a HA the ZA band yields a contribution proportional to N and an exponent α = 1.The result of the HA including the inputs of the three phonon bands ZA, ZO, and ZO' is shown in Fig. 11 as a dotted line.For large N , (∆h) 2 is dominated by atomic displacements associated to ZA modes and we find that it increases linearly with N .For small N , one observes a departure from the linear trend due to contributions of ZO and ZO' modes.
The size dependence of (∆h) 2 obtained from MD simulations for stress-free bilayer graphene can be understood assuming an effective dependence for the frequency of ZA modes as ρ ω 2 = κ k β , where κ is a modified bending constant and β is an exponent controlling the frequency of long-wavelength (small frequency) modes.This expression assumes an effective shape for the ZA band at finite temperatures, as a consequence of anharmonicity in the vibrational modes, and is similar to that considered earlier for monolayer graphene [26,72,94].Assuming such a dispersion for ZA modes in the bilayer, we have for large N : where the sum is extended to k points in the 2D Brillouin zone.
Taking into account the relation between the minimum wavenumber k min and the size N , and replacing the sum in Eq. ( 22) by an integral, one finds a size dependence where D is an integration constant.This means that our exponent α can be related to β as α = β/2 − 1, which yields for bilayer graphene β = 3.38.Similar effective exponents can be derived from MD simulation results for monolayer graphene and graphite, for which we find β = 3.56 and 3.12, respectively.
The 2D modulus of hydrostatic compression B p introduced in Sec.VI also displays finite-size effects.In Fig. 12 we show in a logarithmic plot the dependence of B p on N at T = 300 K. Open squares represent results obtained for bilayer graphene from MD simulations.For comparison, we also display data for monolayer graphene (circles), as well as for graphite (diamonds).For N > 500, B p can be fitted for the bilayer to an expression B p ∼ N −ζ , with an exponent ζ = 0.086.From similar fits for monolayer graphene and graphite, we find the exponents 0.159 and 0.033, respectively.Note that the exponent ζ for the bilayer is about one half of that corresponding to the monolayer.This indicates that the size effect is less important for the former than for the latter, as visualized in Fig. 12.
Looking at Eq. ( 16), and taking into account that S p changes slowly with N , we have for the area fluctuation a size dependence: , where the subscripts B and M refer to bilayer and monolayer, respectively.For small N , the area fluctuation is similar for bilayer and monolayer graphene, but they become comparatively smaller for the bilayer as the size increases.Our exponent for the monolayer can be translated to a dependence on the cell side length B p ∼ L −2ζ , with 2ζ = 0.318, close to the exponent 0.323 found by Los et al. [26] for the dependence of B p on L.
The size dependence of the critical pressure P c introduced in Sec.VI is shown in Fig. 13, where we have plotted values of P c derived from MD simulations at T = 300 K for several cell sizes.One observes at first sight a linear dependence of P c with the inverse cell size, N −1 .This dependence may be understood by considering the effect of a compressive stress on ZA vibrational modes, as follows.
For a single graphene layer under a 2D hydrostatic pressure P , the dispersion relation of ZA modes may be written as where σ = −P [18,81].For increasing compressive stress, ω is reduced, i.e. dω/dP < 0. Thus, for system size N , a graphene layer becomes unstable when the frequency of the ZA mode with wavenumber k min vanishes.This occurs for which yields a critical stress As indicated above, the minimum wavenumber k present for cell size N is k min = 2π/(N S p ) 1/2 .For bilayer graphene, the in-plane critical pressure is given by P c = −2σ c , from where we have: The solid line was calculated using Eq. ( 27).
Putting κ = 1.49eV and S p = 2.617 Å2 /atom, we find P c N = 44.95eV/ Å2 , which is the line displayed in Fig. 13.This line matches well the results of our MD simulations (solid circles), with the exception of the result for the largest size presented in the figure.This deviation from the general trend of smaller sizes may be due to three reasons.First, the presence in the graphene layers of a residual (small) intrinsic stress at T = 300 K, which is not detected for small N , due to the larger values of the corresponding pressure P c [46,93].Second, the graphene bilayer can remain in a metastable state along millions of MD simulation steps for large cell size.This means that observation of the true transition (spinodal) point could require much longer simulations, not available at present for such large system sizes.Third, the dispersion relation for the ZA band in Eq. ( 24), utilized to obtain Eq. ( 27), may be modified for small k (long wavelength), so that the exponent 4 on the r.h.s.could be renormalized in a similar way to the exponent β found for the size dependence of (∆h) 2 .Calculations based on the SCSA predict a universal behavior for scaling exponents in the thermodynamic limit (N → ∞) [82][83][84].This means that the exponents presented above should coincide for 2D systems (including monolayers and bilayers) in the large-size limit.According to such calculations, universality is approached for system size larger than a crossover scale given by the so-called Ginzburg length, L G .This temperaturedependent length can be estimated for graphene (using the bending constant κ and Young's modulus E p ) to be around 40−50 Å at T = 300 K [95,96].This corresponds in our notation to a system size N G ∼ 900.For bilayer graphene, we have considered here simulation cells with length sides up to 150 Å, well above those values of L G .One can, however, understand L G as a reference length for the crossover to a regime where universality is ap-proached, and a direct detection of this universality could be only found for lengths clearly larger than L G .In any case, from the results of our simulations for bilayer (with L up to 150 Å) and monolayer graphene presented here, we do not find any evidence or trend indicating that such a kind of universality will appear for results derived from simulations using larger cells.Thus, if such kind of universality is in fact a physical aspect of 2D crystalline membranes, as predicted by the SCSA, it has not yet been observed from atomistic simulations with interaction potentials mimicking those of actual materials, as graphene.

VIII. SUMMARY
MD simulations allow us to gain insight into elastic properties of 2D materials, as well as on their stability under external stress.We have presented here the results of extensive simulations of bilayer graphene using a well-checked interatomic potential, for a wide range of temperatures, in-plane stresses, and system sizes.We have concentrated on physical properties such as the excess area, interlayer spacing, interatomic distance, elastic constants, in-plane compression modulus, and atomic MSF in the out-of-plane direction.
The elastic constants are found to appreciably change as a function of temperature, especially c 12 .This causes a reduction of the Poisson's ratio for rising T .The inplane compression modulus B p has been obtained from the fluctuations of the in-plane area, a procedure which yields results consistent with those derived from the elastic constants of the bilayer.
For bilayer graphene under in-plane stress, we find a divergence of the MSF (∆h) 2 for an in-plane pressure P c (N ), which corresponds to the limit of mechanical stability of the material.This divergence is accompanied by a vanishing of the in-plane compression modulus, or a divergence of the compressibility ξ p = 1/B p .
Finite-size effects are found to be important for several properties of bilayer graphene.The spinodal pressure P c is found to scale with system size as 1/N .A similar scaling with the inverse size is obtained for the MSF of the interlayer spacing: (∆d) 2 ∼ N −1 .The atomic outof-plane MSF also follows a power law (∆h) 2 ∼ N α with an exponent α = 0.69.For B p , we find for N > 500 at T = 300 K a dependence B p ∼ N −ζ , with an exponent ζ = 0.086.
Comparing the simulation results with those obtained from a HA gives insight into finite-temperature anharmonic effects.Thus, for the atomic MSF in the out-ofplane direction, a HA predicts a linear dependence of (∆h) 2 with system size N , to be compared with the sublinear dependence obtained from the simulations.The change with system size N of (∆h) 2 and the modulus B p for bilayer graphene is slower (i.e., less slope in Figs.11 and 12) than for the monolayer.This is indeed due to interlayer interactions, which manifests themselves in the presence of the layer-breathing ZO' phonon branch.According to calculations based on the SCSA, the size dependence of physical observables such as (∆h) 2 or the in-plane modulus B p should be controlled, for large N , by universal exponents independent of the particular details of the considered 2D system (monolayer or bilayer graphene in our case).We have not observed this universality from our MD simulations for cell size up to 150 Å, and a clarification of this question remains as a challenge for future research.
We finally note that MD simulations as those presented here can give information on the properties of graphene multilayers under stress.This may yield insight into the relative stability of such multilayers in a pressure-temperature phase diagram.Moreover, nuclear quantum effects can affect the mechanical properties of graphene bilayers and multilayers at low temperatures, as shown earlier for graphite.This question can be addressed using atomistic simulations with techniques such as path-integral molecular dynamics.

FIG. 1 :
FIG.1: Top (upper) and side (lower) views of an atomic configuration of bilayer graphene obtained from MD simulations at T = 800 K.In the top view, red and yellow spheres represent carbon atoms in the upper and lower graphene sheets, respectively.

FIG. 2 :
FIG.2: Phonon dispersion bands of bilayer graphene, obtained from diagonalization of the dynamical matrix for the LCBOPII potential model.

FIG. 4 :
FIG.4: MSF of the interlayer distance as a function of the inverse cell size for T = 100 K (circles) and 300 K (squares), as derived from MD simulations of bilayer graphene.

FIG. 6 :
FIG.6: In-plane area per atom vs 2D hydrostatic pressure P at T = 300 K. Shown are results for cell sizes N = 308 (traingles), 448 (squares), 960 (diamonds), and 8400 (circles).Error bars, when not displayed, are in the order of the symbol size.Dashed lines are guides to the eye.

FIG. 7 :
FIG.7: Excess area per atom vs 2D hydrostatic pressure P for cell size N = 960 at T = 300 K (circles), 600 K (squares), and 1000 K (diamonds).Open symbols are data points derived from MD simulations of bilayer graphene.Dashed lines are guides to the eye.

FIG. 8 :
FIG. 8: Temperature dependence of the elastic stiffness constants of bilayer graphene, as derived from MD simulations for N = 960 (open circles): (a) c11, (b) c12.Dashed lines are guides to the eye.Solid squares at T = 0, indicated by arrows, show the values of c11 and c12 derived from the HA using Eqs.(1) and (2).

FIG. 9 :
FIG. 9: 2D modulus of hydrostatic compression, Bp, as a function of temperature for stress-free bilayer graphene (P = 0) and N = 960.Solid circles represent results derived from fluctuations of the in-plane area, using Eq.(16).Open squares are data points obtained from the elastic constants s11 and s12.The dashed line is a guide to the eye.

FIG. 10 :
FIG. 10: Modulus Bp as a function of 2D hydrostatic pressure at T = 300 K, for various cell sizes.Symbols are data points derived from MD simulations for various cell sizes.From left to right: N = 8400, 960, 448, 308, and 240.Dashed lines are guides to the eye.

FIG. 11 :
FIG.11: Atomic MSF in the z direction vs cell size for T = 300 K in a logarithmic plot.Data derived from MD simulations are given for monolayer graphene (circles), bilayer graphene (squares), and graphite (diamonds).Error bars are in the order of the symbol size.Dashed lines are least-square fits to the data points.The dotted line represents the MSF corresponding to a HA.

FIG. 12 :
FIG.12: Modulus Bp as a function of cell size N at T = 300 K on a logarithmic plot.Symbols represent data derived from MD simulations for a graphene monolayer (circles), bilayer (squares), and graphite (diamonds).Dashed lines are guides to the eye.

2 FIG. 13 :
FIG.13: Critical pressure Pc as a function of the inverse cell size at T = 300 K. Solid circles represent data derived from MD simulations of bilayer graphene under in-plane stress.The solid line was calculated using Eq.(27).
using the interlayer distance d 0 of