Quantum effects in graphene monolayers: Path-integral simulations

In this section we present results for interatomic distances in graphene. The temperature dependence of the equilibrium C–C distance, $d_{\rm C-C}$, as derived from our PIMD simulations at zero applied stress is displayed in Fig. 1 (squares). In the low-temperature limit, we find an interatomic distance of 1.4287(1) Å, which increases for increasing temperature, as could be expected from the thermal expansion of the graphene sheet. We note that the size effect of the finite simulation cell is negligible for our purposes. In fact, for a given temperature, we found no difference between the results for $d_{\rm C-C}$ obtained for the considered cell sizes, i.e., differences between data for $N>200$ were similar to the error bars obtained for each cell size (much less than the symbol size in Fig. 1).

For comparison with the results of the quantum simulations, we also present in Fig. 1 the temperature dependence of $d_{\rm C-C}$ in the classical limit with the same LCBOPII potential, as derived from MD simulations (circles). These results show at low $T$ a nearly linear increase, similar to that found for lattice parameters of crystalline solids in a classical approximation, which is known to violate the third law of thermodynamicsCallen (1960); Wallace (1972) (i.e., the thermal expansion coefficient should vanish for $T\to 0$). This anomaly of the classical model is remedied in the quantum simulations, which yield a vanishing slope for the $d_{\rm C-C}$ vs $T$ plot in the low-temperature limit. The results of the classical simulations converge at low $T$ to an interatomic distance of 1.4199 Å, corresponding to the minimum energy for a flat graphene sheet. This value is close to the distance $d_{\rm C-C}$ derived from ab-initio simulations in the limit $T\to 0$.Pozzo et al. (2011)

The quantum simulations predict an interatomic distance larger than the classical calculation, due basically to the fact that zero-point motion of the carbon atoms in the quantum model detects the anharmonicity of the interatomic potential even for $T\to 0$. In this limit, we find a zero-point expansion of the C–C distance of $8.8\times 10^{-3}$ Å, i.e., the bond distance increases by a 0.6% respect the classical prediction. This is close to a zero-point expansion of 0.53% derived by Brito et al.Brito et al. (2015) from path-integral Monte Carlo simulations. This increase in mean bond length is much larger than the precision currently achieved in the determination of 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 rises, as should happen, since in the high-$T$ limit the classical and quantum predictions have to converge one to the other. We observe, however, that at a temperature of 2000 K the C–C distance obtained from the quantum simulations is still larger than the classical prediction (the difference is much larger than the error bars).

In the quantum results, the increase in interatomic distance from $T$ = 0 K to room temperature ($T$ = 300 K) is small, and amounts to $\sim 5\times 10^{-4}$ Å, more than 15 times smaller than the zero-point expansion. Note that the bond expansion associated to zero-point motion is in the order of the thermal expansion predicted by the classical model from $T$ = 0 to 800 K. These results (classical and quantum) for a single layer of graphene are qualitatively similar to those found earlier from simulations of carbon-based materials. For diamond, in particular, it was found a zero-point expansion of the lattice parameter $\Delta a=1.7\times 10^{-2}$ Å, a 0.5% of the classical prediction.Herrero and Ramírez (2000)

The thermal bond expansion as well as the zero-point bond dilation discussed here are a signature of the anharmonicity of the interatomic potential, similar to 3D crystalline solids (see above). In the case of graphene this effect is basically due to anharmonicity in the stretching vibrations of the sp² C-C bond. A more involved anharmonicity appears in the description of the thermal variation of the graphene area, due to the coupling between in-plane and out-of-plane modes, as shown in the following section.

The simulations (both classical MD and PIMD) presented here were carried out in the isothermal-isobaric ensemble, as explained in Sec. II. This means that in a simulation run we fix the number of atoms $N$, the temperature $T$, and the applied stress in the $(x,y)$ plane ($P=0$ in our simulations), thus allowing changes in the area of the simulation cell on which periodic boundary conditions are applied. Since carbon atoms are free to move in the out-of-plane direction (coordinate $z$), it is obvious that in general any measure of the “real” surface of the graphene sheet will give a value larger than the area of the simulation cell in the $(x,y)$ plane. Taking into account that the actual simulations yield directly the mean in-plane area $A_{\|}$ for given $N$, $T$, and $P$, and that this quantity has been discussed earlier in the literature, we present it for our classical and quantum simulations.

In Fig. 2 we display the temperature dependence of the in-plane area $A_{\|}$ associated to the graphene layer, as derived from (a) classical MD and (b) PIMD simulations. In each plot, results for three cell sizes are given: $N$ = 240 (circles), 960 (squares), and 8400 (diamonds) atoms. Let us look first at the results of the classical simulations in Fig. 2(a). For the larger sizes (960 and 8400 atoms) one observes first a decrease in $A_{\|}$ for rising $T$, and at higher temperatures $A_{\|}$ increases with $T$. For the smallest size presented here, the low-temperature decrease in $A_{\|}$ is also present, although almost unobservable in the figure. Moreover, one sees that the normalized in-plane surface per atom is smaller for larger simulation cells. These results are similar to those found in earlier classical Monte Carlo and MD simulations of graphene single layers.Zakharchenko et al. (2009); Gao and Huang (2014); Brito et al. (2015) There are two competing effects which explain this temperature dependence of $A_{\|}$, as discussed below.

In Fig. 2(b) we present the temperature dependence of $A_{\|}$, as derived from PIMD simulations, for the same cell sizes as in Fig. 2(a). This dependence is qualitatively similar to that corresponding to the classical results, but in the quantum simulations the low-temperature decrease of the in-plane area is more pronounced. This is particularly visible for the smallest size presented here ($N$ = 240), for which such a decrease in $A_{\|}$ for rising $T$ is almost inappreciable in the classical results, whereas it is clearly visible in the results of the quantum simulations. It is important to note that, in spite of the large differences in the in-plane area per atom for the different system sizes, in each case (classical or quantum) all sizes converge at low $T$ to a single value. This could be expected because for $T\to 0$ the graphene sheet becomes strictly planar in the classical case and close to planar in the quantum model (a zero-point effect).

For a direct comparison of the in-plane area $A_{\|}$ obtained from classical and quantum simulations, we show in Fig. 3 both sets of results for system size $N$ = 960. In the limit $T\to 0$, the difference between them converges to 0.022 Ų. This difference decreases for rising temperature, as nuclear quantum effects become less relevant. For this system size, the overall behavior of $A_{\|}$ vs $T$ is similar for both classical and quantum models, in the sense that $dA_{\|}/dT<0$ at low temperature and this derivative becomes positive at high $T$. However, the slope of the $A_{\|}$ vs $T$ curve at low temperature is much larger for the quantum than for the classical model.

In view of the results shown in Fig. 3, it seems that $dA_{\|}/dT$ could not vanish for $T\to 0$, in disagreement with the third law of thermodynamicsCallen (1960); Wallace (1972) (This apparent inconsistency is however not real, as shown below in Sec. IV). Moreover, one can ask if $A_{\|}$ can be treated as a true extensive variable, given that it shows a strong size effect, as shown in Fig. 2. In any case, one can use $A_{\|}$ as a helpful variable to carry out simulations of layered systems such as graphene, and indeed it is $A_{\|}$ which defines the area of the simulation cell in the $(x,y)$ plane. This in-plane area has been studied earlier in several worksGao and Huang (2014); Brito et al. (2015); Zakharchenko et al. (2009); Los et al. (2016); Chacón et al. (2015) as a function of temperature and external stress, but one can argue that it is not necessarily a variable to which one can associate properties of a real material surface.

This can be further illustrated by considering a “real” surface $A$ (in 3D space) for graphene, as can be defined from the actual bonds between carbon atoms. Thus, we define a mean area per atom from the mean C–C distance derived from the simulations, as $A=3\sqrt{3}\,d_{\rm C-C}^{2}/4$ (see Fig. 7 in Ref. Ramírez et al., 2016). A similar expression has been used by Hahn et al.Hahn et al. (2016) to study nanocrystalline graphene from atomistic simulations. As the actual interatomic distance changes from bond to bond in a simulation step and for each bond it fluctuates along a simulation run, $A$ is a mean value for a given temperature, and fluctuations of the instantaneous area depend on the fluctuations of the interatomic distance. One could employ alternative definitions for the real area, based for example on the areas of the hexagons that make up the graphene sheet, or on a triangulation based on the atomic positions. We have checked that such definitions yield area values slightly smaller than the area $A$ defined above, but in the present conditions (zero applied stress and not very large out-of-plane fluctuations) this difference is not relevant for our current purposes.

In Fig. 4 we show the areas $A$ and $A_{\|}$ vs temperature for a simulation cell including 960 atoms. In both cases, we present results from classical and PIMD simulations. The surface $A$ is larger than $A_{\|}$, and the difference between both increases with temperature. In fact, $A_{\|}$ is the 2D projection of $A$, and the actual surface becomes increasingly bent as temperature is raised and out-of-plane atomic displacements are larger. For the area $A$ we do not observe the anomalous decrease displayed by $A_{\|}$ in both classical and quantum cases at low temperature. Moreover, the area $A$ derived from quantum simulations shows a temperature derivative $dA/dT$ approaching zero as $T\to 0$, in agreement with low-temperature thermodynamics.

For the classical results, we find that both surfaces $A$ and $A_{\|}$ converge one to the other in the low-$T$ limit, as expected for a convergence of the actual graphene surface to a plane when out-of-plane atomic displacements go to zero as $T\to 0$. In the quantum case, $A$ and $A_{\|}$ become closer as temperature is reduced, but in the low-$T$ limit $A$ is still larger than $A_{\|}$, and the extrapolated difference $A-A_{\|}$ amounts to $1.1\times 10^{-2}$ Ų/atom for $T=0$. This difference is due to the fact that even for $T\to 0$ the graphene layer is not strictly planar due to zero-point motion of the carbon atoms in the out-of-plane direction.

We now discuss the behavior of $A_{\|}$ as a function of $T$, with $dA_{\|}/dT<0$ at relatively low temperature. There appears a competition between two opposite factors. On one side, the actual area $A$ increases as temperature is raised, as shown in Fig. 4 for both classical and quantum simulations. On the other side, bending (rippling) of the surface causes a decrease in its 2D projection, i.e., in $A_{\|}$. At low $T$, the decrease associated to out-of-plane motion dominates the thermal expansion of the actual surface, and thus $dA_{\|}/dT<0$. This is particularly marked for the quantum results, because in this case the thermal expansion at low $T$ is very small. At high temperature, the increase in $A$ dominates the reduction in the projected area associated to out-of-plane atomic displacements. For example, for the classical results, $A$ increases roughly linearly with $T$, but the mean-square atom displacements (causing the reduction of $A_{\|}$) grow up slower than linearly, as a consequence of anharmonicity. This is in line with the discussion given by Gao and HuangGao and Huang (2014) for the thermal behavior of $A_{\|}$ observed in classical MD simulations of a single graphene layer, as well as with the trends for $dA_{\|}/dT$ calculated by Michel et al.Michel et al. (2015a)

The discussion of a real area $A$ and a projected area $A_{\|}$ presented here is reminiscent of the analysis carried out earlier in the context of biological membranes, where consideration of a “true” area is important to describe some properties of those membranes.Chacón et al. (2015) In fact, it has been shown that values of the membrane compressibility may be very different when they are related to $A$ or to $A_{\|}$. Something similar can happen for graphene, and this topic requires further research.

To summarize the results of this section, we emphasize that changes in interatomic distances and in the graphene area (both $A$ and $A_{\|}$) are important anharmonic effects. On one side, changes in the distance $d_{\rm C-C}$ and in the real area $A$ are mainly due to anharmonicity of the C-C stretching vibration. The dilation of the C-C bond displayed in Fig. 1 equally appears in a strictly planar graphene, and shows a temperature dependence similar to that of interatomic distances in solids. In real graphene, $d_{\rm C-C}$ will be also affected by out-of-plane vibrations, which in fact couple with in-plane vibrational modes, but this coupling changes $d_{\rm C-C}$ less than the only anharmonicity of the in-plane modes, at least for the conditions considered here. Something different occurs for the projected area $A_{\|}$, as manifested by the strong size effect which it displays. In this case, the coupling between in-plane and out-of-plane modes is crucial to understand the temperature dependence of $A_{\|}$ in the absence of an external stress ($P=0$). The amplitude of out-of-plane vibrations strongly depend on the size $N$,Gao and Huang (2014) causing the large size effect in $A_{\|}$. At low temperatures, all these anharmonic effects are clearly enhanced by quantum motion, as shown for example in Fig. 3 for the projected area. This is due to the fact that anharmonicity is manifested in the quantum model even in the limit $T\to 0$, contrary to the classical case, where it gradually appears as temperature is raised. This is further quantified and discussed in the next section.

In the zero-temperature limit we find in the classical approach a strictly planar graphene surface with an interatomic distance of 1.4199 Å. This corresponds to fixed atomic nuclei on the equilibrium sites without spatial delocalization, defining a minimum energy $E_{0}$, taken as a reference for our calculations. In the quantum approach, the low-temperature limit includes out-of-plane atomic fluctuations due to zero-point motion, so that the graphene layer is not strictly planar, as commented above. Moreover, anharmonicity of in-plane vibrations cause a slight zero-point bond expansion, giving an interatomic distance of 1.4287 Å, as shown in Fig. 1.

We now turn to the results of our simulations at finite temperatures, and will discuss the resulting energy for graphene. The internal energy, $E(T)$, at temperature $T$ can be written as:

where $E_{\rm el}(A)$ is the elastic energy corresponding to an area $A$, and $E_{\rm vib}(A,T)$ is the vibrational energy of the system. Our simulations directly yield $E(T)$ in both classical and quantum approaches. We can then split the internal energy $E(T)-E_{0}$ into an elastic and a vibrational part.

To define the elastic energy corresponding to an area $A$, we have carried out calculations of the (classical) energy of a supercell including 960 atoms, expanding it isotropically and keeping it flat. This strict 2D atomic layer gives us a reference to evaluate the vibrational energy. The elastic energy $E_{\rm el}(A)$ increases with $A$, and for small expansion it can be approximated as $E_{\rm el}(A)\approx K(A-A_{0})^{2}$, with $K$ = 2.41 eV/Ų. For the area $A$ derived from our classical and PIMD simulations at 300 K, we found elastic energies of 0.4 and 2.8 meV/atom, respectively. This much larger value for the quantum approach is a consequence of the larger surface $A$, as compared to the classical limit at $T$ = 300 K. Those values of the elastic energy increase to 20.8 and 22.9 meV/atom at 2000 K.

Thus, we obtain the vibrational energy $E_{\rm vib}(A,T)$ from the results of our simulations by subtracting in each case the elastic energy from the internal energy $E(T)$. In Fig. 5 we show the temperature dependence of $E_{\rm vib}(A,T)$, derived from simulations for a supercell including 960 carbon atoms. At 300 K, the vibrational energy amounts to 182 meV/atom, somewhat smaller than the values reported for diamond of 195 and 210 meV/atom, yielded by path-integral simulations with Tersoff-type and tight-binding potentials, respectively.Herrero and Ramírez (2000); Ramírez et al. (2006) This has to be compared to the value expected in a classical harmonic approximation: $E_{\rm vib}^{\rm cl}(T)=3k_{B}T$ = 77.6 meV/atom ($k_{B}$ is Boltzmann’s constant). Note also that the elastic energy $E_{\rm el}$ in the quantum approach at 300 K is 1.5% of the internal energy $E(T)-E_{0}$, most of this energy corresponding to the vibrational energy $E_{\rm vib}$.

PIMD simulations allow us to separate the potential ($E_{p}$) and kinetic ($E_{k}$) contributions to the vibrational energy,Herman et al. (1982); Giansanti and Jacucci (1988); Gillan (1990) so that $E_{\rm vib}=E_{k}+E_{p}$. In the classical limit, the kinetic energy amounts to $3k_{B}T/2$ per atom. In the quantum case, however, $E_{k}$ depends on the spatial delocalization $Q_{r}^{2}$ of the quantum paths, which can be obtained from the simulations as indicated above [see Eq. (3)]. In our simulations of graphene, both kinetic and potential energy were found to slightly increase with system size, but their convergence is rather fast. In Fig. 6, we present $E_{k}$ and $E_{p}$ derived from PIMD simulations as a function of cell size at $T$ = 300 K. Symbols indicate results of our simulations for $E_{k}$ [panel (a)] and $E_{p}$ [panel (b)], whereas dashed lines are guides to the eye. In both $E_{k}$ and $E_{p}$ there appears a shift of less than 1 meV/atom when increasing the cell size from 40 atoms to the largest size considered here (about 30,000 atoms). For cells in the order of 1000 atoms the size effect is almost inappreciable when compared to the largest cell. The potential energy is found to be smaller than the kinetic energy, indicating an appreciable anharmonicity of the lattice vibrations (see below).

Differences between the kinetic and potential contribution to the vibrational energy can be used for a quantification of the anharmonicity in condensed matter, as has been discussed earlier from path-integral simulations, e.g., for H impurities in siliconHerrero and Ramírez (1995) and van der Waals solids.Herrero (2002) In Fig. 7 we display the ratio $E_{k}/E_{p}$ between kinetic and potential parts of the vibrational energy of graphene as a function of temperature. Results are shown for classical (circles) and PIMD (squares) simulations carried out for a system size $N$ = 960 atoms. The effect of system size on the energy ratio is a slight shift in both cases, classical and quantum, as can be derived from results such as those presented in Fig. 6.

For a purely harmonic model for the vibrational modes, one expects $E_{k}=E_{p}$ (virial theoremLandau and Lifshitz (1980); Feynman (1972)), i.e., an energy ratio equal to unity at any temperature in both classical and quantum approaches. The ratio $E_{k}/E_{p}=1$ is found in the classical model for the low-temperature limit, since in this case the atomic motion does not feel the anharmonicity of the interatomic interactions, due to the vanishingly small vibrational amplitudes. This is not the case of the quantum results for $T\to 0$, because now the vibrational amplitudes remain finite and feel the anharmonicity. In particular, we find $E_{k}>E_{p}$, and at low temperature this ratio amounts to about 1.05. In the classical model, $E_{k}/E_{p}$ increases at low $T$, due to an enhancement of vibrational amplitudes, as in this approach the mean-square atomic displacements increase at low $T$ linearly with rising temperature. This increase in the energy ratio continues until about $T$ = 1000 K, and at higher temperatures $E_{k}/E_{p}$ is found to decrease, close to the quantum result. This slow decrease for rising $T$ obtained for the energy ratio from both classical and PIMD simulations does not mean that the difference $E_{k}-E_{p}$ decreases (in fact it rises with $T$), but is due to the large increase in the denominator $E_{p}$. Also, for $T\gtrsim$ 1000 K the ratio $E_{k}/E_{p}$ turns out to be larger for the classical case than in the quantum approach, even though the difference $E_{k}-E_{p}$ is similar in both cases, because in the classical approach $E_{p}$ is smaller.

Turning to the energy results for the quantum approach at $T\to 0$, it is interesting to note that earlier analysis of anharmonicity in solids, based on quasiharmonic approximations and perturbation theory indicate that low-temperature changes in the vibrational energy with respect to a harmonic calculation are mainly due to the kinetic energy. In fact, assuming perturbed harmonic oscillators (with perturbations of $x^{3}$ or $x^{4}$ type) at $T=0$, the first-order change in the energy is caused by changes in $E_{k}$, the potential energy remaining unaltered with respect to its unperturbed value.Landau and Lifshitz (1965); Herrero and Ramírez (1995)

In this section we present results for the mean-square displacements of carbon atoms in graphene, as obtained from our PIMD simulations. We mainly focus on the character of these atomic displacements, i.e., if they can be described by a classical model, or the C atoms mainly behave as quantum particles. One obviously expects that the quantum description will become more relevant as the temperature is lowered, but we show that the system size also plays an important role in this question. We use the notation introduced in Sec. II.

In Fig. 8 we display results for the motion in the out-of-plane direction $z$, derived for a cell including 96 atoms. Diamonds represent the mean-square displacement $(\Delta z)^{2}$ as a function of temperature [see Eq. (2)]. As explained above, this displacement can be split into two parts: one of them properly quantum in nature, corresponding to the spread of the quantum paths, $Q_{z}^{2}$ (squares), and another of semiclassical character, taking into account the centroid motion, i.e., global displacements of the paths, $C_{z}^{2}$ (circles). In the limit $T\to 0$, $C_{z}^{2}$ vanishes and $Q_{z}^{2}$ converges to a value of about $6\times 10^{-3}$ Ų. $Q_{z}^{2}$ decreases for increasing temperature, whereas $C_{z}^{2}$ increases almost linearly with $T$. For the system size shown in Fig. 8 (N = 96), both contributions to $(\Delta z)^{2}$ are nearly equal at $T$ = 50 K. At higher temperatures, the semiclassical part $C_{z}^{2}$ is the main contribution to the atomic displacements on the out-of-plane direction. Values of $C_{z}^{2}$ obtained from our PIMD simulations coincide within error bars with the mean-square atomic displacement in the $z$ direction obtained from classical MD simulations. The proper quantum delocalization can be estimated from the mean extension of the quantum paths in the $z$ direction, i.e., from $Q_{z}^{2}$. At 25 and 300 K we find an average extension $(\Delta z)_{Q}=(Q_{z}^{2})^{1/2}\approx$ 0.06 and 0.03 Å, respectively.

The vibrational amplitudes in the layer plane $(x,y)$ are smaller than in the $z$ direction. This is true for both, quantum and classical contributions. In the limit $T\to 0$, zero-point motion yields a mean-square displacement on the plane $Q_{\|}^{2}=3.6\times 10^{-3}$ Ų, which means $Q_{x}^{2}$ = $Q_{y}^{2}$ = $1.8\times 10^{-3}$ Ų. Thus, for each direction in the $(x,y)$ plane we find a value about three times smaller than the zero-point mean-square displacement in the $z$ direction. This is due, of course, to the low-frequency out-of-plane modes, which cause a larger vibrational amplitude. At finite temperatures, the difference between vibrational amplitudes in-plane and out-of-plane are larger than at $T=0$. At 300 K, we find $(\Delta z)^{2}/(\Delta x)^{2}$ = 11.9 and 61.8 for simulation cells with 96 and 960 atoms, respectively. This is a consequence of the appearance of long-wavelength modes (with low frequency) with larger vibrational amplitudes in the larger cells, as explained below.

The picture shown in Fig. 8 for the atom displacements in the $z$ direction is qualitatively the same for different system sizes, but the temperature region where $Q_{z}^{2}$ or $C_{z}^{2}$ is the main contribution to $(\Delta z)^{2}$ is strongly dependent on the system size $N$. This is mainly due to the enhancement of the classical-like contribution for increasing size, a fact already observed and described earlier from results of classical MD simulations.Los et al. (2009); Gao and Huang (2014); Ramírez et al. (2016) In Fig. 9 we present mean-square displacements $Q_{z}^{2}$ and $C_{z}^{2}$ in the out-of-plane direction at $T$ = 25 K for different $N$. One observes that the quantum contribution $Q_{z}^{2}$ converges fast as $N$ is increased, and in fact for $N>50$, changes in $Q_{z}^{2}$ are much smaller than the symbol size. Clear finite-size effects are only observed for very small simulation cells. The semiclassical contribution $C_{z}^{2}$ is small for small $N$, and grows with system size, so that at $T$ = 25 K it becomes similar to $Q_{z}^{2}$ for $N\sim 200$. Calling $N_{c}$ the system size at which $C_{z}^{2}$ and $Q_{z}^{2}$ are equal, we found that $N_{c}$ decreases for rising temperature, and for 50 K we have $N_{c}\sim 100$ [see Fig. 8].

On the other side, for a given system size $N$, the ratio $Q_{z}^{2}/C_{z}^{2}$ decreases for increasing $T$, and there is a crossover temperature $T_{c}$ for which this ratio is unity. For $T>T_{c}$ classical-like motion dominates the atomic mean-square displacement in the $z$ direction. In Fig. 10 we show $T_{c}$ as a function of the system size $N$. Symbols are results derived from our PIMD simulations. The dashed line through the data points was obtained as a linear fit of $\log T$ vs $\log N$ for the crossover temperatures derived from the simulations for several cell sizes. This is consistent with a power-law dependence of the crossover temperature on system size, i..e, $T_{c}=a\,N^{-b}$ with a factor $a=1002$ K and an exponent $b=0.67$.

This means that for a given temperature (even very low), there is a system size for which classical-like motion along the out-of-plane direction dominates over atomic quantum delocalization. Thus, extrapolating the dashed line in Fig. 10, we find for $T$ = 1 K and 0.1 K crossover sizes $N_{c}\approx 3\times 10^{4}$ and $9\times 10^{5}$, respectively. These values are much larger than the simulation cells which can be dealt with in PIMD simulations at low temperatures. We note, however, that the dependence of $T_{c}$ on $N$ may depart from the power-law given above (dashed line in Fig. 10) for larger system sizes. This is related to the size dependence of $C_{z}^{2}$, which has been proposed to increase logarithmically with $N$ (from an atomistic model for the spectral amplitudesGao and Huang (2014); Ramírez et al. (2016)) or as a power-law (based on a scaling theory in the continuum medium approachLos et al. (2009)). In any case, this difference would be detected for system sizes much larger than those considered here in the PIMD simulations.

This can be rationalized in terms of vibrational modes in the $z$ direction appearing for different system sizes. The mean-square displacements $C_{z}^{2}$ and $Q_{z}^{2}$ can be estimated in a harmonic approximation in the following way. The classical part is given by:

where $m$ is the atomic mass and the sum is extended to the wavevectors ${\bf k}=(k_{x},k_{y})$ in the reciprocal lattice corresponding to the two-dimensional simulation cell.Ramírez et al. (2016) This expression is valid for relatively low temperatures, as for high $T$ anharmonic effects cause $C_{z}^{2}$ to rise sublinearly with $T$.Gao and Huang (2014) On the other side, in the limit $T\to 0$, $Q_{z}^{2}$ converges to

Increasing the system size $N$ causes the appearance of vibrational modes with longer wavelength $\lambda$. In fact, one has an effective cut-off $\lambda_{max}\approx l$, with $l=(NA_{\|})^{1/2}$. Calling $k=|{\bf k}|$, this translates into $k_{min}=2\pi/\lambda_{max}$, i.e., $k_{min}\sim N^{-1/2}$. We consider a dispersion relation $\omega({\bf k})$ of the form $\rho\,\omega^{2}=\sigma k^{2}+\kappa k^{4}$, consistent with an atomic description of grapheneRamírez et al. (2016) ($\rho$, surface mass density; $\sigma$, effective stress; $\kappa$, bending modulus). Then, the sum in Eq. (7) diverges in the large-size limit due to the $\omega({\bf k})^{2}$ term in the denominator and the two-dimensional character of the wavevector ${\bf k}$Gao and Huang (2014); Ramírez et al. (2016) (see results for $C_{z}^{2}$ vs $N$ at 25 K in Fig. 9). However, the sum in Eq. (8) converges to a finite value for large $N$, in agreement with the results of our PIMD simulations at low temperature (see results for $Q_{z}^{2}$ vs $N$ in Fig. 9).

Given a temperature $T$, modes with frequency $\omega\lesssim\omega_{c}(T)=k_{B}T/\hbar$ may be considered in the classical regime. Then, increasing $N$ one reaches a system size $N_{c}(T)$ for which any further increase introduces new modes with frequency $\omega<\omega_{c}$, and therefore contributing to increase the classical displacement $C_{z}^{2}$ more than the quantum one $Q_{z}^{2}$. Thus, at any finite temperature classical-like displacements will dominate over quantum delocalization in the out-of-plane direction, provided that the system size is larger than the corresponding $N_{c}$.

All this can be visualized in terms of the harmonic model. Taking into account that at low temperatures $\sigma\approx 0$, the dispersion relation in the $z$ direction is $\omega({\bf k})^{2}\approx\kappa k^{4}/\rho$ (see Ref. Ramírez et al., 2016). At temperature $T$, the mean-square displacement can be approximated as

Then, from Eqs. (7) and (9), and remembering that $Q_{z}^{2}=(\Delta z)^{2}-C_{z}^{2}$, we calculate for each $N$ the temperature at which $Q_{z}^{2}=C_{z}^{2}$. This yields the curve shown as a dashed-dotted line in Fig. 10, displaying a dependence of $T_{c}$ on $N$ very similar to that derived from PIMD simulations for system sizes accessible in our calculations. In fact, for $N\lesssim 1000$, the results can be approximated in both cases by a power-law with nearly the same exponent. For larger sizes, the harmonic calculation deviates to temperatures somewhat smaller than the extrapolation of the simulation data. At this point, we cannot assure if this deviation is a real trend of the physical system or only a consequence of the harmonic approach. The main limitations of this approximation are the neglect of anharmonicity in the vibrational modes (expected to be reasonably small at low $T$) and the overestimation of frequencies in the ${\bf k}$-space region far from the $\Gamma$ point (k = 0). In that region, $\omega$ increases with $k$ slower than $k^{2}$ (see Ref. Mounet and Marzari, 2005; Michel et al., 2015b; Ramírez et al., 2016). Even with these limitations the harmonic model captures qualitatively, and almost quantitatively, the basic aspects of the competition between classical-like and proper quantum dynamics of the carbon atoms in the $z$ direction.

We finally note that the competition between $C_{z}^{2}$ and $Q_{z}^{2}$ as a function of $N$ discussed here appears only in the $z$ direction. For motion in the $(x,y)$ plane, both $C_{x}^{2}$ and $Q_{x}^{2}$ converge fast with increasing system size to their asymptotic value, similarly to the behavior known for vibrational motion in 3D solids.Herrero and Ramírez (2014)