# Path-integral simulation of graphene monolayers under tensile stress

###### Abstract

Finite-temperature properties of graphene monolayers under tensile stress have been studied by path-integral molecular dynamics (PIMD) simulations. This method allows one to consider the quantization of vibrational modes in these crystalline membranes and to analyze the influence of anharmonic effects in the membrane properties. Quantum nuclear effects turn out to be appreciable in structural and thermodynamic properties of graphene at low temperature, and they can even be noticeable at room temperature. Such quantum effects become more relevant as the applied stress is increased, mainly for properties related to out-of-plane atomic vibrations. The relevance of quantum dynamics in the out-of-plane motion depends on the system size, and is enhanced by tensile stress. For applied tensile stresses, we analyze the contribution of the elastic energy to the internal energy of graphene. Results of PIMD simulations are compared with calculations based on a harmonic approximation for the vibrational modes of the graphene lattice. This approximation describes rather well the structural properties of graphene, provided that the frequencies of ZA (flexural) acoustic modes in the transverse direction include a pressure-dependent correction.

###### pacs:

61.48.Gh, 65.80.Ck, 63.22.Rc## I Introduction

In recent years there has been a surge of interest on two-dimensional materials, and graphene in particular, due to their unusual electronic, elastic and thermal properties.Geim and Novoselov (2007); Castro Neto et al. (2009); Flynn (2011); Roldan et al. (2017) In fact, graphene displays high values of thermal conductivity,Ghosh et al. (2008); Nika et al. (2009); Balandin (2011) as well as large in-plane elastic constants.Lee et al. (2008) Its mechanical properties are also important for possible applications, such as cooling of electronic devices.Seol et al. (2010); Prasher (2010)

The structural arrangement for pure defect-free graphene corresponds to a planar honeycomb lattice. At finite temperatures, there appear out-of-plane displacements of the C atoms, and for , quantum fluctuations related to zero-point motion give rise to a departure of strict planarity of the graphene sheet.Herrero and Ramírez (2016) In particular, one has low-lying vibrational excitations associated to large-scale ripples perpendicular to the plane.Fasolino et al. (2007) Moreover, a graphene sheet can actually bend and depart from planarity for other reasons, such as the presence of defects and external stresses.Meyer et al. (2007); de Andres et al. (2012a)

A thin membrane crumples in the presence of a compressive stress. This behavior has been investigated during the last three decades in lipid membranesEvans and Rawicz (1990); Safran (1994) and polymer films.Cerda and Mahadevan (2003); Wong and Pellegrino (2006) In graphene, crumpling originates from out-of-plane phonons as well as from static wrinkling, and has been observed in both supported and freestanding samples.Kirilenko et al. (2011); Nicholl et al. (2015) Mechanical properties such as stiffness and bending rigidity can be renormalized due to crumpling.Ruiz-Vargas et al. (2011); Kosmrlj and Nelson (2013, 2014) For graphene, it has been found that the maximum compressive stress that a freestanding sheet can sustain without crumpling decreases with system size, and has been estimated to be about 0.1 N/m at room temperature in the thermodynamic limit.Ramírez and Herrero (2017)

A tensile stress applied in the graphene plane does not affect the planarity of the sheet, but causes appreciable changes in the elastic properties of the material.López-Polín et al. (2017) Thus, it has been observed that the in-plane Young modulus is increased by a factor of three when applying a stress of 1 N/m.Ramírez and Herrero (2017) The bending rigidity does also change with the tensile stress, and in fact it decreases but not so critically as the in-plane elastic constants. In this context, it is important to note that the actual area per atom, , is not readily measurable, and the accessible observable is usually its projection, , onto the mean plane of the membrane (). Thus, one may refer the elastic properties of graphene either to the area or to , which may behave in very different ways. For example, a negative thermal expansion coefficient is found for graphene when one refers to , but the thermal expansion associated to the area is positive.Pozzo et al. (2011); Herrero and Ramírez (2016)

Recent experimental and theoretical work has shown the influence of strain in several characteristics of graphene, such as electronic transport, optical properties, and the formation of moiré patterns.Woods et al. (2014) Similar properties have been also studied in other two-dimensional materials, as metallic dichalcogenides.Amorim et al. (2016)

Equilibrium and dynamical properties of graphene have been studied earlier by using Monte Carlo and molecular dynamics simulations. These simulations were based on ab-initio,Shimojo et al. (2008); Pozzo et al. (2011); de Andres et al. (2012b); Chechin et al. (2014) tight-binding,Herrero and Ramírez (2009); Cadelano et al. (2009); Akatyeva and Dumitrica (2012); Lee et al. (2013) and empirical interatomic potentials.Fasolino et al. (2007); Shen et al. (2013); Los et al. (2009); Ramírez et al. (2016); Magnin et al. (2014); Los et al. (2016) In most of these simulations, C atoms were treated as classical particles, which is accurate at relatively high temperatures but is not suitable to study thermodynamic variables at low temperature. The quantum character of the atomic motion can be taken into account by employing path-integral simulations, which allow to consider quantum and thermal fluctuations in many-body systems at finite temperatures.Gillan (1988); Ceperley (1995) Path-integral simulations of a single graphene sheet have been lately performed to study equilibrium properties of this material.Brito et al. (2015); Herrero and Ramírez (2016) Moreover, nuclear quantum effects have been studied by a combination of density-functional theory and a quasi-harmonic approximation for vibrational modes in this crystalline membrane.Mounet and Marzari (2005); Shao et al. (2012)

In this paper, we employ path-integral molecular dynamics (PIMD) simulations to study structural and vibrational properties of graphene under tensile stress. We consider different sizes for the simulation cell, as finite-size effects are known to be important for some properties of graphene.Gao and Huang (2014); Los et al. (2016); Herrero and Ramírez (2016) The magnitude of nuclear quantum effects in the graphene properties is assessed by comparing the results of PIMD simulations with data obtained from classical simulations. We find that quantum effects are relevant to describe the temperature and pressure dependence of graphene’s real and in-plane areas, as well as to describe the amplitude of the out-of-plane motion, especially at low temperatures. Our data indicate that the relevance of nuclear quantum effects increases as tensile stress is raised. Results of PIMD simulations are compared with data derived from a harmonic approximation for the out-of-plane vibrations. This approximation turns out to be rather accurate, provided that the vibrational frequencies of ZA acoustic modes are conveniently renormalized for different applied stresses.

The paper is organized as follows. In Sec. II, we present the computational method used in the simulations. Structural properties such as in-plane and real area are given in Sec. III as a function of applied stress. Results for the internal, vibrational, and elastic energy of graphene are discussed in Sec. IV. In Sec. V we study the out-of-plane atomic motion, with emphasis on the competition between classical-like and quantum dynamics. In Sec. VI we summarize the main results.

## Ii Computational Method

We use the PIMD method to obtain equilibrium properties of graphene under tensile stress. This procedure is based on the Feynman path-integral formulation of statistical mechanics, a nonperturbative technique to study many-body quantum systems at finite temperatures.Feynman (1972) The implementation of this formulation for numerical simulations is based on an isomorphism between the quantum system and a fictitious classical system, in which each quantum particle is described by a ring polymer (corresponding to a cyclic quantum path) composed of (Trotter number) beads.Kleinert (1990) This becomes exact in the limit . Details on this simulation technique can be found elsewhere.Chandler and Wolynes (1981); Gillan (1988); Ceperley (1995); Herrero and Ramírez (2014) The dynamics in PIMD is artificial, since it does not correspond to the actual dynamics of the real quantum particles. However, it is useful for sampling the many-body configuration space, yielding accurate results for time-independent equilibrium properties of the actual quantum system.

The Born-Oppenheimer surface for the nuclear dynamics is derived here from an effective empirical potential, developed for carbon-based systems, namely the so-called LCBOPII.Los et al. (2005) This is a long-range carbon bond order potential, which was previously used to perform classical simulations of diamond,Los et al. (2005) graphite,Los et al. (2005) liquid carbon,Ghiringhelli et al. (2005a) as well as graphene sheets.Zakharchenko et al. (2009); Fasolino et al. (2007); Los et al. (2016) A relevant application of this effective potential was the calculation of the carbon phase diagram including diamond, graphite, and the liquid, and showing its precision by comparison of the predicted diamond-graphite line with experimental results.Ghiringhelli et al. (2005b)

The LCBOPII potential has been more recently employed to study graphene, giving a good description of elastic properties such as the Young’s modulus.Zakharchenko et al. (2009); Politano et al. (2012) According to previous simulations,Ramírez et al. (2016); Herrero and Ramírez (2016); Ramírez and Herrero (2017) the original LCBOPII parameterization has been slightly modified to increase the zero-temperature bending constant of graphene from 1.1 eV to a value of 1.49 eV, more consistent with experimental data.Lambin (2014) This effective potential was lately used to perform PIMD simulations, allowing to assess the extent of quantum effects in graphene sheets from a comparison with results of classical simulations.Herrero and Ramírez (2016)

Other effective interatomic potentials have been employed in recent years to study various properties of graphene. In particular, the AIREBO potential model has been widely used in this field.Sfyris et al. (2015); Memarian et al. (2015); Zou et al. (2016); Anastasi et al. (2016); Ghasemi and Rajabpour (2017) Comparing the LCBOPII and AIREBO models, we find that they yield very similar equilibrium C–C distance and in-plane thermal expansion coefficient, as derived from classical molecular dynamics simulations.Memarian et al. (2015); Ghasemi and Rajabpour (2017); Herrero and Ramírez (2016) Results for the Young’s modulus of graphene derived from the LCBOPII potential are closer to those given by ab initio calculations.Memarian et al. (2015)

The calculations presented here were carried out in the isothermal-isobaric ensemble, where we fix the number of carbon atoms (), the applied stress (), and the temperature (). We employed effective algorithms for carrying out PIMD simulations in this statistical ensemble, as those presented in the literature.Tuckerman et al. (1992); Tuckerman and Hughes (1998); Martyna et al. (1999); Tuckerman (2002) Specifically, we used staging variables to define the bead coordinates, and the constant-temperature ensemble was achieved by coupling chains of four Nosé-Hoover thermostats. A supplementary chain of four barostats was coupled to the area of the simulation box to give the required pressure .Tuckerman and Hughes (1998); Herrero and Ramírez (2014) The kinetic energy has been calculated by using the so-called virial estimator, which has a statistical uncertainty smaller than the potential energy, , of the system.Herman et al. (1982); Tuckerman and Hughes (1998) Other technical details about the simulations presented here are the same as those given elsewhere.Herrero et al. (2006); Herrero and Ramírez (2011); Ramírez et al. (2012)

Atomic forces were analytically derived from position derivatives of the instantaneous potential energy (note that ). The estimator of the two-dimensional (2D) stress tensor is the same as that employed in previous works,(Ramírez et al., 2016; Ramírez and Herrero, 2017) and its formulation for PIMD simulations of graphene is similar to that given earlier for three-dimensional (3D) solids.Ramírez et al. (2008); Herrero and Ramírez (2014) Details on the pressure estimator employed here are presented in Appendix A. The mechanical stress in the plane of graphene is obtained from the trace of the tensor :

(1) |

Note that in the case of applying a large compressive stress (, not considered here), one may have severe bending or crumpling of the graphene sheet. In this case the in-plane stress may appreciably differ from the actual stress felt by the real graphene surface (related to the area ).Ramírez and Herrero (2017); Fournier and Barbetta (2008)

We consider rectangular simulation cells with similar side lengths and in the and directions of the reference plane, and periodic boundary conditions were assumed. Sampling of the configuration space has been carried out at temperatures between 12 K and 2000 K. For comparison with results of PIMD simulations, some classical molecular dynamics (MD) simulations have been also performed. In our context this is achieved by by setting = 1. For the quantum simulations, was taken proportional to the inverse temperature: = 6000 K, which roughly gives a constant precision in the PIMD results at different temperatures.Herrero et al. (2006); Herrero and Ramírez (2011); Ramírez et al. (2012) Cells of size up to 8400 and 14720 atoms were considered for PIMD and classical MD simulations, respectively. For a given temperature, a typical simulation run consisted of PIMD steps for system equilibration, followed by steps for the calculation of ensemble average properties.

## Iii Structural properties

The simulations presented here were performed in the isothermal-isobaric ensemble, as explained above in Sec. II. Thus, in a simulation run we fix the number of carbon atoms , the temperature , and the applied stress in the plane, allowing for changes in the in-plane area of the simulation cell for which periodic boundary conditions are assumed. Carbon atoms are free to move in the out-of-plane direction ( coordinate), and in general any measure of the real surface of a graphene sheet at should give a value larger than the area of the simulation cell in the plane. In this line, there has appeared in recent years a discussion in the context of biological membranes, dealing with the question whether it is more convenient to describe the properties of those membranes using the concept of a real surface rather than a projected (in-plane) surface.Imparato (2006); Waheed and Edholm (2009); Chacón et al. (2015) The same question has been also recently raised for crystalline membranes such as graphene.Pozzo et al. (2011); Herrero and Ramírez (2016); Nicholl et al. (2017); Ramírez and Herrero (2017) This can be important for addressing the calculation of thermodynamic properties, because the in-plane area is the variable conjugate to the stress used in our simulations, and the real area (also called effective, true, or actual area in the literatureImparato (2006); Fournier and Barbetta (2008); Waheed and Edholm (2009); Chacón et al. (2015)) is conjugate to the usually-called surface tension.Safran (1994) The difference has been recently denoted as hidden area by Nicholl et al.Nicholl et al. (2017)

The real area in 3D space is calculated here by a triangulation based on the atomic positions along a simulation run. is obtained from the areas associated to structural hexagons. Each hexagon contributes by a sum of six triangles, each one formed by the positions of two neighboring carbon atoms and the barycenter of the hexagon (mean point of its six vertices).Ramírez and Herrero (2017) There are other similar definitions that can be employed for the area , as those based on the interatomic distance C–C.Hahn et al. (2016); Herrero and Ramírez (2016) The area based on triangulation employed here seems more precise to deal with the 2D nature of a graphene layer in 3D space. It has been shown earlier that has a very small size effect, in fact negligible in comparison with that appearing for .Ramírez and Herrero (2017)

In Fig. 1 we show the temperature dependence of the in-plane area , obtained from classical MD (open circles) and PIMD simulations (solid squares) for a supercell with = 960 atoms. Results are given for = 0, , and eV Å. Tensile stress causes not only an increase in , but its temperature dependence also changes. For each considered value of the stress, the curve derived from quantum simulations displays a minimum, that shifts to lower temperatures as the tensile stress is increased. Thus, such a minimum evolves from 1000 K for to 400 K for eV Å. In the classical simulations, however, one finds a shallow minimum for , that is absent for the tensile stresses shown in Fig. 1 (in fact we did not observe it for eV Å either, not shown in the figure). The classical results for 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)

At low the results of PIMD simulations verify , i.e., the corresponding curves shown in Fig. 1 (solid symbols) become flat close to , as required by the third law of thermodynamics. In the limit , the difference between quantum and classical results converges to 0.022 Å/atom in the absence of applied stress (). This difference decreases for rising temperature, as nuclear quantum effects become less important. For eV Å, we find for a difference of 0.031 Å/atom. The increase in at low temperature is due to zero-point motion associated to in-plane acoustic modes (LA and TA). The frequency of these modes decreases for increasing (i.e., when tensile stress is increased, according to positive Grüneisen parameters), and therefore their vibrational amplitudes are larger. This causes a larger zero-point expansion of for larger tensile stress.

The presence of a minimum in the curves derived from PIMD simulations is due to two competing effects, as discussed earlier for graphene without stress.Gao and Huang (2014); Michel et al. (2015); Herrero and Ramírez (2016) On one hand, the area increases as temperature is raised, and on the other hand, surface bending gives rise to a decrease in its 2D projection, i.e., . At low , this decrease associated to out-of-plane motion dominates the thermal expansion of the real surface, and . For the quantum results, the thermal expansion at low is very small compared to the classical calculations for which , thus causing a more appreciable decrease in for raising in the quantum case. At high temperatures, the increase in predominates over the contraction in the projected area due to out-of-plane motion.

For unstressed graphene it has been indicated that finite-size effects can be important for several structural properties of the crystalline membrane.Herrero and Ramírez (2016); Ramírez and Herrero (2017) It is now worthwhile to consider finite-size effects for the in-plane area of graphene under stress. In Fig. 2 we present the size dependence of for several tensile stresses. In all cases, decreases for increasing , and reaches a well-defined plateau for large sizes. One observes that the convergence to the large-size value is faster for larger tensile stress. Moreover, the difference between the large-size limit and the value corresponding to = 24 (the smallest supercell considered here) appreciably decreases from Å/atom for to Å/atom for eV Å.

For comparison, we also present in Fig. 2 results for derived from classical MD simulations at 300 K. The difference between quantum and classical results for = 0 amounts to 0.017 Å/atom, and it is nearly constant for the system sizes considered here. This difference increases to 0.022 Å/atom at in the low-temperature limit, as indicated above (see Fig. 1). For eV Å and K, our classical simulations yield an curve similar to the quantum one (not shown in Fig. 2 to avoid overcrowding). In particular, for a system size , we found an in-plane area = 2.6423 Å/atom, so that the difference between classical and quantum results at this tensile stress is similar to that found for . It is interesting to note that the increase in area due to quantum nuclear motion at 300 K is the same as that caused by a relatively large tensile stress of about eV Å ( N/m).

We now turn to the real surface of graphene and its measure through the area . It was shown earlier from simulations at that the surface is larger than , and the difference between both increases with temperature. This is clear from the fact that is a 2D projection of , and the actual surface becomes increasingly bent as temperature is raised and the amplitude of out-of-plane atomic vibrations becomes larger. An important difference between the temperature dependence of and is that the latter first decreases for increasing and then it increases at higher , with a minimum at a temperature . For rising tensile stress, the vibrational amplitude in the direction decreases (see below), so that the temperature of minimum is lowered. This becomes even clearer in the results of classical simulations, for which the shallow minimum in the curve disappears at relatively low pressures, and it is not observed in the data presented for = and eV Å in Fig. 1. For the area one does not observe the decrease displayed by in both classical and quantum simulations at low temperatures (see Ref. Herrero and Ramírez, 2016 for results at ).

In Fig. 3 we present the areas and vs. tensile stress for a simulation cell including 960 atoms. In both cases, we present results from classical (open symbols) and PIMD (solid symbols) simulations. Circles correspond to the in-plane area , whereas squares represent data for the real area . One notices that quantum effects are appreciable at room temperature. The main aspects of this figure are the following. Tensile stress causes an increase of about 5% in both and from to eV Å. Moreover, quantum nuclear effects cause in both cases a surface expansion of about 0.02 Å/atom, which increases slightly as the tensile stress is raised.

To make connection of our results derived from atomistic simulations with an analytical formulation of crystalline membranes, we note that the relation between and can be expressed in the continuum limit (macroscopic view) asImparato (2006); Waheed and Edholm (2009); Ramírez and Herrero (2017)

(2) |

where is the height of the membrane surface, i.e. the distance to the reference plane. The difference can be calculated in a classical approach by Fourier transformation of the r.h.s. of Eq. (2).Safran (1994); Chacón et al. (2015); Ramírez and Herrero (2017) This requires the introduction of a dispersion relation for out-of-plane modes (ZA band), where are 2D wavevectors. The frequency dispersion in this acoustic (flexural) band can be well approximated by the expression , consistent with an atomic description of grapheneRamírez et al. (2016) (; , surface mass density; , effective stress; , bending modulus). The effective stress can be written as , with a term that appears at finite temperature even in the absence of an applied stress () due to out-of-plane motion (at 300 K, eV Å).Ramírez et al. (2016)

After Fourier transformation one has for the area per atom:Safran (1994); Chacón et al. (2015)

(3) |

For large the sum in Eq. (3) can be approximated by an integral:Ramírez and Herrero (2017)

(4) |

The limits in the integral are the cut-off and the size-dependent minimum wavevector , with . The integral in Eq. (4) converges provided that , which is the case here. It allows us to explicitly write the size-dependent ratio as

(5) |

with the large-size limit ( or ):

(6) |

Eq. (5), although in principle not very accurate for small system size, yields for , , and = 300 K ( eV Å, = 1.7 eV; see Ref. Ramírez and Herrero, 2017) a shift in of , which translates into an increase in of Å/atom with respect to the large-size limit. From the results of our simulations we find a size effect in of Å/atom for . Note that, apart from the replacement of the sum in Eq. (3) by an integral, the above expressions assume harmonic vibrations for out-of-plane motion, which becomes less accurate as temperature increases for the onset of larger anharmonicity. Note also the appearance of the stress () in the logarithmic term in Eq. (5), which causes that an increase in tensile stress ( more negative) gives rise to a faster convergence of the area with system size, according to the results shown in Fig. 2.

## Iv Internal energy

At and zero applied stress we find with the LCBOPII potential in a classical approach a strictly planar graphene surface with an interatomic distance = 1.4199 Å, i.e., an area of 2.6189 Å per atom, which we call . This corresponds to a graphene sheet with fixed atomic nuclei on their equilibrium sites without spatial delocalization, giving the minimum energy , taken as a reference for our calculations at nonzero temperature and applied stress. In a quantum approach, the limit includes out-of-plane atomic fluctuations associated to zero-point motion, and the graphene sheet is not strictly planar. In addition, anharmonicity of in-plane vibrations gives rise to a zero-point lattice expansion (increase in the area , see Sec. III), which for yields an interatomic distance = 1.4287 Å, around 1% larger than the classical minimum.

The internal energy is calculated as a sum of the kinetic and potential energy obtained from the simulations at a given temperature. In our simulations of graphene, and have been found to slightly increase with system size, and their convergence is rather fast. Thus, for cells in the order of 200 atoms the size effect in the internal energy is almost inappreciable when compared to the largest cells.Herrero and Ramírez (2016) The kinetic energy is associated to vibrational motion of carbon atoms (in-plane and out-of-plane), but the potential energy includes contributions due to atomic vibrations and to the elastic energy due to changes in the area of graphene at finite temperatures and applied stresses (see below).

In Fig. 4 we display the stress dependence of the internal energy, , as derived from classical and PIMD simulations for system size = 960. Results are shown for = 25 K (circles) and 300 K (squares). Open and solid symbols correspond to classical and PIMD simulations, respectively. At , the classical energy per atom is basically given by the vibrational energy , as follows from the equipartition theorem in a harmonic approximation (HA). As the tensile stress is increased ( more negative), the classical internal energy increases for a given temperature, due to the contribution of the elastic energy associated to a finite strain in the graphene lattice. The behavior of the quantum results shown in Fig. 4 is similar to the classical ones. The main difference is the increase in internal energy caused by quantization of the nuclear motion. For zero stress, this increase amounts to 165 meV/atom at 25 K and 107 meV/atom at 300 K. These shifts do not appreciably change in the stress region considered here, and in fact the difference between quantum and classical results appears to be nearly constant in the results shown in Fig. 4.

To analyze the different contributions to the internal energy, , we write

(7) |

In this expression is the elastic energy corresponding to an area , and is the vibrational energy of the system. Although not explicitly indicated, the area is a function of the stress and temperature . Our simulations give , and using Eq. (7) we can then split the internal energy into an elastic and a vibrational part. The vibrational contribution can in turn be split into kinetic and potential energy parts: .

The elastic energy is defined here as the increase in energy corresponding to a strictly planar graphene layer with area respect to the minimum energy . Thus, we have calculated for a supercell including 960 carbon atoms, expanding it isotropically and keeping it flat. Finite-size effects on the elastic energy are very small, and in fact negligible for our current purposes, as happens for the size dependence of the area . For , the elastic energy increases with , and for small lattice expansion it can be approximated as , with = 2.41 eV Å. At room temperature ( 300 K) and for small stresses ( close to ), the elastic energy is much smaller than the vibrational energy , but this can be different for low and/or large applied stresses (see below). Once calculated the elastic energy for the area resulting from the simulations at given and , we obtain the vibrational energy by subtracting the elastic energy from the internal energy : (see Eq. (7)). Then, the potential energy corresponding to vibrational motion, , is found as .

In Fig. 5 we present the different contributions to the internal energy of graphene vs. tensile stress, as derived from PIMD simulations for = 960 atoms and = 25 K. In this figure, diamonds represent the elastic energy, and circles and squares indicate and , respectively. At , is close to zero, but slightly positive, as a consequence of the zero-point expansion of the graphene lattice, which causes that . For increasing tensile stress, rises and becomes similar to and for eV Å. At still larger tensile stress, the elastic contribution is the largest one, as shown in Fig. 5. For = 300 K the picture is qualitatively the same. The elastic energy increases roughly a constant value (26 meV/atom) with respect to the results at 25 K in the whole stress range shown in Fig. 5. The same happens for the kinetic energy, with a rise of 6 meV/atom. As a result, the crossing of and at 300 K occurs for a tensile stress eV Å.

For a purely harmonic model for the vibrational modes, one expects (virial theoremLandau and Lifshitz (1980); Feynman (1972)), i.e., an energy ratio at any temperature in both classical and quantum approaches. This is not strictly the case for the results of our PIMD, because the vibrational amplitudes are finite, even at low temperatures, and feel the anharmonicity of the interatomic potential. In particular, we find , for all temperatures and tensile stresses considered here. As displayed in Fig. 5 for = 25 K, the difference increases as the tensile stress is raised, so that is about 5% larger than for small stress, and around 9% larger for a stress of eV Å. Differences between the kinetic and potential contribution to the vibrational energy have been used for a quantification of the anharmonicity in condensed matter, as discussed earlier from path-integral simulations, e.g., for van der Waals solidsHerrero (2002) and H impurities in silicon.Herrero and Ramírez (1995)

Concerning the energy results for our quantum approach at low , we note that analyses 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 mostly due to the kinetic energy. Thus, considering perturbed harmonic oscillators with perturbations of type or at , first-order changes in the energy are due to changes in , and the potential energy stays unshifted in its unperturbed value.Landau and Lifshitz (1965); Herrero and Ramírez (1995)

## V Out-of-plane motion

In this section we study the mean-square displacements of carbon atoms in the direction, normal to the graphene sheet, as obtained from our PIMD simulations. We mostly concentrate on the nature of these atomic displacements, i.e., if they can be well described by a classical model, or the C atoms appreciably behave as quantum particles. We expect of course that a classical description will lose accuracy as the temperature is reduced, but in the case of graphene it has been shown earlier that other factors such as the system size play also an important role in this question.Herrero and Ramírez (2016) Moreover, an external stress modifies the vibrational frequencies in the material, thus causing a change in the vibrational amplitudes and in the accuracy of a classical description at a given temperature.

PIMD simulations can be used to study vibrational amplitudes or atomic delocalization at finite temperatures. This includes a thermal (classical-like) motion, as well as a delocalization due to the quantum nature of the atomic nuclei, which can be quantified by the spacial extension of the paths associated to a given atomic nucleus. For each quantum path, we define the center-of-gravity (centroid) as

(8) |

where is the 3D position of bead in the ring polymer associated to nucleus . For the out-of-plane motion, we focus on the -coordinate of the polymer beads. Then, the mean-square displacement of the atomic nucleus in the direction along a PIMD simulation run is defined as

(9) |

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)

(10) |

The total spatial delocalization of atomic nucleus in the direction at a finite temperature includes, in addition to , another contribution which accounts for the classical-like motion of the centroid coordinate , i.e.

(11) |

with

(12) |

This term converges at high to the mean-square displacement given by a classical model, since in this limit each quantum path collapses onto a single point (). In the results presented below, we will show data for , , and , calculated as averages for the atoms in a simulation cell. For example:

(13) |

To connect the results of our simulations with the out-of-plane displacements corresponding to vibrational modes of the graphene sheet, we recall that the atomic mean-square displacement at temperature is given in a HA by

(14) |

where the index ( = 1, 2) refers to the phonon bands ZA and ZO, with atomic displacements along the direction.Mounet and Marzari (2005); Karssemeijer and Fasolino (2011); Wirtz and Rubio (2004) The sum in is extended to wavevectors in the hexagonal Brillouin zone, with discrete points spaced by and .Ramírez et al. (2016) Eq. (14) has been used here to calculate in a harmonic approach. Increasing the system size causes the appearance of vibrational modes with longer wavelength . In fact, one has for the phonons an effective wavelength cut-off , with , and the minimum wavevector is , i.e., . For the calculations presented below, based on the formula in Eq. (14), we have used the vibrational frequencies in the phonon branches ZA and ZO obtained from diagonalization of the dynamical matrix corresponding to the LCBOPII potential employed here.

For an applied stress , the most important effect in the ZA and ZO bands is a change in frequency of the ZA modes in the low-frequency region, for which

(15) |

The zero-stress band calculated for the minimum-energy structure (area ), verifies for small : . Then, for the small- region is dominated by the quadratic term (linear in ) in Eq. (15), so that for 1 Å.

In Fig. 6 we show results for the motion in the out-of-plane direction, obtained for a cell including 96 atoms. The use of a relatively small simulation cell is convenient to visualize the behavior of in the low-temperature region, where quantum effects are prominent. For larger cell sizes, these effects appear only at lower temperatures, which turns out to be difficult to observe from PIMD simulations.Herrero and Ramírez (2016) Solid symbols represent the mean-square displacement derived from the simulations, as a function of temperature for three different applied stresses: (circles), (squares), and eV Å (diamonds). Lines were calculated from a harmonic approximation based on the ZA and ZO phonon bands of graphene, as indicated above.

One observes in Fig. 6 that the vibrational amplitude decreases as the tensile stress increases, mainly due to an increase in vibrational frequencies of ZA modes with low ( 1 Å). This becomes important as the temperature is raised, but it is also appreciable in the low-temperature region, as shown in the figure. The lines derived from a HA are close to the results of the PIMD simulations at low temperature, but both sets of results depart progressively one from the other as temperature is raised. For , , but the opposite happens for 100 K. For relatively large tensile stresses of and eV/Å, we find in the whole temperature range presented in Fig. 6. For K, the difference between and steadily increases.

In Fig. 7 we display the mean-square displacements for as a function of applied stress for three temperatures: 100, 300, and 500 K. Symbols are results of PIMD simulations, whereas the lines correspond to the HA based on ZA and ZO vibrational modes. One observes first an important decrease in as the tensile stress is raised. This decrease is most appreciable for stresses in the range from 0 to eV Å. For larger stresses, the reduction of becomes slower. One also notices that the largest difference between and occurs for , and it becomes smaller for larger pressure. For eV Å, both sets of results cross each other and the data derived from PIMD become slightly larger than those corresponding to the HA.

To get better insight into the influence of nuclear quantum effects on for different pressures, we have plotted in Fig. 8 the mean-square displacements as derived from classical MD (open circles) and PIMD (solid circles) simulations at 25 K. The lines were obtained from the HA described above for the quantum case [see Eq. (14)], and for the classical calculation we used the expression

(16) |

One sees that the relevance of quantum effects increases for rising tensile stress, as can be measured from the ratio between the mean-square displacements in the quantum and classical case. In fact, this ratio goes from 1.3 for to 3.7 for eV Å. For a given stress, the difference between the quantum results for and the classical ones, , decreases as temperature is raised. In fact, at K the mean-square displacement derived from PIMD simulations is about 10% larger than the classical result in the stress range displayed in Fig. 8.

As indicated above, the mean-square displacement can be divided into two parts, and , the first one properly quantum in nature, measuring the extension of the quantum paths, and the second of a classical-like character, taking account of the centroid motion, i.e., global displacements of the paths. For the sake of comparing with the results of PIMD simulations, we have also calculated in the harmonic approximation , using Eqs. (14) and (16). To visualize the evolution of both terms as a function of temperature, we have plotted in Fig. 9 along with its two contributions and for = 96 and a stress eV Å. In the limit , vanishes and converges to a value of about Å. decreases for increasing temperature, as nuclear quantum effects become less relevant. On the contrary, the classical contribution increases with , linearly in the HA [see Eq. (16)], and almost linearly for the results of PIMD simulations.

For the system size shown in Fig. 9 (N = 96), both contributions to are equal at 64 K. At higher temperatures, the classical-like part is the main contribution to the atomic displacements in the direction. Values of given by PIMD simulations coincide within error bars with the mean-square atomic displacement obtained from classical MD simulations. The actual quantum delocalization can be estimated from the mean extension of the quantum paths in the direction, i.e., from . For eV Å, we find at 25 and 300 K an average extension 0.06 and 0.03 Å, respectively.

The picture displayed in Fig. 9 for the atom displacements in the out-of-plane direction is qualitatively the same for different system sizes and applied stresses, but the temperature region where or is the main contribution to as well as the crossing point greatly depend on both variables and . The dependence on is mainly caused by the enhancement of the classical-like contribution for increasing size, as observed earlier for results of classical MD simulations.Los et al. (2009); Gao and Huang (2014); Ramírez et al. (2016) The quantum contribution has a small size effect and converges fast as is increased.Herrero and Ramírez (2016) For a given system size , the ratio decreases for increasing (see Fig. 9), and there is a crossover temperature for which this ratio is unity. For classical-like motion dominates the atomic motion in the direction.

In Fig. 10 we present as a function of the applied stress for two system sizes: = 96 and 960. Symbols are results derived from PIMD simulations and lines correspond to crossover temperatures derived from the HA based on the ZA and ZO phonon bands. Below the lines we have , i.e. the quantum contribution dominates in , and the opposite happens above the lines with classical-like motion dominating the out-of-plane displacements. The crossover temperature depends on the system size , since the effective low-energy cutoff scales as . This means that for a given stress , an increase in causes the appearance of ZA modes with lower frequencies, which contribute to reduce (for a given , the classical behavior dominates for ).

The main limitations of the HA are the neglect of anharmonicity in the transverse vibrational modes, which is expected to be reasonably small at low , and the use of Eq. (15) for ZA mode frequencies under stress. This equation, giving as a sum of a zero-stress term, , and a term linear in the applied stress, , is expected to be valid for small- ZA modes, which dominate the mean-square displacement of C atoms in the transverse direction to the graphene plane. Taking into account these limitations, the harmonic model captures qualitatively, and almost quantitatively, the basic aspects of the competition between classical-like and quantum dynamics of the C atoms in the direction.

## Vi Summary

We have presented results of PIMD simulations of a graphene monolayer at several temperatures and tensile stresses. The importance of quantum effects has been quantified by comparing results of such quantum simulations with those obtained from classical MD simulations. Structural variables are found to change when quantum nuclear motion is taken into account, especially at low temperatures. Thus, the sheet area and interatomic distances change appreciably in the range of temperatures and stresses considered here.

The LCBOPII potential model was shown earlier to give a reliable description of structural and thermodynamic properties of graphene. We have investigated here its reliability to describe nuclear quantum effects in the presence of a tensile stress. The results obtained in the simulations have allowed us to analyze both the in-plane and real area of graphene as functions of and . The difference grows (decreases) as temperature (stress) is raised. The thermal contraction of becomes less important as the tensile stress increases and the amplitude of out-of-plane vibrations decreases. We emphasize that PIMD simulations yield a negative thermal expansion of at low for pressures so high as eV Å ( N/m). On the contrary, in classical calculations this thermal expansion becomes positive at much smaller stresses.

Zero-point expansion of the graphene layer due to nuclear quantum effects is not negligible, and it amounts to an increase of about 1% in the area . This zero-point effect is reduced for rising tensile stress. Moreover, the temperature dependence of the in-plane area is qualitatively different when derived from PIMD or classical simulations, even at temperatures between 300 and 1000 K. Such a difference appears for all tension stresses considered here, up to eV Å.

Atomic vibrations in the out-of-plane direction have been particularly considered, as they are important for the area and for the relative stability of the planar graphene layer vs. crumpling. Quantum effects are dominant for these vibrational modes provided that . However, the actual atomic motion at any finite temperature, resulting from the sum of mode contributions, is dominated by the classical thermal contribution as soon as the system size is large enough. This size effect appears in the quantum simulations at low temperatures, as a result of the presence of vibrational modes in the ZA band with smaller wavenumbers (frequencies) for larger graphene cells. In this respect, an interesting result is that the temperature region where quantum motion is dominant is enhanced by an external tensile stress, as shown in Fig. 10, i.e., for a given system size increases as the stress is raised.

An important point related to the consistency of the simulation results is their agreement with the principles of thermodynamics, in particular with the third law. This means that thermal expansion coefficients should converge to zero for . We have found that this requirement is verified by both the in-plane area and the real area obtained from PIMD simulations for and (tensile stress).

An analysis of graphene under tensile stresses larger than those considered here may be interesting for its effect on mechanical, electronic, and optical properties. Such a study can be hindered by limitations associated to effective potentials at large strains, and would require the use of ab-initio methods (e.g., density-functional theory). An efficient combination of these methods with path-integral simulations is still restricted to cell sizes relatively small, a limitation that is expected to be progressively overcome in the forthcoming years.

###### Acknowledgements.

The authors acknowledge the support of J. H. Los in the implementation of the LCBOPII potential. This work was supported by Dirección General de Investigación, MINECO (Spain) through Grants FIS2012-31713 and FIS2015-64222-C2.## Appendix A Pressure estimator

The Cartesian coordinates of the atomic nuclei in the crystalline membrane are denoted as , where the indicates the atom, and refers to the bead number. The staging coordinates are defined by a linear transformation of that diagonalizes the harmonic energy associated to the effective interactions between neighboring beads:Tuckerman and Hughes (1998); Tuckerman (2002)

(17) |

(18) |

The estimator of the 2D stress tensor is similar to that employed in previous work for 3D systems.(Ramírez et al., 2008; Herrero and Ramírez, 2014) For PIMD simulations, its components are given by expressions such asTuckerman and Hughes (1998)

(19) |

where the brackets indicate an ensemble average, is the dynamic mass associated to the staging coordinate , , are components of its corresponding velocity, and is an element of the 2D strain tensor. The masses are given byTuckerman and Hughes (1998); Tuckerman (2002); Herrero and Ramírez (2014)

(20) |

(21) |

where is the nuclear mass. The constant is given for each bead by for and ().

The stress used in our PIMD simulations is

(22) |

or with the estimator

(23) |

is the in-plane harmonic energy, that in terms of staging coordinates can be written as

(24) |

The last term in Eq. (23) is a sum of derivatives of the potential energy at different imaginary times (beads) . In the present case of graphene, these derivatives were carried out analytically.

## References

- Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nature Mater. 6, 183 (2007).
- Castro Neto et al. (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- Flynn (2011) G. W. Flynn, J. Chem. Phys. 135, 050901 (2011).
- Roldan et al. (2017) R. Roldan, L. Chirolli, E. Prada, J. Angel Silva-Guillen, P. San-Jose, and F. Guinea, Chem. Soc. Rev. 46, 4387 (2017).
- Ghosh et al. (2008) S. Ghosh, I. Calizo, D. Teweldebrhan, E. P. Pokatilov, D. L. Nika, A. A. Balandin, W. Bao, F. Miao, and C. N. Lau, Appl. Phys. Lett. 92, 151911 (2008).
- Nika et al. (2009) D. L. Nika, E. P. Pokatilov, A. S. Askerov, and A. A. Balandin, Phys. Rev. B 79, 155413 (2009).
- Balandin (2011) A. A. Balandin, Nature Mater. 10, 569 (2011).
- Lee et al. (2008) C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008).
- Seol et al. (2010) J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, et al., Science 328, 213 (2010).
- Prasher (2010) R. Prasher, Science 328, 185 (2010).
- Herrero and Ramírez (2016) C. P. Herrero and R. Ramírez, J. Chem. Phys. 145, 224701 (2016).
- Fasolino et al. (2007) A. Fasolino, J. H. Los, and M. I. Katsnelson, Nature Mater. 6, 858 (2007).
- Meyer et al. (2007) J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, Nature 446, 60 (2007).
- de Andres et al. (2012a) P. L. de Andres, F. Guinea, and M. I. Katsnelson, Phys. Rev. B 86, 144103 (2012a).
- Evans and Rawicz (1990) E. Evans and W. Rawicz, Phys. Rev. Lett. 64, 2094 (1990).
- Safran (1994) S. A. Safran, Statistical Thermodynamics of Surfaces, Interfaces, and Membranes (Addison Wesley, New York, 1994).
- Cerda and Mahadevan (2003) E. Cerda and L. Mahadevan, Phys. Rev. Lett. 90, 074302 (2003).
- Wong and Pellegrino (2006) Y. W. Wong and S. Pellegrino, J. Mech. Materials Struct. 1, 3 (2006).
- Kirilenko et al. (2011) D. A. Kirilenko, A. T. Dideykin, and G. Van Tendeloo, Phys. Rev. B 84, 235417 (2011).
- Nicholl et al. (2015) R. J. T. Nicholl, H. J. Conley, N. V. Lavrik, I. Vlassiouk, Y. S. Puzyrev, V. P. Sreenivas, S. T. Pantelides, and K. I. Bolotin, Nature Commun. 6, 8789 (2015).
- Ruiz-Vargas et al. (2011) C. S. Ruiz-Vargas, H. L. Zhuang, P. Y. Huang, A. M. van der Zande, S. Garg, P. L. McEuen, D. A. Muller, R. G. Hennig, and J. Park, Nano Lett. 11, 2259 (2011).
- Kosmrlj and Nelson (2013) A. Kosmrlj and D. R. Nelson, Phys. Rev. E 88, 012136 (2013).
- Kosmrlj and Nelson (2014) A. Kosmrlj and D. R. Nelson, Phys. Rev. E 89, 022126 (2014).
- Ramírez and Herrero (2017) R. Ramírez and C. P. Herrero, Phys. Rev. B 95, 045423 (2017).
- López-Polín et al. (2017) G. López-Polín, M. Jaafar, F. Guinea, R. Roldán, C. Gómez-Navarro, and J. Gómez-Herrero, Carbon 124, 42 (2017).
- Pozzo et al. (2011) M. Pozzo, D. Alfè, P. Lacovig, P. Hofmann, S. Lizzit, and A. Baraldi, Phys. Rev. Lett. 106, 135501 (2011).
- Woods et al. (2014) C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, et al., Nature Phys. 10, 451 (2014).
- Amorim et al. (2016) B. Amorim, A. Cortijo, F. de Juan, A. G. Grushine, F. Guinea, A. Gutierrez-Rubio, H. Ochoa, V. Parente, R. Roldan, P. San-Jose, et al., Phys. Rep. 617, 1 (2016).
- Shimojo et al. (2008) F. Shimojo, R. K. Kalia, A. Nakano, and P. Vashishta, Phys. Rev. B 77, 085103 (2008).
- de Andres et al. (2012b) P. L. de Andres, F. Guinea, and M. I. Katsnelson, Phys. Rev. B 86, 245409 (2012b).
- Chechin et al. (2014) G. M. Chechin, S. V. Dmitriev, I. P. Lobzenko, and D. S. Ryabov, Phys. Rev. B 90, 045432 (2014).
- Herrero and Ramírez (2009) C. P. Herrero and R. Ramírez, Phys. Rev. B 79, 115429 (2009).
- Cadelano et al. (2009) E. Cadelano, P. L. Palla, S. Giordano, and L. Colombo, Phys. Rev. Lett. 102, 235502 (2009).
- Akatyeva and Dumitrica (2012) E. Akatyeva and T. Dumitrica, J. Chem. Phys. 137, 234702 (2012).
- Lee et al. (2013) G.-D. Lee, E. Yoon, N.-M. Hwang, C.-Z. Wang, and K.-M. Ho, Appl. Phys. Lett. 102, 021603 (2013).
- Shen et al. (2013) H.-S. Shen, Y.-M. Xu, and C.-L. Zhang, Appl. Phys. Lett. 102, 131905 (2013).
- Los et al. (2009) J. H. Los, M. I. Katsnelson, O. V. Yazyev, K. V. Zakharchenko, and A. Fasolino, Phys. Rev. B 80, 121405 (2009).
- Ramírez et al. (2016) R. Ramírez, E. Chacón, and C. P. Herrero, Phys. Rev. B 93, 235419 (2016).
- Magnin et al. (2014) Y. Magnin, G. D. Foerster, F. Rabilloud, F. Calvo, A. Zappelli, and C. Bichara, J. Phys.: Condens. Matter 26, 185401 (2014).
- Los et al. (2016) J. H. Los, A. Fasolino, and M. I. Katsnelson, Phys. Rev. Lett. 116, 015901 (2016).
- Gillan (1988) M. J. Gillan, Phil. Mag. A 58, 257 (1988).
- Ceperley (1995) D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- Brito et al. (2015) B. G. A. Brito, L. Cândido, G.-Q. Hai, and F. M. Peeters, Phys. Rev. B 92, 195416 (2015).
- Mounet and Marzari (2005) N. Mounet and N. Marzari, Phys. Rev. B 71, 205214 (2005).
- Shao et al. (2012) T. Shao, B. Wen, R. Melnik, S. Yao, Y. Kawazoe, and Y. Tian, J. Chem. Phys. 137, 194901 (2012).
- Gao and Huang (2014) W. Gao and R. Huang, J. Mech. Phys. Solids 66, 42 (2014).
- Feynman (1972) R. P. Feynman, Statistical Mechanics (Addison-Wesley, New York, 1972).
- Kleinert (1990) H. Kleinert, Path Integrals in Quantum Mechanics, Statistics and Polymer Physics (World Scientific, Singapore, 1990).
- Chandler and Wolynes (1981) D. Chandler and P. G. Wolynes, J. Chem. Phys. 74, 4078 (1981).
- Herrero and Ramírez (2014) C. P. Herrero and R. Ramírez, J. Phys.: Condens. Matter 26, 233201 (2014).
- Los et al. (2005) J. H. Los, L. M. Ghiringhelli, E. J. Meijer, and A. Fasolino, Phys. Rev. B 72, 214102 (2005).
- Ghiringhelli et al. (2005a) L. M. Ghiringhelli, J. H. Los, A. Fasolino, and E. J. Meijer, Phys. Rev. B 72, 214103 (2005a).
- Zakharchenko et al. (2009) K. V. Zakharchenko, M. I. Katsnelson, and A. Fasolino, Phys. Rev. Lett. 102, 046808 (2009).
- Ghiringhelli et al. (2005b) L. M. Ghiringhelli, J. H. Los, E. J. Meijer, A. Fasolino, and D. Frenkel, Phys. Rev. Lett. 94, 145701 (2005b).
- Politano et al. (2012) A. Politano, A. R. Marino, D. Campi, D. Farías, R. Miranda, and G. Chiarello, Carbon 50, 4903 (2012).
- Lambin (2014) P. Lambin, Appl. Sci. 4, 282 (2014).
- Sfyris et al. (2015) D. Sfyris, E. N. Koukaras, N. Pugno, and C. Galiotis, J. Appl. Phys. 118, 075301 (2015).
- Memarian et al. (2015) F. Memarian, A. Fereidoon, and M. D. Ganji, Superlattices Microstruct. 85, 348 (2015).
- Zou et al. (2016) J.-H. Zou, Z.-Q. Ye, and B.-Y. Cao, J. Chem. Phys. 145, 134705 (2016).
- Anastasi et al. (2016) A. A. Anastasi, K. Ritos, G. Cassar, and M. K. Borg, Mol. Simul. 42, 1502 (2016).
- Ghasemi and Rajabpour (2017) H. Ghasemi and A. Rajabpour, J. Phys. Conf. Series 785, 012006 (2017).
- Tuckerman et al. (1992) M. E. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys. 97, 1990 (1992).
- Tuckerman and Hughes (1998) M. E. Tuckerman and A. Hughes, in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B. J. Berne, G. Ciccotti, and D. F. Coker (Word Scientific, Singapore, 1998), p. 311.
- Martyna et al. (1999) G. J. Martyna, A. Hughes, and M. E. Tuckerman, J. Chem. Phys. 110, 3275 (1999).
- Tuckerman (2002) M. E. Tuckerman, in Quantum Simulations of Complex Many–Body Systems: From Theory to Algorithms, edited by J. Grotendorst, D. Marx, and A. Muramatsu (NIC, FZ Jülich, 2002), p. 269.
- Herman et al. (1982) M. F. Herman, E. J. Bruskin, and B. J. Berne, J. Chem. Phys. 76, 5150 (1982).
- Herrero et al. (2006) C. P. Herrero, R. Ramírez, and E. R. Hernández, Phys. Rev. B 73, 245211 (2006).
- Herrero and Ramírez (2011) C. P. Herrero and R. Ramírez, J. Chem. Phys. 134, 094510 (2011).
- Ramírez et al. (2012) R. Ramírez, N. Neuerburg, M. V. Fernández-Serra, and C. P. Herrero, J. Chem. Phys. 137, 044502 (2012).
- Ramírez et al. (2008) R. Ramírez, C. P. Herrero, E. R. Hernández, and M. Cardona, Phys. Rev. B 77, 045210 (2008).
- Fournier and Barbetta (2008) J.-B. Fournier and C. Barbetta, Phys. Rev. Lett. 100, 078103 (2008).
- Imparato (2006) A. Imparato, J. Chem. Phys. 124, 154714 (2006).
- Waheed and Edholm (2009) Q. Waheed and O. Edholm, Biophys. J. 97, 2754 (2009).
- Chacón et al. (2015) E. Chacón, P. Tarazona, and F. Bresme, J. Chem. Phys. 143, 034706 (2015).
- Nicholl et al. (2017) R. J. T. Nicholl, N. V. Lavrik, I. Vlassiouk, B. R. Srijanto, and K. I. Bolotin, Phys. Rev. Lett. 118, 266101 (2017).
- Hahn et al. (2016) K. R. Hahn, C. Melis, and L. Colombo, J. Phys. Chem. C 120, 3026 (2016).
- Michel et al. (2015) K. H. Michel, S. Costamagna, and F. M. Peeters, Phys. Status Solidi B 252, 2433 (2015).
- Landau and Lifshitz (1980) L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, Oxford, 1980), 3rd ed.
- Herrero (2002) C. P. Herrero, Phys. Rev. B 65, 014112 (2002).
- Herrero and Ramírez (1995) C. P. Herrero and R. Ramírez, Phys. Rev. B 51, 16761 (1995).
- Landau and Lifshitz (1965) L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Pergamon, Oxford, 1965), 2nd ed.
- Gillan (1990) M. J. Gillan, in Computer Modelling of Fluids, Polymers and Solids, edited by C. R. A. Catlow, S. C. Parker, and M. P. Allen (Kluwer, Dordrecht, 1990), p. 155.
- Karssemeijer and Fasolino (2011) L. J. Karssemeijer and A. Fasolino, Surf. Sci. 605, 1611 (2011).
- Wirtz and Rubio (2004) L. Wirtz and A. Rubio, Solid State Commun. 131, 141 (2004).