Nuclear quantum effects in graphene bilayers

TL;DR

This study uses path-integral molecular dynamics to analyze nuclear quantum effects in graphene bilayers, revealing significant quantum contributions to structural properties and vibrational modes across a wide temperature range.

Contribution

It provides the first detailed quantification of nuclear quantum effects in graphene bilayers, including zero-point expansion and anharmonic vibrational behavior, using advanced simulation techniques.

Findings

Nuclear quantum effects cause a 0.015 Å zero-point expansion in interlayer spacing.

Out-of-plane compressibility increases with temperature, with a 6% rise at low temperatures due to quantum effects.

Quantum effects are significant in vibrational modes but diminish with system size.

Abstract

Graphene bilayers display peculiar electronic and mechanical characteristics associated to their two-dimensional character and relative disposition of the sheets. Here we study nuclear quantum effects in graphene bilayers by using path-integral molecular dynamics simulations, which allow us to consider quantization of vibrational modes and study the effect of anharmonicity on physical variables. Finite-temperature properties are analyzed in the range from 12 to 2000~K. Our results for graphene bilayers are compared with those found for graphene monolayers and graphite. Nuclear quantum effects turn out to be appreciable in the layer area and interlayer distance at finite temperatures. Differences in the behavior of in-plane and real areas of the graphene sheets are discussed. The interlayer spacing has a zero-point expansion of $1.5 \times 10^{-2}$ \AA\ with respect to the classical minimum. The compressibility of graphene bilayers in the out-of-plane direction is found to be similar to that of graphite at low temperature, and increases faster as temperature is raised. The low-temperature compressibility increases by a 6\% due to zero-point motion. Especial emphasis is laid upon atomic vibrations in the out-of-plane direction. Quantum effects are present in these vibrational modes, but classical thermal motion becomes dominant over quantum delocalization for large system size. The significance of anharmonicities in this atomic motion is estimated by comparing with a harmonic approximation for the vibrational modes in graphene bilayers.

Tables

Table 1. Table 1: Contributions to the internal energy of bilayer graphene at different temperatures, as derived from PIMD simulations for N = 960 𝑁 960 N=960 . E 0 subscript 𝐸 0 E_{0} is the energy minimum for an isolated graphene sheet, taken as reference for the internal energy. Units of energy are meV/atom. Error bars are in all cases less than ± 0.05 plus-or-minus 0.05 \pm 0.05  meV/atom.

$T$ (K)	$E_{\mathrm{int}}$	$E_{\mathrm{el}}$	$E_{\mathrm{vib}}$	$E-E_0$
50	-25.0	1.5	170.8	147.3
100	-25.0	1.5	171.5	148.0
300	-25.0	1.6	182.8	159.4
500	-24.9	1.9	206.9	183.9
750	-24.7	2.7	249.2	227.2
1000	-24.4	3.8	299.6	279.0
1500	-23.9	7.6	413.0	396.7
2000	-23.1	13.5	534.7	525.1

Equations

$$A=⟨\frac{1}{N_{\mathrm{Tr}}}\sum _{j=1}^{N_{\mathrm{Tr}}}{A}^j⟩,$$

$${\overline{\mathbf{r}}}_i=\frac{1}{N_{\mathrm{Tr}}}\sum _{j=1}^{N_{\mathrm{Tr}}}{\mathbf{r}}_{ij},$$

$$(\Delta z{)}_i^2=\frac{1}{N_{\mathrm{Tr}}}\sum _{j=1}^{N_{\mathrm{Tr}}}⟨(z_{ij}-⟨{\overline{z}}_i⟩{)}^2⟩,$$

$$Q_{z,i}^2=\frac{1}{N_{\mathrm{Tr}}}\sum _{j=1}^{N_{\mathrm{Tr}}}⟨(z_{ij}-{\overline{z}}_i{)}^2⟩.$$

$$(\Delta z{)}_i^2=C_{z,i}^2+Q_{z,i}^2,$$

$$C_{z,i}^2=⟨{\left({\overline{z}}_i-⟨{\overline{z}}_i⟩\right)}^2⟩=⟨{\overline{z}}_i^2⟩-⟨{\overline{z}}_i{⟩}^2.$$

$$(\Delta z{)}^2=\frac{1}{2N}\sum _{i=1}^{2N}(\Delta z{)}_i^2,$$

$$(\Delta z{)}^2=\frac{1}{2N}\sum _{i,\mathbf{k}}\frac{\hslash }{2m{\omega }_i(\mathbf{k})}\coth \left(\frac{\hslash {\omega }_i(\mathbf{k})}{2k_{B}T}\right),$$

$$C_z^2=\frac{1}{2N}\sum _{i,\mathbf{k}}\frac{k_{B}T}{m{\omega }_i(\mathbf{k}{)}^2}.$$

$$(\Delta z{)}_0^2=Q_{z,0}^2=\frac{1}{2N}\sum _{i,\mathbf{k}}\frac{\hslash }{2m{\omega }_i(\mathbf{k})}.$$

$$E=E_0+E_{\mathrm{int}}+E_{\mathrm{el}}+E_{\mathrm{vib}},$$

$${\chi }_z=-\frac{1}{⟨V⟩}\frac{\partial ⟨V⟩}{\partial P_z},$$

$${\chi }_z=\frac{(\Delta V{)}^2}{k_{B}T⟨V⟩},$$

$${\chi }_z=\frac{L_x{L}_y}{k_{B}T}\frac{(\Delta c{)}^2}{⟨c⟩}.$$

$$E_{\mathrm{int}}=E_{\mathrm{int}}^0+\frac{1}{2}k(c-c_0{)}^2$$

$${\chi }_0=-\frac{1}{V_0}{\left(\frac{\partial V}{\partial P_z}\right)}_0=\frac{1}{N{A}_0{c}_0}{\left(\frac{{\partial }^2{\bar{E}}_{\mathrm{int}}}{\partial V^2}\right)}_0^{-1},$$

$${\chi }_0=\frac{A_0}{2c_{0}k},$$

$${\chi }_z=\frac{(\Delta V{)}^2}{k_{B}T⟨V⟩},$$

$$(\Delta c{)}^2=\frac{k_{B}T}{2Nk}$$

Full text

I Introduction

In recent years there has been amid scientists a surge of interest in two-dimensional (2D) materials, which display thicknesses in the atomic scale, reaching in some cases the level of one atomic layer. Among these materials, graphene is one of the most studied due to its remarkable electronic,Geim and Novoselov (2007); Castro Neto et al. (2009) mechanical,Lee et al. (2008); Gao et al. (2018) and thermal properties.Ghosh et al. (2008); Seol et al. (2010); Balandin (2011) An additional advantage of graphene, as compared with other layered systems, is that monolayers, bilayers and multilayers can be produced in a controlled manner,Zhang et al. (2016) which allows one to carry out detailed studies of their physical properties.

Bilayer graphene, in particular, presents peculiar electronic properties different from those of monolayer graphene and graphite.Novoselov et al. (2006) It is known that it displays ripples and out-of-plane deformations like the monolayers,Meyer et al. (2007) which are thought to be a relevant scattering mechanism for electrons in this material.Gibertini et al. (2010) Recently, graphene bilayers have attracted increased interest after the discovery that they display unconventional superconductivity by stacking both sheets twisted relative to each other by a small angle.Cao et al. (2018) Moreover, it has been lately found the appearance of Mott-like insulator states, due to the presence of electrons localized in the superlattice associated to a moiré pattern,Cao et al. (2018) as well as uncommon magnetic properties.Yndurain (2019)

A deep comprehension of structural and thermodynamic properties of 2D systems in three-dimensional space has been for many years a persistent goal in statistical physics.Safran (1994); Nelson et al. (2004) This question has been largely discussed in the context of biological membranes and soft condensed matter.Chacón et al. (2015); Ruiz-Herrero et al. (2012) Nevertheless, the complexity of these systems makes it difficult to work out microscopic approaches based on realistic interatomic interactions. Graphene bilayers constitute a singular realization of crystalline membranes composed of two atomic layers, where an atomistic description is feasible, thus paving the way for a deeper insight into the physical properties of this type of systems.Fasolino et al. (2007); Balandin (2011); Alofi and Srivastava (2013); Amorim et al. (2014); Mann et al. (2016) Additionally, graphene can give information on the thermodynamic stability of 2D crystals in general, which has been debated and may be related to the anharmonic coupling between in-plane and out-of-plane vibrational modes.Amorim et al. (2014); Coquand and Mouhanna (2016)

Atomistic simulations have been used to understand finite-temperature properties of graphene.Fasolino et al. (2007); Cadelano et al. (2009); Akatyeva and Dumitrica (2012); Chechin et al. (2014); Magnin et al. (2014); Los et al. (2016) In general, carbon atoms were treated in those simulations as classical particles. However, the Debye temperature of graphene for out-of-plane vibrational modes is $\Theta_{D}^{\rm out}\gtrsim$ 1000 K,Politano et al. (2011) which suggests that quantum fluctuations should be relevant up to relatively high temperatures. The quantum character of atomic motion can be taken into account by using path-integral simulations, which allow one to consider thermal and quantum fluctuations at finite temperatures.Gillan (1988); Ceperley (1995) In these methods the nuclear degrees of freedom can be quantized in an efficient way, thus permitting to carry out quantitative analyses of anharmonic effects in condensed matter. This kind of simulations were carried out for graphene monolayers in recent years to study structural and thermodynamic properties of this material.Brito et al. (2015); Herrero and Ramírez (2016, 2018a); Hasik et al. (2018)

Here we extend those studies to graphene bilayers, for which new features appear as a consequence of the interaction between layers and the resultant coupling between atomic displacements in the out-of-plane direction. We use the path-integral molecular dynamics (PIMD) method to investigate the properties of graphene bilayers at temperatures between 12 and 2000 K. Simulation cells of different sizes are considered, as finite-size effects are known to be important for some properties of graphene.Gao and Huang (2014); Herrero and Ramírez (2016); Los et al. (2016) The thermal behavior of the graphene surface is studied, considering the difference between in-plane and real areas. Moreover, the interlayer distance and its thermal fluctuations yield information about the compressibility of the bilayer in the out-of-plane direction. Our results for the bilayer are compared with those obtained for graphene monolayers and graphite, which provides information on the transition of physical properties from an isolated graphene sheet to the bulk material. Results of the simulations are also compared with predictions based on harmonic vibrations of the crystalline sheets. This approximation happens to be noticeably accurate at low temperatures, once the frequencies of out-of-plane modes in graphene bilayers are properly described.

Path-integral methods similar to that used here have been previously applied to investigate nuclear quantum effects in pure and doped carbon-based materials, such as diamondRamírez et al. (2006) and graphite.Herrero and Ramírez (2010) More recently, these computational techniques have been also employed to study adsorption of helium and hydrogen on graphene.Kwon and Ceperley (2012); Herrero and Ramírez (2009); Davidson et al. (2014)

The paper is organized as follows. In Sec. II, we describe the computational methods and details of the calculations. In Sec. III we present our results for the in-plane and real area of bilayer graphene. In Sec. IV we present data of the internal energy and the different contributions to it. The interlayer spacing and compressibility of the bilayer in the out-of-plane direction is discussed in Sec. V, and the character of the out-of-plane atomic motion (classical vs quantum) is dealt with in Sec. VI. Finally we summarize the main results in Sec. VII.

II Computational Method

II.1 Path-integral molecular dynamics

We use PIMD simulations to study equilibrium properties of graphene bilayers at various temperatures. This procedure is based on the path-integral formulation of statistical mechanics, which is a suitable nonperturbative approach to study finite-temperature properties of many-body quantum systems.Feynman (1972) In actual applications of this method to numerical simulations, each quantum particle is described as a set of $N_{\rm Tr}$ beads (the so-called Trotter number), which behave as classical-like particles arranged to build a ring polymer.Gillan (1988); Ceperley (1995) This representation becomes formally exact in the limit $N_{\rm Tr}\to\infty$. Details on this kind of simulation techniques can be found in Refs. Gillan, 1988; Ceperley, 1995; Herrero and Ramírez, 2014; Cazorla and Boronat, 2017.

Here we use molecular dynamics simulations to sample the configuration space of the classical isomorph of our quantum system ($2N$ carbon atoms). The dynamics in this computational technique is artificial, as it does not represent the real quantum dynamics of the actual particles under consideration. However, it turns out to be very effective to sample the many-body configuration space, thus yielding accurate results for time-independent equilibrium properties of the quantum system.

In the context of PIMD simulations, one needs a Born-Oppenheimer surface for the nuclear dynamics, which should be derived from an adequate description of the interatomic interactions. Employing an ab-initio method would largely restrict the size of the manageable simulation cells, so we use the LCBOPII effective potential, a long-range bond order potential, mainly employed to carry out classical simulations of carbon-based systems.Los et al. (2005) It has been used, in particular, to study the phase diagram of carbon (diamond, graphite, and liquid carbon), displaying its reliability by predicting rather accurately the diamond-graphite transition line.Ghiringhelli et al. (2005) In the last few years, the LCBOPII potential was also employed to reliably describe various properties of graphene,Fasolino et al. (2007); Los et al. (2016) in particular its Young’s modulus.Ramírez and Herrero (2017); Zakharchenko et al. (2009); Politano et al. (2012); Ramírez and Herrero (2018) This potential has been lately used to perform PIMD simulations, aiming at assessing the magnitude of nuclear quantum effects in graphene monolayers,Herrero and Ramírez (2016) and to study thermodynamic properties of this material.Herrero and Ramírez (2018a) In this paper, according to previous simulations,Ramírez et al. (2016); Herrero and Ramírez (2016); Ramírez and Herrero (2017) the earliest parameterization of the LCBOPII potential has been slightly modified to increase the zero-temperature bending constant $\kappa$ of a graphene monolayer from 0.82 to 1.49 eV, a value closer to experimental data and ab-initio calculations.Lambin (2014)

Calculations have been carried out in the isothermal-isobaric ensemble, where we fix the number of C atoms ($2N$), the in-plane applied stress (here $P_{xy}=0$), and the temperature ($T$). We have employed effective algorithms for performing PIMD simulations in this ensemble, such as those described in the literature for this kind of simulations.Tuckerman and Hughes (1998); Martyna et al. (1999); Tuckerman (2002) Staging variables were used to define the bead coordinates, and the constant-temperature ensemble was accomplished by coupling chains of four Nosé-Hoover thermostats to each staging variable. An additional chain of four barostats was coupled to the in-plane area of the simulation box ($xy$ plane) to give the constant pressure $P_{xy}=0$.Tuckerman and Hughes (1998); Herrero and Ramírez (2014) The kinetic energy $E_{k}$ has been calculated by using the virial estimator, which displays a statistical uncertainty much smaller than the primitive estimator, especially at high temperatures.Herman et al. (1982); Tuckerman and Hughes (1998) This means that the error bar associated to the kinetic energy in our calculations at $T>100$ K is smaller than that corresponding to the potential energy of the system. The relative accuracy of the kinetic energy, as compared with the potential energy, increases for rising temperature. Both error bars are similar at $T=$ 50 K, whereas at 1000 K and 2000 K the error bar of the kinetic energy is 8 and 25 times less than that of the potential energy, respectively. Other technical details about this type of simulations can be found elsewhere.Herrero et al. (2006); Herrero and Ramírez (2011); Ramírez et al. (2012)

We have carried out PIMD simulations of graphene bilayers with AB stacking employing simulations boxes with $2N$ carbon atoms, $N$ going from 24 to 8400. The considered cells had similar side lengths in the $x$ and $y$ directions ($L_{x}\approx L_{y}$), for which periodic boundary conditions were assumed. Carbon atoms can unrestrictedly move in the out-of-plane direction, thus simulating a free-standing graphene bilayer, i.e., we have free boundary conditions in the $z$ coordinate. The temperature $T$ was in the range from 12.5 to 2000 K. For a given temperature, a typical simulation run included $2\times 10^{5}$ PIMD steps for system equilibration, followed by $8\times 10^{6}$ steps for the calculation of average properties. The Trotter number $N_{\rm Tr}$ was taken proportional to the inverse temperature, as $N_{\rm Tr}T$ = 6000 K, which roughly keeps a constant precision for the PIMD results at different temperatures. The time step $\Delta t$ for the molecular dynamics has been taken as 0.5 fs, which is adequate for the C atomic mass and the temperatures considered here. To compare with the results of these simulations, some classical molecular dynamics simulations of graphene bilayers were also performed. In our context this corresponds to setting $N_{\rm Tr}$ = 1.

To compare with our results for graphene bilayers, we have also carried out some PIMD simulations of graphite with the same interatomic potential LCBOPII. In this case we employed cells consisting of $4N$ carbon atoms, i.e., four graphene sheets, and periodic boundary conditions were applied in the three space directions. We used cells with $N$ = 240 and 960.

II.2 Layer area

In the isothermal-isobaric ensemble employed here one fixes the applied stress $P_{xy}$ in the $xy$ plane, as indicated above, which allows for fluctuations in the in-plane area of the simulation cell, given by $L_{x}L_{y}$. Carbon atoms can unrestrictedly move in the $z$ direction (out-of-plane direction), so a measure of the real surface of a graphene sheet at finite temperatures will give a value larger than the in-plane area.

The concept of a real surface has been discussed for biological membranes as an interesting tool to describe some of their physical properties, instead of the projected in-plane surface.Waheed and Edholm (2009); Chacón et al. (2015) Something similar has been proposed in recent years for crystalline membranes such as graphene.Pozzo et al. (2011); Herrero and Ramírez (2016); Nicholl et al. (2017); Ramírez and Herrero (2017) For biological membranes, it was shown that values of the compressibility may appreciably differ when they are related to the real area $A$ or in-plane area $A_{p}$, and something analogous has been noticed recently for the elastic properties of graphene monolayers, as derived from classical molecular dynamics simulations.Ramírez and Herrero (2017) The real area has been also called true, actual, or effective area in the literature.Fournier and Barbetta (2008); Waheed and Edholm (2009); Chacón et al. (2015) A clear distinction between both areas is basic to understand some thermodynamic properties of 2D materials. In fact, the area $A_{p}$ is the variable conjugate to the in-plane stress $P_{xy}$ in the isothermal-isobaric ensemble employed here, whereas the real area $A$ is conjugate to the usually-called surface tension.Safran (1994) Nicholl et al.Nicholl et al. (2015, 2017) have shown that certain experimental techniques are sensitive to properties related to the real area $A$, whereas other methods may be suitable to quantify variables associated to the in-plane area $A_{p}$. The difference between both areas, $A-A_{p}$, has been named hidden area for graphene in Ref. Nicholl et al., 2017, and excess area in the context of biological membranes.Helfrich and Servuss (1984); Fournier and Barbetta (2008)

For each graphene sheet, we calculate the real area $A$ by a triangulation based on the positions of the C atoms along a simulation run. In this procedure, $A$ is obtained as a sum of areas of the hexagons in the graphene structure.Ramírez and Herrero (2017); Herrero and Ramírez (2018b) Each hexagon contributes as a sum of six triangles, each one formed by the positions of two adjacent C atoms and the barycenter of the hexagon.Ramírez and Herrero (2017); Herrero and Ramírez (2018b) In our path-integral method, the area $A$ is given as an average over the $N_{\rm Tr}$ beads associated to the atomic nuclei:

[TABLE]

where $A^{j}$ is the instantaneous area per atom for imaginary time (bead) $j$, and the brackets indicate an ensemble average (a mean value for a simulation run).

In general, $A\geq A_{p}$, and both areas coincide for strictly planar graphene layers, as occurs in a classical calculation at $T=0$. When one takes into account nuclear quantum effects, $A$ and $A_{p}$ are not exactly equal, even for $T\to 0$, because of the zero-point motion of C atoms in the transverse $z$ direction. For graphene monolayers it was found that both areas display temperature dependencies qualitatively different: while $A_{p}$ shows a negative thermal expansion in a wide temperature range, $A$ does not present such a behavior.Zakharchenko et al. (2009); Herrero and Ramírez (2016) This may be different for graphene bilayers, as discussed below. In the following, $A$ and $A_{p}=L_{x}L_{y}/N$ will refer to the real and in-plane area per atom, respectively.

II.3 Mean-square displacements

For each quantum path of a particle (here atomic nucleus), we define the centroid (center of mass) as

[TABLE]

where ${\bf r}_{ij}\equiv(x_{ij},y_{ij},z_{ij})$ is the three-dimensional position of bead $j$ in the ring polymer associated to nucleus $i$. For the out-of-plane motion we consider the $z$-coordinate of the polymer beads. Then, the mean-square displacement $(\Delta z)^{2}_{i}$ of atomic nucleus $i$ ($i=1,...,2N$) in the $z$ direction along a PIMD simulation run is defined as

[TABLE]

Here $\overline{z}_{i}$ is the instantaneous $z$-coordinate of the centroid of atom $i$, and $\left<\overline{z}_{i}\right>$ is its mean value along a simulation run. Hence, $\left<(z_{ij}-\left<\overline{z}_{i}\right>)^{2}\right>$ is the MSD of the z-coordinate $z_{ij}$ of bead $j$ with respect to the average centroid $\left<\overline{z}_{i}\right>$. Then, $(\Delta z)^{2}_{i}$ is the mean of those displacements for the beads corresponding to atomic nucleus $i$ (j = 1, …, $N_{\rm Tr}$).

The kinetic energy of a particle is related to its quantum delocalization, or in the present context, to the spread of the paths associated to it. This can be measured by the mean-square radius-of-gyration of the ring polymers, with an out-of-plane component:Gillan (1988, 1990)

[TABLE]

Note the difference between the r.h.s. of Eqs. (3) and (4): in the former there appears an average of the centroid position over the whole trajectory, $\left<\overline{z}_{i}\right>$, whereas in the latter we have the instantaneous $\overline{z}_{i}$ for each configuration. Then, $Q_{z,i}^{2}$ is the $z$-component of the mean-square ”radius-of-gyration” of the paths corresponding to atom $i$. The total spatial delocalization $(\Delta z)^{2}_{i}$ of atomic nucleus $i$ in the $z$ direction at a finite temperature includes, in addition to $Q_{z,i}^{2}$, another contribution which accounts for classical-like motion of the centroid coordinate $\overline{z}_{i}$, i.e.

[TABLE]

with

[TABLE]

Thus, in our context of the path-integral formulation, the term $C_{z,i}^{2}\textordmasculine$ is the mean-square displacement (MSD) of the centroid of atomic nucleus $i$, and the quantum component $Q_{z,i}^{2}$ is the average MSD of the path (here beads) with respect to the instantaneous centroid. $C_{z,i}^{2}$ behaves as a semiclassical thermal contribution to $(\Delta z)^{2}_{i}$, since at high temperature it converges to the mean-square displacement corresponding to a classical model, and in this limit the quantum paths collapse onto single points ($Q_{z,i}^{2}\to 0$).

In the other limit, for $T\to 0$, $C_{z,i}^{2}$ vanishes and $Q_{z,i}^{2}$ corresponds to zero-point motion of atomic nucleus $i$. In the results presented below, we will show data for $(\Delta z)^{2}$ calculated as an average for $2N$ atoms in the simulation cell:

[TABLE]

and similarly for $Q_{z}^{2}$ and $C_{z}^{2}$.

II.4 Harmonic approximation

For out-of-plane vibrations in graphene bilayers, the behavior of $C_{z}^{2}$ and $Q_{z}^{2}$ may be explained in terms of a harmonic model for the vibrational modes. In a quantum harmonic approximation (HA), the mean-square displacement at temperature $T$ is given by

[TABLE]

where $m$ is the carbon atomic mass, $k_{B}$ is Boltzmann’s constant, and the index $i$ ($i$ = 1, …, 4) refers to the phonon bands with atomic displacements along the $z$ direction (ZA, ZO’, and two ZO bands).Karssemeijer and Fasolino (2011); Yan et al. (2008); Singh and Hennig (2013); Koukaras et al. (2015) The sum in ${\bf k}$ is extended to wavevectors ${\bf k}=(k_{x},k_{y})$ in the 2D hexagonal Brillouin zone, with ${\bf k}$ points spaced by $\Delta k_{x}=2\pi/L_{x}$ and $\Delta k_{y}=2\pi/L_{y}$.Ramírez et al. (2016) For rising system size $N$, there appear vibrational modes with longer wavelength $\lambda$. We have an effective wavelength cut-off $\lambda_{max}\approx L$, with $L=(NA_{p})^{1/2}$, so that the minimum wavevector accessible is $k_{0}=2\pi/\lambda_{max}$, which scales as $k_{0}\sim N^{-1/2}$.

The classical contribution to the MSD at temperature $T$ is given by

[TABLE]

This expression is accurate for classical motion at relatively low temperatures, and at high $T$ anharmonicity causes $C_{z}^{2}$ to increase sublinearly with $T$, as shown below. From $(\Delta z)^{2}$ and $C_{z}^{2}$, we calculate $Q_{z}^{2}$ as the difference $Q_{z}^{2}=(\Delta z)^{2}-C_{z}^{2}$. For $T\to 0$ the MSD converges in the harmonic approximation to

[TABLE]

For the calculations presented below, based on the HA, we have used the vibrational frequencies of the four phonon branches in the $z$ direction, derived from a diagonalization of the dynamical matrix corresponding to the LCBOPII potential employed here. The main difference with the phonon bands in monolayer graphene is the appearance of the layer-breathing ZO’ band, which is nearly flat in the 2D ${\bf k}$-space region close to the $\Gamma$ point ($|{\bf k}|=k=0$), with a frequency for small $k$: $\omega_{0}$ = 92 cm*-1*. This value is close to that obtained earlier from ab-initio calculations for graphene bilayers.Yan et al. (2008)

III In-plane and real area

For an isolated graphene layer we find in a classical calculation in the $T=0$ limit a planar sheet with an interatomic distance $d_{\rm C-C}$ = 1.4199 Å. In a quantum approach, the low-temperature limit includes out-of-plane zero-point motion, so that the graphene sheet is not strictly plane even at $T=0$. In addition to this, anharmonicity of in-plane vibrations gives rise to a zero-point bond expansion, yielding an interatomic distance of 1.4287 Å, i.e., an increase of $8.8\times 10^{-3}$ Å with respect to the classical value.Herrero and Ramírez (2016)

For the graphene bilayer we find at the minimum-energy configuration $d_{\rm C-C}$ = 1.4193 Å, a little smaller than for the monolayer, similarly to the results presented in Ref. Zakharchenko et al., 2010. Our PIMD simulations give for the bilayer a low-temperature C–C distance of 1.4281 Å, i.e., the same zero-point expansion of the interatomic distance as that found for the graphene monolayer. Even if this increase in bond length may seem small, it is much larger than the precision reached for determining cell parameters from diffraction techniques.Yamanaka et al. (1994); Ramdas et al. (1993); Kazimorov et al. (1998) The bond expansion due to nuclear quantum effects decreases as temperature is raised, as in the high-$T$ limit the classical and quantum predictions have to converge one to the other. However, at temperatures so high as 2000 K the value of $d_{\rm C-C}$ derived from PIMD simulations is still distinctly larger than the classical prediction.

As indicated in Sec. II, in the isothermal-isobaric ensemble employed here we allow for fluctuations in the size of the simulation cell in the $xy$ plane ($L_{x}$ and $L_{y}$ lengths). This means that the in-plane area of the graphene sheets changes along a simulation run, controlled by the external in-plane stress $P_{xy}$ (here $P_{xy}=0$). Moreover, C atoms can move in the $z$ coordinate (out-of-plane direction), so that in general any measure of the real surface of a graphene sheet at $T>0$ will give a value larger than the in-plane area.

In Fig. 1 we present the temperature dependence of the in-plane area $A_{p}$ and real area $A$ of a graphene bilayer, as derived from both classical molecular dynamics (open circles) and PIMD simulations (solid squares). These results were obtained for a simulation cell with $N=960$. We discuss first the results of the classical simulations. For the in-plane area $A_{p}$, one observes a slow increase at low $T$, with a slope that becomes larger as the temperature is raised. For larger simulation cells, the low-temperature slope is found to be smaller, and for the largest cells considered here ($N=8400$), it is still positive in the whole temperature range. This contrasts with the results obtained in Ref. Zakharchenko et al., 2010 from classical Monte Carlo simulations, where a slight decrease of the in-plane area was found at low temperature. This is due to the difference in the interatomic potentials, which yields in our case a low-temperature bending constant $\kappa=1.49$ eV vs a value of 0.82 eV derived from the earlier version of the potential for a graphene sheet.Ramírez et al. (2016) A smaller $\kappa$ causes an enhancement of out-of-plane vibrations, which favors a decrease in $A_{p}$ (see below).

For the classical results, we find that both surfaces $A$ and $A_{p}$ converge one to the other in the low-$T$ limit, as expected for a strictly planar configuration when out-of-plane atomic displacements vanish as $T\to 0$. In this classical limit, $A$ and $A_{p}$ go to a value of 2.6169 Å2/atom, which is a little smaller than that obtained for an isolated graphene sheet, i.e. 2.6190 Å2/atom. The real area $A$ derived from classical molecular dynamics simulations increases as a function of temperature much faster than the in-plane area $A_{p}$, and in fact the expansion of $A$ from $T=0$ to 2000 K is a factor of 3.5 larger than that of $A_{p}$.

In Fig. 1 we also show the temperature dependence of $A$ and $A_{p}$ derived from PIMD simulations. This dependence is qualitatively different from that obtained in the classical simulations. For the real area we find a behavior similar to that traditionally observed for the crystal volume in most three-dimensional solids,Kittel (1966) i.e., a rather flat region at low $T$ and a smooth increase for higher temperatures. This is in fact an anharmonic effect, which in 2D materials such as graphene is combined with out-of-plane displacements that grow with rising temperature and contribute to increase the interatomic distance $d_{\rm C-C}$. At low $T$, $A$ converges to 2.6438 Å2/atom. This means a zero-point expansion of $2.5\times 10^{-2}$ Å2/atom with respect to the classical minimum. This difference between classical and quantum results decreases for rising $T$, since nuclear quantum effects become less important.

In contrast to the classical results, the in-plane area derived from PIMD simulations decreases for increasing $T$ in a wide temperature range and then it increases at higher temperatures after reaching a minimum at $T_{m}\approx 800$ K. In both cases, classical and quantum, the difference between real and projected areas increases with temperature. In fact, $A_{p}$ is a measure of a 2D projection of the real graphene surface, and ripples of the actual surface have larger amplitudes as temperature is raised. In the quantum simulations, the in-plane area $A_{p}$ converges at low $T$ to 2.6388 Å2/atom, a value somewhat smaller than that corresponding to the real area. For $T\to 0$, both derivatives $dA/dT$ and $dA_{p}/dT$ converge to zero, as should happen according to the third law of Thermodynamics. Note that this does not happen in the classical approach, a well-known physical inconsistency for this kind of calculations at low temperature.Kittel (1966); Ashcroft and Mermin (1976)

There appear two competing effects which explain the temperature dependence of $A_{p}$. First, the real area $A$ grows as temperature increases in both cases, classical and quantum. Second, the presence of ripples in the graphene sheets gives rise to a decrease in their 2D projection, i.e., in $A_{p}$. At low temperature, the reduction associated to out-of-plane motion in the quantum approach predominates over the thermal expansion of the actual surface, and consequently we have $dA_{p}/dT<0$. In the classical calculation, however, motion in the $z$ direction at low temperature is not enough to compensate for the increase in the area $A$, so that $dA_{p}/dT>0$. At high $T$, the growth of $A$ dominates the decrease in projected area caused by out-of-plane atomic motion. This is in line with the analysis presented by Gao and HuangGao and Huang (2014) for the thermal behavior of $A_{p}$ observed in classical molecular dynamics simulations of a graphene monolayer, and with the calculations by Michel et al.Michel et al. (2015)

It has been shown earlier from the results of classical and quantum simulationsHerrero and Ramírez (2016); Ramírez and Herrero (2017); Gao and Huang (2014) that the size of the simulation cell employed in the calculations can appreciably affect some magnitudes of graphene monolayers, such as the in-plane area per atom $A_{p}$. In Fig. 2 we present the areas $A$ and $A_{p}$ vs temperature, as derived from PIMD simulations of graphene bilayers with several sizes of the simulation cell. Data for $A_{p}$ correspond, from top to bottom, to $N$ = 240 (triangles), 960 (squares), 3840 (diamonds), and 8400 (circles). The latter size is similar to that employed in earlier classical Monte Carlo simulations of graphene bilayers by Zakharchenko at al.Zakharchenko et al. (2010)

In the results of our quantum simulations one observes that $A_{p}(T)$ displays a minimum in all considered cases. This minimum becomes deeper and slightly shifts to higher temperature as the system size increases, converging to a value $T_{m}=850(\pm 50)$ K for the largest cells discussed here. For the real area $A$ (open circles) the size effect is very small, and in fact it is unobservable at the scale of Fig. 2. $A$ and $A_{p}$ become closer to each other as temperature is lowered, but in the low-$T$ limit $A$ is still larger than $A_{p}$, and we find a difference $A-A_{p}=5\times 10^{-3}$ Å2/atom for $T\to 0$. We note that, in spite of the differences in the in-plane area per atom for the different system sizes, the results for all system sizes converge at low $T$ to a single value. This is due to the fact that the graphene sheet becomes totally planar for $T\to 0$ in the classical case and close to planar in the quantum model (see above).

It is interesting to compare the thermal behavior of $A$ and $A_{p}$ for graphene bilayers with that of monolayer graphene. This is presented in Fig. 3 for the quantum case with $N$ = 960. For the real area we find a similar behavior for monolayer and bilayer, with $A$ a little displaced to lower values for the bilayer. The difference between both cases appears as a rigid shift in the temperature region shown in the figure up to 2000 K.

The behavior of the in-plane area $A_{p}$ is somewhat more complex, since at low $T$, $A_{p}$ is larger for the monolayer than for the bilayer (as the real area $A$), but they become equal at $T\approx 200$ K, and at higher temperatures $A_{p}$ for the bilayer is clearly larger. This is a consequence of the competition between the thermal expansion of the real area and the contraction of the projected in-plane area associated to out-of-plane atomic displacements, as indicated above. For the bilayer these displacements (or surface ripples) are smaller (see Sec. VI), and therefore the contraction of $A_{p}$ is smaller.

We have also plotted in Fig. 3 results for $A_{p}$ of graphite derived from PIMD simulations with $N=960$. These results follow the trend displayed when passing from the monolayer to the bilayer, i.e., at low temperature $A_{p}$ for graphite is smaller than that corresponding to bilayer graphene, and becomes larger than the latter for $T>300$ K. The tendency of the in-plane area of graphite to increase with respect to the monolayer at low temperature coincides with the results obtained from density-functional perturbation theory by Mounet and Marzari.Mounet and Marzari (2005) For graphite we obtained a real area $A$ slightly smaller than that of the bilayer. It converges to 2.6419 Å2/atom at low temperature and increases parallel to the real area of the bilayer for rising temperature (not shown in Fig. 3 for simplicity).

IV Internal energy

In this section we present and discuss the contributions to the internal energy of graphene bilayers, corresponding to our isothermal-isobaric ensemble with $P_{xy}=0$ and different temperatures $T$. In Fig. 4 we show the internal energy (kinetic plus potential energy) of a graphene bilayer as a function of temperature, derived from our PIMD simulations for $N=960$ (solid circles). For comparison, we also display results of classical simulations for the bilayer (solid squares), as well as quantum and classical data for a graphene monolayer (open circles and squares). For the zero of energy we have taken the energy $E_{0}$ of an isolated flat graphene sheet, i.e., the energy minimum in a classical calculation at $T=0$. One first observes that the zero-temperature classical limit for the bilayer is shifted by $E_{\rm int}^{0}=-25$ meV/atom, associated to the stabilization energy due to the interaction between layers. This value agrees with that given by Zakharchenko at al.Zakharchenko et al. (2010) from classical Monte Carlo simulations of graphene bilayers, and lies in the intermediate range of binding energies derived from different ab-initio calculations for the AB stacking of the bilayer.Mostaani et al. (2015)

From PIMD simulations of the bilayer, we find that the internal energy $E-E_{0}$ converges at low $T$ to a value of 147 meV/atom, which translates into a zero-point energy of 172 meV/atom. This value is similar to that found for the monolayer (171 meV/atom), which is not strange since the main contribution to the zero-point energy comes from high-frequency in-plane vibrational modes, which are not affected by the interaction between layers. At high temperature, the energy results derived from PIMD simulations converge to those of classical simulations, but at $T=1500$ K we still observe an appreciable difference between classical and quantum results for both monolayer and bilayer.

The results shown in Fig. 4 correspond to $N=960$, as indicated above. For other cell sizes we obtained results for the internal energy very close to those shown in the figure, and in fact indistinguishable from them at the scale of Fig. 4. This does not happen for other properties of the bilayer, as indicated below.

To obtain insight into the origin of the changes in $E$ as a function of temperature, we may decompose the internal energy in different contributions as:

[TABLE]

where $E_{\rm el}$ is the elastic energy corresponding to an area $A$, $E_{\rm vib}$ is the vibrational energy of the system, and $E_{\rm int}$ is the stabilization energy due to the interaction between layers. All three contributions change with the temperature. The expression for the internal energy of the bilayer in Eq. (11) is similar to that corresponding to an isolated graphene monolayer,Herrero and Ramírez (2016) the only difference being the absence of the term $E_{\rm int}$ in the latter case. This term corresponds to the binding energy of the graphene layers. At low temperature $E_{\rm int}$ is indistinguishable from the classical minimum $E_{\rm int}^{0}=-25$ meV/atom, and it slightly changes as a function of the interlayer distance, which increases slowly for rising $T$ (see below). At the highest temperature considered here ($T=2000$ K), the interlayer distance amounts to 3.575 Å, which gives $E_{\rm int}=-23.1$ meV/atom (see Table I).

PIMD simulations directly give $E(T)$, and the elastic energy can be calculated from the resulting real area $A$. To obtain the elastic energy corresponding to a given area $A$, we have calculated the classical energy of a flat graphene sheet with that area, obtained by isotropically expanding or contracting the minimum-energy configuration in the layer plane.Herrero and Ramírez (2016) For small changes in the real area, $E_{\rm el}$ is found to change as $E_{\rm el}(A)\approx K(A-A_{0})^{2}$, with $K$ = 2.41 eV/Å2. At low temperatures most of the energy $E-E_{0}$ corresponds to the vibrational energy, and the contribution of the elastic energy increases as the real area grows up.

For a graphene monolayer, the elastic energy $E_{\rm el}$ corresponding to the area $A$ derived from PIMD simulations at 300 K amounts to 1.5% of the internal energy $E-E_{0}$, most of this energy corresponding to the vibrational energy $E_{\rm vib}$. For the graphene bilayer at 300 K, we have $E_{\rm el}$ = 1.6 meV/atom, to be compared with the internal energy $E-E_{0}$ = 159.4 meV/atom (see Table I). Then, in this case $E_{\rm el}$ amounts to 1.0% of the internal energy. This value increases to a 2.6% at 2000 K, since at this temperature the elastic and internal energies take values of 13.5 and 525.1 meV/atom, respectively. We find a reduction of the relative contribution of the elastic energy to the internal energy of the bilayer, as compared to monolayer graphene. This is mainly caused by the smaller value of the real area $A$ for the bilayer (see Fig. 3).

We calculate the vibrational energy $E_{\rm vib}$ from the results of PIMD simulations by subtracting the elastic and interaction energy, $E_{\rm el}$ and $E_{\rm int}$, from the internal energy $E$ at each temperature [see Eq. (11)]. At 50 and 300 K, $E_{\rm vib}$ amounts to 171 and 183 meV/atom, respectively (see Table I). The latter value is somewhat smaller than those found for diamond from path-integral simulations at 300 K, namely 195 and 210 meV/atom, obtained using Tersoff-type and tight-binding potentials, respectively.Herrero and Ramírez (2000); Ramírez et al. (2006) For comparison, in a classical harmonic approximation one has: $E_{\rm vib}^{\rm cl}(T)=3k_{B}T$ = 77.6 meV/atom.

V Interlayer spacing and compressibility

In this section we study the interlayer distance, $c$, and its fluctuations in bilayer graphene. In Fig. 5 we present the temperature dependence of the mean equilibrium distance $\langle c\rangle$, as derived from our PIMD simulations (solid circles). For comparison, we also show the results of classical molecular dynamics simulations (open circles). The classical data converge at low $T$ to an interlayer spacing $c_{0}$ = 3.3372 Å, which corresponds to the minimum-energy configuration, i.e., totally planar graphene sheets in AB stacking. These classical calculations yield a linear increase of the mean interlayer distance $\langle c\rangle$ for rising temperature, with a derivative $d\langle c\rangle/dT=1.2\times 10^{-4}$ Å/K. This linear interlayer expansion is similar to that found for lattice parameters of crystalline solids in a classical approximation, which is known to violate the third law of thermodynamics at low temperature,Callen (1960) since thermal expansion coefficients should vanish for $T\to 0$. This anomaly of the classical model is remedied in the quantum simulations, which yield a vanishing derivative $d\langle c\rangle/dT$ in the low-temperature limit, similarly to the behavior of real and in-plane areas (see Sec. III).

PIMD simulations give an interlayer spacing larger than classical calculations, mainly due to zero-point motion of the C atoms in the quantum model, which detects anharmonicities in the interatomic potential even at low temperature. For $T\to 0$ the interlayer distance converges to 3.3521 Å. Thus, we find a zero-point expansion of $1.5\times 10^{-2}$ Å, i.e., the mean spacing between layers increases by about 0.5% with respect to the classical prediction. At room temperature ($T$ = 300 K), the difference between classical and quantum results amounts to $3.2\times 10^{-3}$ Å, about five times less than in the low-temperature limit. Size effects due to the finite simulation cells are negligible for the interlayer spacing. In fact, for a given temperature, we did not find any difference between the results for $\langle c\rangle$ obtained for the considered cell sizes, i.e., differences were in the order of the error bars found for each cell size (less than the symbol size in Fig. 5).

In Fig. 5 we also display the interlayer spacing of graphite derived from PIMD simulations (solid squares). At low temperature these results converge to a value very close to the mean distance $\langle c\rangle$ for the bilayer. For rising $T$, the mean spacing $\langle c\rangle$ for the bilayer becomes progressively larger than that of graphite, which agrees with larger thermal fluctuations of $c$ in the bilayer (see below). Diamonds in Fig. 5 represent data points derived for graphite from x-ray diffraction experiments.Baskin and Meyer (1955) Comparing these experimental results with those of our PIMD simulations for graphite, we find that the interatomic potential LCBOPII employed here overestimates the interlayer spacing by a 0.5%.

From the analysis of the interlayer spacing and its fluctuations one can study the compressibility of bilayer graphene in the out-of-plane direction. The isothermal compressibility $\chi_{z}$ in the $z$ direction is defined as

[TABLE]

where $V=L_{x}L_{y}c$ and $P_{z}$ is a uniaxial stress in the $z$ direction. At a finite temperature $T$, $\chi_{z}$ of bilayer graphene can be conveniently calculated from PIMD simulations with $P_{z}=0$ by using the fluctuation formulaLandau and Lifshitz (1980); Herrero (2008)

[TABLE]

where the volume mean-square fluctuations due to changes in the interlayer spacing $c$ are given by $(\Delta V)^{2}=L_{x}^{2}L_{y}^{2}(\Delta c)^{2}$, with $(\Delta c)^{2}=\langle c^{2}\rangle-\langle c\rangle^{2}$. Thus, we calculate here the isothermal compressibility $\chi_{z}$ by using the expression

[TABLE]

In Fig. 6 we display the dependence of $\chi_{z}$ upon temperature, as derived from our PIMD simulations (solid squares). For comparison we also present data for $\chi_{z}$ obtained from classical simulations using Eq. (14) (open symbols). The classical compressibility $\chi_{0}$ for $T\to 0$ can be calculated from the dependence of the system energy on the interlayer spacing $c$ close to $c_{0}$. For small variations in $c$ we have an interaction energy $E_{\rm int}=E_{\rm int}^{0}+k(c-c_{0})^{2}/2$, with a constant $k$ = 0.093 eV Å*-2*/atom. Then, the classical zero-temperature compressibility is given by $\chi_{0}=A_{0}/2c_{0}k$, where $A_{0}$ is the area per atom (see Appendix A). We find $\chi_{0}=2.63\times 10^{-12}$ cm2 dyn*-1* or 0.0263 GPa*-1*. Our classical results converge at low $T$ to this limit, which is indicated in Fig. 6 with a horizontal arrow.

The compressibilities derived from PIMD simulations at $T\gtrsim 300$ K are slightly higher that the classical results, and at lower temperatures they depart progressively one from the other. At low $T$ the quantum results converge to a value of $2.79(2)\times 10^{-2}$ GPa*-1*. This means an appreciable increase in $\chi_{0}$ of a 6% with respect to the classical value $\chi_{0}$.

In connection with the interlayer coupling, the main difference between the phonon spectrum of monolayer and bilayer graphene is the appearance in the latter of the so-called ZO’ vibrational band, which is nearly flat in a region of 2D ${\bf k}$-space close to the $\Gamma$ point. The frequency of this band at the $\Gamma$ point ($\omega_{\rm ZO^{\prime}}^{\Gamma}$, called here $\omega_{0}$ for simplicity) is related to the coupling constant $k$ as $\omega_{0}=(k_{N}/M_{\rm red})^{1/2}$, with $k_{N}=2Nk$ and the reduced mass $M_{\rm red}=\frac{1}{2}Nm$ ($m$: carbon atomic mass), so that $\omega_{0}=2(k/m)^{1/2}$. Using the coupling constant $k$ = 0.093 eV Å*-2*/atom we find a frequency $\omega_{0}$ = 92 cm*-1*. This corresponds to the layer-breathing $A_{2g}$ Raman-active mode, for which a frequency of 89 cm*-1* has been reported.Lin et al. (2018)

The interlayer coupling has been studied earlier from classical Monte Carlo simulations in Ref. Zakharchenko et al., 2010. These authors employed a parameter $\gamma$ to describe the low-frequency part of the ZO’ band, which is given in our terminology by $\gamma=\rho\,\omega_{0}^{2}/4$, $\rho$ being the surface mass density. This gives a value $\gamma$ = 0.035 eV Å*-4*, in agreement with the low-temperature results of Monte Carlo simulations.Zakharchenko et al. (2010)

For comparison with the results for the graphene bilayer, we also present in Fig. 6 some data points for the compressibility of graphite, as derived from PIMD simulations and the fluctuation formulas presented above. At low $T$ these results converge within error bars to the same value as the bilayer compressibility, since in both cases the MSD $(\Delta c)^{2}$ are nearly identical. At higher temperatures, however, $(\Delta c)^{2}$ is smaller for graphite, so that its compressibility is lower than that of bilayer graphene.

It is important to note that the compressibility $\chi_{z}$ defined here coincides in the case of graphite with the inverse of the elastic constant $C_{33}$ of this material, as this constant relates strain and stress in the $z$ direction. In Fig. 6 we present the inverse elastic constant $1/C_{33}$ determined for pyrolytic graphite from ultrasonic test methodsBlakslee et al. (1970) (triangle up) and from neutron diffraction results combined with a force modelNicklow et al. (1972) (triangle down). These data, obtained at room temperature, are slightly displaced horizontally around 300 K for clarity. The compressibility $\chi_{z}$ for graphite is overestimated by our simulation results about a 5% with respect to these data derived form experiment. For natural graphite, KomatsuKomatsu (1964) obtained a low-temperature value $C_{33}=3.56\times 10^{11}$ dyn/cm2 from specific-heat measurements, which corresponds to $\chi_{z}$ = 0.0282 GPa*-1*. This value (not shown in the figure) is closer to our results for graphite, but there is no available error bar for it.

VI Out-of-plane motion

Graphene, as a 2D material in three-dimensional space, has peculiar properties due to out-of-plane motion, such as the decrease in in-plane area $A_{p}$ for rising temperature discussed in Sec. III. For monolayer graphene, motion in the $z$ direction has been discussed earlier as a function of temperature, applied stress, and system size.Ramírez et al. (2016); Ramírez and Herrero (2017); Gao and Huang (2014) In this context, one expects that nuclear quantum effects should be important at relatively low temperatures, due to the low frequency of vibrational modes in the $z$ direction, especially those corresponding to the ZA flexural band. These effects have been studied earlier for a graphene monolayer, and in particular the competition between classical thermal motion and quantum delocalization. In that case it was found that the outcome of this competition is not trivial, as it depends significantly not only on temperature, but also on the applied stress and system size.Herrero and Ramírez (2016, 2017)

In this section we analyze the mean-square displacements of C atoms in the out-of-plane direction, as obtained from PIMD simulations of bilayer graphene. We will mainly focus on the character of these atomic displacements, in order to find whether they can be well explained by a classical model, or the carbon atoms noticeably behave as quantum particles. We put particular emphasis on the question whether the system size $N$ plays or not a relevant role in this problem. In this line, PIMD simulations allow us to split vibrational amplitudes or atomic delocalization into two parts: one component associated to thermal (classical-like) motion and another corresponding to a proper quantum contribution, which can be quantified by the mean size of the quantum paths at a given temperature (see Sec. II.C).

We will first present results for the MSD obtained in classical simulations, where the statistics can be more readily improved because less computational resources are required. Quantum effects are expected to be noticeable mainly for $T$ lower than room temperature, due to the relatively low frequency of the vibrational modes that give the main contribution to the MSD in the out-of-plane direction. Note that this does not happen for in-plane modes, with an average frequency larger than out-of-plane vibrations, so that their associated zero-point energy is clearly larger.

In Fig. 7 we display results for the MSD in the $z$ direction for C atoms in bilayer graphene, as derived for classical molecular dynamics simulations for three system sizes (solid symbols). From top to bottom: $N$ = 3840 (diamonds), 960 (squares), and 240 (circles). For comparison we also show results for monolayer graphene with $N$ = 960 (open squares, labeled as “ML”). The vibrational amplitude in the $z$ direction increases with system size, and displays an appreciable anharmonic effect, since the data obtained from the simulations largely depart from linearity. In fact, in a classical harmonic approximation one expects for a given system size a MSD increasing linearly with temperature, as shown by the solid line corresponding to $N=3840$. From the results for monolayer and bilayer with $N=960$, we observe an appreciable decrease in the MSD for the bilayer, as compared with the monolayer. The ratio between MSD of monolayer and bilayer is 2.3 at 50 K and decreases for rising $T$, taking a value of $\approx 1.6$ in the region between 1000 and 2000 K. The main difference between monolayer and bilayer in this respect is the less relative importance of ZA modes in the bilayer, due to the appearance in this case of ZO’ and two ZO vibrational bands.

Let us now turn to the results of our quantum simulations. In Fig. 8 we display data for the motion in the out-of-plane direction, corresponding to a cell with $N$ = 240. Shown are data for the MSD $(\Delta z)^{2}$ obtained for graphene monolayers (open symbols) and bilayers (solid symbols). In both cases we present results from classical (squares, labeled as “cl”) and PIMD simulations (circles, labeled as “q”). At high $T$, classical and quantum data converge one to the other in both cases, molayer and bilayer graphene. One observes that the MSD derived from classical simulations goes to zero in the low-temperature limit in each case, whereas PIMD simulations yield a finite MSD caused by zero-point motion. We find $(\Delta z)^{2}_{0}$ = 5.9 and $4.9\times 10^{-3}$ Å2 for the monolayer and bilayer, respectively. The lower zero-point vibrational amplitude for the bilayer is in line with an average increase in the frequency of out-of-plane modes caused by interlayer interactions (less relative importance of the ZA phonon band). Comparing the quantum results for monolayer and bilayer at different temperatures, we find a ratio between MSD of monolayer and bilayer which increases from 1.2 at low $T$ to 1.6 at high temperature, as in the classical limit. For comparison with the data found from our simulations, we also present in Fig. 8 the result of a quantum HA (solid line), as derived from Eq. (8). At low temperature, the results of HA and PIMD simulations coincide, and progressively depart one from the other as temperature is raised.

The image displayed in Fig. 8 for the atom displacements in the $z$ direction is qualitatively similar for different system sizes $N$. The main difference appears in the relative contribution of $Q_{z}^{2}$ and $C_{z}^{2}$ to the total MSD $(\Delta z)^{2}$. This is caused by the enhancement of the classical MSD for increasing $N$, a fact observed also in data derived from classical molecular dynamics simulations of graphene monolayers.Gao and Huang (2014); Ramírez et al. (2016)

As noted above, an interesting point that can be studied from PIMD simulations is the competition between classical-like and quantum motion as a function of temperature and system size. One expects that quantum motion should dominate at relatively low temperatures, but this turns out to be highly dependent on the size $N$. In Fig. 9 we present results for $C_{z}^{2}$ and $Q_{z}^{2}$ as a function of temperature for several system sizes, from $N$ = 24 to 960. Data points for the quantum delocalization $Q_{z}^{2}$ for different cell sizes are indistinghuisable at the scale of the figure. Clear size effects are only found for very small simulation cells ($N<20$). On the contrary, the results for $C_{z}^{2}$ change very much with system size. In the temperature region shown in Fig. 9 the classical MSD $C_{z}^{2}$ is nearly linear with $T$ for the considered system sizes. The values of $C_{z}^{2}$ presented here, corresponding to the MSD of the path centroids in PIMD simulations, coincide within error bars with the atomic MSD obtained from classical simulations. Note that at higher $T$, $C_{z}^{2}$ clearly departs from linearity, as shown in Fig. 7. For comparison we also show in Fig. 9 results of the classical HA given by Eq. (9) for $N=240$ (dashed-dotted line). From the data presented in this figure, it is also worthwhile noting that for system sizes $N<1000$, $C_{z}^{2}$ obtained from the simulations scales as $N^{\epsilon}$, with an exponent $\epsilon=0.67$.

Given the increasing slope $dC_{z}^{2}/dT$ for rising $N$, the crossing of the curves corresponding to $C_{z}^{2}$ and $Q_{z}^{2}$ moves to lower temperatures. This means that for large $N$, classical-like motion becomes dominant over quantum delocalization for displacemtns in the $z$ direction. The quantum contribution $Q_{z}^{2}$ converges for $T\to 0$ to the value given above $(\Delta z)^{2}_{0}=5\times 10^{-3}$ Å2 for the zero-point delocalization in the $z$ direction. The nuclear quantum delocalization may be estimated from the mean extension of the quantum paths, i.e., from $Q_{z}^{2}$. At 12.5, 50, and 300 K we find an average extension $(\Delta z)_{Q}=(Q_{z}^{2})^{1/2}$ = 0.066, 0.055, and 0.032 Å, respectively.

For a system size $N$, the ratio $Q_{z}^{2}/C_{z}^{2}$ decreases for increasing $T$, and there appears a crossover temperature $T_{c}$ for which this ratio equals unity. This corresponds to the crossing of the lines of classical and quantum displacements, $C_{z}^{2}$ and $Q_{z}^{2}$, in Fig. 9. For temperatures $T>T_{c}$ classical-like motion is dominant for the atomic motion in the $z$ direction. In Fig. 10 we display $T_{c}$ as a function of the system size $N$, as calculated from PIMD simulations of bilayer graphene. For comparison, we also present data for a graphene monolayer, as well as for graphite. Symbols are results derived from our simulations and dashed lines are guides to the eye. For a given system size $N$, the crossover temperature for the bilayer is higher than that corresponding to monolayer graphene, and lower than $T_{c}$ for graphite. This means that the quantum behavior of out-of-plane motion is more relevant for the bilayer than for an isolated monolayer. It is even more aprreciable for graphite. For small $N$ we find similar values of $T_{c}$ in the three cases, and the difference between them increases as the system size is raised. For $N=960$, we have $T_{c}$ = $10.1(\pm 0.6)$ K for the monolayer, $17.8(\pm 0.8)$ K for the bilayer, and $32.3(\pm 0.9)$ K for graphite. Thus, for this system size, $T_{c}$ for the bilayer is a factor 1.8 higher than for the monolayer.

On the other side, for a given temperature $T$ the ratio $C_{z}^{2}/Q_{z}^{2}$ grows for increasing $N$, and there is a system size $N_{c}$ for which classical motion along the $z$ direction becomes dominant over quantum delocalization. The origin of this behavior is the following. For a temperature $T$, vibrational modes with frequency $\omega<\omega_{c}(T)=k_{B}T/\hbar$ can be considered in the classical regime. The main contribution to $C_{z}^{2}$ and $Q_{z}^{2}$ comes from ZA flexural modes, whose frequency behaves for small $k$ as $\rho\,\omega_{\rm ZA}^{2}=\sigma k^{2}+\kappa k^{4}$, and $\omega_{\rm ZA}$ decreases as $k$ is reduced ($\rho$, surface mass density; $\sigma$, effective stress; $\kappa$, bending modulus).Ramírez et al. (2016) For a size $N$ we have an effective minimum wavevector $k_{0}=2\pi/\lambda_{max}\approx 2\pi/(NA_{p})^{1/2}$. Then, increasing $N$ we attain a size $N_{c}(T)$ for which additional increase in size introduces new modes (ZA modes in particular) with frequency $\omega<\omega_{c}$, contributing to rise $C_{z}^{2}$ more than $Q_{z}^{2}$. Hence, for any $T>0$ classical-like motion dominates over quantum delocalization, when the system size is larger than the corresponding $N_{c}$.

For $N>100$ the data points corresponding to the bilayer in Fig. 10 can be fitted according to a power-law dependence of the crossover temperature on system size, i..e, $T_{c}=b\,N^{-a}$ with a prefactor $b=735$ K and an exponent $a=0.54$. Extrapolating this dependence to low temperatures, we have for $T$ = 1 K and 0.1 K crossover sizes $N_{c}\approx 2\times 10^{5}$ and $1.4\times 10^{7}$, respectively. These sizes are much larger than those actually manageable in our simulations at these temperatures.

The solid line in Fig. 10 corresponds to the HA model which takes into account the four vibrational bands of bilayer graphene in the $z$ direction. For $N\gtrsim 100$, this model predicts a dependence of the crossover temperature $T_{c}$ upon $N$ similar to that derived from PIMD simulations. It yields $T_{c}$ values smaller than the simulations, but for large system size we find a power law dependence with the same exponent in both cases, i.e., $T_{c}\sim N^{-a}$ with $a=0.54(\pm 0.02)$. The main limitations of the HA discussed here are the neglect of anharmonicity in the out-of-plane vibrational modes (expected to be reasonably small at low $T$) and coupling with in-plane modes, which is also expected to increase for rising $T$. Taking into account these limitations, the harmonic approach is still able to reproduce qualitatively, and nearly quantitatively, the main features of the competition between classical-like and quantum dynamics of C atoms in the out-of-plane motion.

VII Summary

We have presented results of PIMD simulations of graphene bilayers at zero in-plane stress in a wide range of temperatures. The importance of nuclear quantum effects has been assessed by comparing the results of these simulations with those obtained from classical molecular dynamics simulations. Structural variables such as interatomic distances, as well as in-plane and real areas are found to increase when quantum nuclear motion is considered, even at temperatures in the order of 1000 K. Zero-point expansion of the graphene layers amounts to about 1% of the areas $A$ and $A_{p}$. This is much larger than the precision currently reached by diffraction techniques in determining structural parameters. The characteristic behavior of the in-plane area $A_{p}$ (decreasing at low $T$ and increasing at high $T$) is a result of the coupling between in-plane and out-of-plane modes. Changes in $A_{p}$ are , however, less than those obtained for isolated graphene layers, as a consequence of interlayer interactions.

We have put particular emphasis on the vibrational motion in the out-of-plane direction. Quantum effects appear in these vibrational modes at low temperatures, but thermal classical motion becomes dominant for large system size. This important size effect is a consequence of the enhancement of classical-like displacements, whereas quantum delocalization is nearly unaffected by the system size. As a result, the crossover temperature $T_{c}$ at which classical motion becomes dominant scales for large size as $T_{c}\sim N^{-a}$ with an exponent $a=0.54$.

Anharmonicity of the vibrational modes in the $z$ direction is clearly observable in the MSDs presented in Figs. 7, 8, and 9. This affects to both classical and quantum results. Such anharmonicity shows up markedly in the temperature dependence of the in-plane area $A_{p}$. The real area $A$ is basically determined by the interatomic distance $d_{\rm C-C}$, so its change with temperature is caused by anharmonicity of the in-plane modes. In this context, it is noteworthy that a pure harmonic approximation yields rather accurate results for MSDs at low temperatures, once the vibrational frequencies have been calculated for the classical equilibrium geometry of the bilayer at $T=0$.

The compressibility $\chi_{z}$ of graphene bilayers in the out-of-plane direction has been obtained from the fluctuations of the interlayer distance. This method has turned out to be accurate enough to follow the increase in $\chi_{z}$ for rising $T$, and also to assess the importance of quantum effects at low temperature. For $T\to 0$ we find an increase in $\chi_{z}$ of a 6% with respect to the classical limit.

A further check of our finite-temperature results would consist in studying structural and thermodynamic properties of graphene multilayers from an ab-initio method. This is, however, not yet possible taking into account the relatively large size of the supercells required to study these properties and the length of the trajectories necessary for a low statistical uncertainty.

Path-integral simulations analogous to those presented in this paper could contribute to understand structural and dynamical properties of light-atom monolayers on graphene. This is the case of graphane, where nontrivial quantum features can appear at low temperature, associated to the small mass of hydrogen.

Acknowledgements.

The authors acknowledge the help of J. H. Los in the implementation of the LCBOPII potential. This work was supported by Ministerio de Ciencia, Innovación y Universidades (Spain) through Grant FIS2015-64222-C2.

Appendix A Compressibility

For an interlayer distance $c$ close to the minimum-energy distance $c_{0}$, the interaction energy per atom in bilayer graphene can be written as

[TABLE]

with $k$ = 0.093 eV/Å2. The classical zero-temperature compressibility in the $z$ direction is given by

[TABLE]

where $P_{z}$ is a uniaxial stress in the $z$-direction, $\bar{E}_{\rm int}$ is the interaction energy per simulation cell, $\bar{E}_{\rm int}=2NE_{\rm int}$, and changes in $V=NA_{p}c$ are associated to the variation of interlayer spacing $c$. We find

[TABLE]

where $A_{0}$ is the area per atom at $T=0$ (minimum-energy configuration).

The isothermal compressibility $\chi_{z}$ at temperature $T$ can be calculated from the fluctuation formula:

[TABLE]

with the volume mean-square fluctuations $(\Delta V)^{2}=N^{2}A_{p}^{2}(\Delta c)^{2}$.

In a harmonic approximation, the MSD of the interlayer spacing, $(\Delta c)^{2}$, can be calculated from the cost of energy associated to changes of $c$ for a size $N$ ($2N$ atoms in bilayer graphene): $\bar{E}_{\rm int}=\bar{E}_{\rm int}^{0}+kN(c-c_{0})^{2}$. This yields at temperature $T$ a classical harmonic MSD:

[TABLE]

which introduced into Eq. (18) gives $\chi_{z}=A_{p}/2ck$. This expression for the compressibility is valid for interlayer distances close to $c_{0}$, where $\bar{E}_{\rm int}$ changes as $(c-c_{0})^{2}$. In general, $\chi_{z}$ is affected by anharmonicity and has different values for classical and quantum calculations. The results presented in Sec. V were obtained from our simulations by using the fluctuation formula given in Eq. (18).