# Thermal properties of graphene under tensile stress

###### Abstract

Thermal properties of graphene display peculiar characteristics associated to the two-dimensional nature of this crystalline membrane. These properties can be changed and tuned in the presence of applied stresses, both tensile and compressive. Here we study graphene monolayers under tensile stress by using path-integral molecular dynamics (PIMD) simulations, which allows one to take into account quantization of vibrational modes and analyze the effect of anharmonicity on physical observables. The influence of the elastic energy due to strain in the crystalline membrane is studied for increasing tensile stress and for rising temperature (thermal expansion). We analyze the internal energy, enthalpy, and specific heat of graphene, and compare the results obtained from PIMD simulations with those given by a harmonic approximation for the vibrational modes. This approximation turns out to be precise at low temperatures, and deteriorates as temperature and pressure are increased. At low temperature the specific heat changes as for stress-free graphene, and evolves to a dependence as the tensile stress is increased. Structural and thermodynamic properties display nonnegligible quantum effects, even at temperatures higher than 300 K. Moreover, differences in the behavior of the in-plane and real areas of graphene are discussed, along with their associated properties. These differences show up clearly in the corresponding compressibility and thermal expansion coefficient.

###### pacs:

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

Two-dimensional materials, and graphene in particular, have attracted in last years a great deal of attention in the scientific community, due to their peculiar electronic, elastic, and thermal properties.Geim and Novoselov (2007); Roldan et al. (2017); Castro Neto et al. (2009) Thus, graphene presents high values of the thermal conductivity,Balandin (2011); Ghosh et al. (2008) and large in-plane elastic constants.Lee et al. (2008) Moreover, its mechanical properties are interesting for possible applications, as cooling of electronic devices.Prasher (2010); Seol et al. (2010)

The ideal structure of pure defect-free graphene is a planar honeycomb lattice, but departures from this flat configuration may appreciably alter its atomic-scale and macroscopic properties.Meyer et al. (2007) Several reasons can cause bending of a graphene sheet, such as the presence of defects and external stresses.Fasolino et al. (2007) In addition, thermal fluctuations at finite temperatures give rise to out-of-plane displacements of the C atoms, and for , zero-point motion yields also a departure of strict planarity of the graphene sheet.Herrero and Ramírez (2016)

The influence of strain in several characteristics of two-dimensional (2D) materials, such as graphene and metallic dichalcogenides, has been emphasized in recent years. This includes electronic transport, optical properties, and the formation of moiré patterns.Woods et al. (2014); Amorim et al. (2016) From a structural viewpoint, external stresses may cause significant changes in crystalline membranes, which can crumple in the presence of a compressive stress, as is well known for lipid membranesEvans and Rawicz (1990); Safran (1994) and polymer films.Cerda and Mahadevan (2003); Wong and Pellegrino (2006) For graphene, crumpling was observed in supported as well as freestanding samples, and has been explained as due to both out-of-plane phonons and static wrinkling.Kirilenko et al. (2011); Nicholl et al. (2015) Recent molecular dynamics simulations indicate that the maximum compressive stress that a freestanding graphene sheet can sustain without crumpling decreases as the system size grows, and was estimated to be about 0.1 N/m at room temperature in the thermodynamic limit.Ramírez and Herrero (2017)

A tensile stress in the graphene plane does not break the planarity of the sheet, but gives rise to significant variations in the elastic properties of the material.López-Polín et al. (2017) For example, the in-plane Young modulus increases by a factor of three for a tensile stress of 1 N/m.Ramírez and Herrero (2017) In this respect, the real area per atom, , can be crucial to understand the elastic properties, which have been in the past usually referred to its projection, , onto the mean plane of the membrane (). This question has been discussed along the years for biological membranes,Helfrich and Servuss (1984); Safran (1994); Fournier and Barbetta (2008); Tarazona et al. (2013); Chacón et al. (2015) and has been recently examined for crystalline membranes such as graphene.Pozzo et al. (2011); Nicholl et al. (2015, 2017); Ramírez et al. (2016); Herrero and Ramírez (2016) In particular, the high-quality data obtained by Nicholl et al.Nicholl et al. (2015, 2017) clearly indicate that some experimental techniques can measure properties related to the real area , whereas other techniques may be used to quantify variables connected to the projected area . The difference has been called hidden area for graphene in Ref. Nicholl et al., 2017, as well as excess area for biological membranes.Helfrich and Servuss (1984); Fournier and Barbetta (2008) A precise knowledge of the behavior of both areas is important to clarify the temperature and stress dependence of structural and thermodynamic properties of graphene. Thus, it is possible to define the elastic properties in relation to the area or to , which may behave very differently. In fact, it is known that referring to the in-plane area , one finds a negative thermal expansion coefficient, but for the real area the thermal expansion is positive.Pozzo et al. (2011); Herrero and Ramírez (2016)

A deep comprehension of the properties of 2D systems has been for many years a persistent goal in statistical physics.Safran (1994); Nelson et al. (2004); Tarazona et al. (2013) This has been in part due to the complexity of the considered systems, such as biological membranes and soft condensed matter.Tarazona et al. (2013); Chacón et al. (2015) In this respect, graphene can be dealt with as a model system where descriptions at an atomic level can be connected with physical properties of the material. Thus, thermal properties of graphene have been investigated in recent years,Amorim et al. (2014); Pop et al. (2012); Wang et al. (2016); Fong et al. (2013); Mann et al. (2016) and in particular its thermal expansion and heat conduction were studied by various theoretical and experimental techniques.Balandin (2011); Pop et al. (2012); Ma et al. (2012); Alofi and Srivastava (2013); da Silva et al. (2014); Burmistrov et al. (2016)

Several theoretical works carried out to study thermodynamic properties of graphene (e.g., specific heat, thermal expansion, …), were based on density-functional-theory (DFT) calculations combined with a quantum quasi-harmonic approximation (QHA) for the vibrational modes.Mounet and Marzari (2005); Sevik (2014); Mann et al. (2017) This is expected to yield reliable results at low temperature, but may be questioned at relatively high temperature, especially for the thermal expansion, due to an important anharmonic coupling between in-plane and out-of-plane modes, not included in the QHA. Moreover, classical Monte Carlo and molecular dynamics simulations based on ab-initio,Shimojo et al. (2008); Pozzo et al. (2011); de Andres et al. (2012); 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 potentialsFasolino et al. (2007); Ramírez et al. (2016); Magnin et al. (2014); Los et al. (2016) can give reliable results at relatively high temperature, but fail to describe thermodynamic properties at , with 1000 K the Debye temperature of the material. This means, in particular, that room-temperature results obtained for some properties of graphene from classical atomistic simulations may be clearly unrealistic. These shortcuts may be overcome by using simulation methods which explicitly include nuclear quantum effects, in particular those based on Feynman path integrals.Gillan (1988); Ceperley (1995); Brito et al. (2015); Herrero and Ramírez (2016)

Here we use the path-integral molecular dynamics (PIMD) method to study thermal properties of graphene under tensile stress at temperatures between 12 and 2000 K. The thermal behavior of the graphene surface is studied, considering the difference between in-plane and real areas. The in-plane thermal expansion coefficient turns out to be negative at low temperatures, with a crossing to positive values at a temperature which decreases fast as tensile stress is raised. Particular emphasis is laid on the temperature dependence of the specific heat at low , for which results of the simulations are compared with predictions based on harmonic vibrations of the crystalline membrane. This approximation happens to be noticeably accurate at low temperatures, once the frequencies of out-of-plane ZA modes are properly renormalized for changing applied stress.

The paper is organized as follows. In Sec. II we describe the computational method employed in the simulations. Results for the internal energy and enthalpy of graphene are given in Sec. III. The thermal expansion is discussed in Sec. IV. In Sec. V we present results for the specific heat, whereas in Sec. VI the compressibility of graphene is discussed. In Sec. VII we summarize the main results.

## Ii Computational Method

### ii.1 Path-integral molecular dynamics

We use PIMD simulations to study equilibrium properties of graphene monolayers as a function of temperature and pressure. The PIMD procedure is based on the Feynman path-integral formulation of statistical mechanics,Feynman (1972) an adequate nonperturbative approach to study finite-temperature properties of many-body quantum systems. In the applications of this method to numerical simulations, each quantum particle is represented by a set of (Trotter number) replicas (or beads), that behave as classical-like particles forming a ring polymer.Gillan (1988); Ceperley (1995) Thus, one deals with a classical isomorph whose dynamics is artificial, since it does not reflect the real quantum dynamics of the actual particles, but is useful for effectively sampling the many-body configuration space, yielding precise results for time-independent equilibrium properties of the quantum system. Details on this type of simulation techniques are given in Refs. Gillan, 1988; Ceperley, 1995; Herrero and Ramírez, 2014; Cazorla and Boronat, 2017.

The interatomic interactions between carbon atoms are described here by using the LCBOPII effective potential, a long-range bond order potential, which has been mainly employed to carry out classical simulations of carbon-based systems.Los et al. (2005) In particular, it was used to study the phase diagram of carbon, including graphite, diamond, and the liquid, and showed its accuracy in predicting rather precisely the graphite-diamond transition line.Ghiringhelli et al. (2005) In recent years, this interatomic potential has been also found to describe well several properties of graphene,Fasolino et al. (2007); Los et al. (2016) and its Young’s modulus in particular.Zakharchenko et al. (2009); Politano et al. (2012); Ramírez and Herrero (2017) This interatomic potential was lately employed to carry out PIMD simulations, which allowed to quantify quantum effects in graphene monolayers by comparing with results of classical simulations,Herrero and Ramírez (2016) as well as to study thermodynamic properties of this 2D material.Herrero and Ramírez (2018) Here, in line with earlier simulations,Ramírez et al. (2016); Herrero and Ramírez (2016); Ramírez and Herrero (2017) the original parameterization of the LCBOPII potential has been slightly changed to rise the zero-temperature bending constant from 1.1 eV to 1.49 eV, a value close to experimental data.Lambin (2014)

The calculations presented here have been performed in the isothermal-isobaric ensemble, where one fixes the number of carbon atoms (), the applied stress (), and the temperature (). The stress in the reference plane of graphene, with units of force per unit length, coincides with the so-called mechanical or frame tension.Ramírez and Herrero (2017); Fournier and Barbetta (2008); Shiba et al. (2016) is obtained in the simulations from the stress tensor , whose components are given by expressions such asTuckerman and Hughes (1998); Herrero and Ramírez (2018)

(1) |

where the brackets indicate an ensemble average and are staging coordinates,Tuckerman et al. (1993) with and . In Eq. (1), is the dynamic mass associated to and , are components of its corresponding velocity. The constant is given by for and . Here , is the instantaneous potential energy, and is an element of the 2D strain tensor. The stress , conjugate to the in-plane area , is given by the trace of the stress tensor:

(2) |

Effective algorithms have been employed for the PIMD simulations, as those described in the literature.Tuckerman et al. (1992); Tuckerman and Hughes (1998); Martyna et al. (1999) In particular, a constant temperature was accomplished by coupling chains of four Nosé-Hoover thermostats, and an additional chain of four barostats was coupled to the area of the graphene simulation box to yield the required stress .Tuckerman and Hughes (1998); Herrero and Ramírez (2014) The equations of motion were integrated by employing the reversible reference system propagator algorithm (RESPA), allowing to define different time steps for the integration of the fast and slow degrees of freedom.Martyna et al. (1996) The kinetic energy has been calculated by using the virial estimator, which displays a statistical uncertainty significantly smaller than the potential energy, .Tuckerman and Hughes (1998) The time step associated to the interatomic forces has been taken as 0.5 fs, which was found to be suitable for the carbon atomic mass and the temperature range considered here. More details on this kind of PIMD simulations can be found elsewhere.Tuckerman and Hughes (1998); Herrero et al. (2006); Herrero and Ramírez (2011)

Rectangular simulation cells with = 960 atoms have been considered, with similar side lengths in the and directions ( Å), and periodic boundary conditions were assumed. Sampling of the configuration space was performed in the temperature range between 12 K and 2000 K. A temperature-dependent Trotter number has been defined as , which yields a roughly constant precision for the PIMD results at different temperatures.Herrero et al. (2006); Herrero and Ramírez (2011); Ramírez et al. (2012) Given a temperature and a stress , a typical simulation run included PIMD steps for system equilibration and steps for calculation of average properties.

### ii.2 In-plane vs real area

As explained above, in the isothermal-isobaric ensemble used here we fix the applied stress in the plane, allowing fluctuations in the in-plane area of the simulation cell for which periodic boundary conditions are applied. Carbon atoms are free to move in the coordinate (out-of-plane direction), and in general any measure of the real surface of a graphene sheet at should give a value larger than the in-plane area. In this respect, it has been discussed for biological membranes that their properties should be described using the notion of a real surface rather than a projected (in-plane) surface.Waheed and Edholm (2009); Chacón et al. (2015) A similar question has been also recently posed 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 may be relevant for addressing the calculation of thermodynamic variables, as the in-plane area is the variable conjugate to the stress used in our simulations. The real area (also called true, actual, or effective area in the literatureFournier and Barbetta (2008); Waheed and Edholm (2009); Chacón et al. (2015)) is conjugate to the usually-called surface tension.Safran (1994)

The in-plane area has been most used to present the results of atomistic simulations of graphene layers.Gao and Huang (2014); Brito et al. (2015); Zakharchenko et al. (2009); Los et al. (2016); Chacón et al. (2015) For biological membranes, however, it has been shown that values of the compressibility may significantly differ when they are related to or to , and something similar has been recently found for the elastic properties of graphene from classical molecular dynamics simulations.Ramírez and Herrero (2017)

Here we calculate a real area in three-dimensional (3D) space by a triangulation based on the actual positions of the C atoms along a simulation run (in fact we use the beads associated to the atomic nuclei). Specifically, is obtained from a sum of areas corresponding to the structural hexagons. 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)

The instantaneous area per atom for imaginary time (bead) is given by

(3) |

where is the area of triangle in hexagon , and the sum in is extended to the hexagons in a cell containing carbon atoms. Here, the triangles are defined with the coordinates corresponding to bead of atom . The area is then calculated as

(4) |

It is clear that coincides with for strictly planar graphene, and in general one has . Both areas display temperature dependencies qualitatively different. In fact, does not present a negative thermal expansion, as happens for the in-plane area in a large temperature rangeZakharchenko et al. (2009); Herrero and Ramírez (2016) (see below). The difference increases with temperature, as bending of the graphene sheet increases, but even for , and are not exactly equal, due to zero-point motion of the C atoms in the transverse direction.

## Iii Energy and enthalpy

### iii.1 Internal energy

The internal energy is obtained from the results of our PIMD simulations as a sum of the kinetic, , and potential energy, , at a given temperature. The kinetic energy has been calculated by employing the virial estimator,Tuckerman and Hughes (1998) which displays a statistical uncertainty smaller than the potential energy. We have included in the center-of-mass translational energy, a classical magnitude amounting to at temperature . This quantity is irrelevant for the energy per atom in large systems, but we have included it to minimize finite size effects.Herrero and Ramírez (2016)

We express the internal energy as , taking as energy reference the value corresponding to the equilibrium configuration of a planar graphene surface in a classical approach at (minimum-energy configuration of the considered LCBOPII potential, without quantum atomic delocalization). This corresponds to an interatomic distance Å, i.e., an area Å per atom. In a quantum approach, out-of-plane atomic fluctuations associated to zero-point motion appear even for , and the graphene layer is not strictly planar. Moreover, anharmonicity of in-plane vibrations gives rise to a zero-point lattice expansion, yielding a distance Å, i.e., around 1% larger than the classical distance at .Herrero and Ramírez (2016)

In Fig. 1 we show the temperature dependence of the internal energy per atom, , as derived from our PIMD simulations in the isothermal-isobaric ensemble for = 0 (circles), (squares), and eV Å (diamonds). For it was shown earlier that the size effect on the internal energy per atom is negligible for = 960 (less than the symbol size in Fig. 1), as compared with the largest cells considered in Ref. Herrero and Ramírez, 2016. The zero-point energy, , is found to be 172 meV/atom for = 0, and slightly higher for eV Å, but it appreciably increases for eV Å to a value of 207 meV/atom. We will show later that this rise in internal energy is basically due to the elastic energy associated to an increase in the area . For comparison with our results, we also present in Fig. 1 the energy obtained by Brito et al.Brito et al. (2015) from path-integral Monte Carlo simulations using the REBO potential (dashed line). These authors called it vibrational energy, and corresponds to our internal energy; our vibrational energy is defined below (see Sec. III.C). The REBO potential yields a zero-point energy of 0.181 eV/atom, i.e., a 5% higher than that found here with the LCBOPII potential.

The rise in internal energy with an applied tensile stress is presented in Fig. 2 for three temperatures: = 25 K (circles), 300 K (squares), and 500 K (diamonds). For zero stress we find at 25 K an internal energy very close to that of the ground state . In the stress range displayed in Fig. 2 the internal energy can be fitted to an expression , with a coefficient of the quadratic term Å eV, nearly independent of the temperature. This is the leading term for the variation of internal energy with the applied pressure, which is nearly parabolic for stresses between 0 and eV Å.

The change in internal energy with temperature and pressure can be described from variations in the elastic and vibrational contributions to . This is analyzed in the following sections.

### iii.2 Elastic energy

An important part of the internal energy corresponds to the elastic energy due to changes in the area of graphene. It is directly related to the actual interatomic distance, i.e., to the real area , rather than to the in-plane area (or in-plane strain ). Then, the internal energy at temperature can be written asHerrero and Ramírez (2016)

(5) |

where is the elastic energy for an area , and is the vibrational energy of the system. The area is a function of the stress and temperature , but this is not explicitly indicated in Eq. (5) for simplicity of the notation. One expects non-zero values of the elastic energy, even in the absence of an externally applied stress, due to thermal expansion at finite temperatures, as well as for zero-point expansion at . is given by contributions of phonons in graphene, both in-plane and out-of-plane vibrational modes.

We define the elastic energy corresponding to an area as the increase in energy of a strictly planar graphene layer with respect to the minimum energy (for an area = 2.61888 Å/atom). By definition , and for small changes of we found that it follows a dependence , with = 2.41 eV Å. This dependence for yields a 2D bulk modulus = , i.e., = 12.6 eV Å, in agreement with the result found earlier in classical calculations at .Ramírez and Herrero (2017)

Our PIMD simulations directly yield , which can be split into an elastic and a vibrational part, as in Eq. (5). 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).

In Fig. 3 we display the temperature dependence of the elastic energy, as derived from our PIMD simulations for , , and eV Å. For a given external stress, increases with , as a consequence of the thermal expansion of the real area , which turns out to be positive for all temperatures (, see below). Note that, in contrast, the in-plane area displays a thermal contraction () in a wide temperature range, so that it is not a suitable candidate for a reliable definition of the elastic energy.

For = 0 we find a positive elastic energy in the zero-temperature limit, = 1.7 meV/atom, due to zero-point lattice expansion, which causes that . This low- limit appreciably increases for due to the stress-induced increase in area , which yields values of 12 and 52 meV/atom for the elastic energy at and eV Å, respectively. For finite temperatures, one observes in Fig. 3 that increases faster with temperature for larger tensile stress. This is due to an increase in the graphene compressibility (reduction of the elastic constants) for rising tensile stress, which in turn causes an increase in the thermal expansion coefficient (see Sec. IV).

### iii.3 Vibrational energy

After calculating the elastic energy for an area resulting from PIMD simulations at given and , we obtain the vibrational energy by subtracting the elastic energy from the internal energy: (see Eq. (5)). In Fig. 4 we present the temperature dependence of the vibrational energy of graphene for (circles), (squares), and eV Å (diamonds), as derived from our simulations. In contrast to Figs. 1 and 3, where one observes that and increase with the applied tensile stress , in Fig. 4 one sees that the vibrational energy is lower for higher tensile stress. This is mainly due to a decrease in the vibrational frequency of in-plane modes for increasing stress (increasing area), corresponding to positive Grünesien parameters.Ashcroft and Mermin (1976); Mounet and Marzari (2005)

The zero-point vibrational energy is found to decrease from 170 meV/atom for to 164 and 155 meV/atom for and eV Å, respectively. This decrease, although clearly noticeable, is smaller than the increase of about 50 meV/atom in the elastic energy for eV Å (see Fig. 3). These zero- energy values per atom correspond in a harmonic approximation to , which is a mean value for the frequencies in the six phonon bands of graphene.Karssemeijer and Fasolino (2011); Wirtz and Rubio (2004) Our results for correspond to cm, to be compared with 833 cm for eV Å.

The three curves for shown in Fig. 4 for different stresses are roughly parallel in the displayed temperature range. Values of presented in this figure are clearly larger than those corresponding to a classical model for the vibrational modes. In this limit, the vibrational energy per atom is , which means 39 and 78 meV/atom for and 600 K, respectively. In contrast, we find = 166 and 206 meV/atom from PIMD simulations for eV Å at those temperatures.

We note that in the language of membranes and 2D elastic media there appears an energy contribution due to bending of the surface, that is usually taken into account through the bending constant , which measures the rigidity of the membrane. In our present formulation, the bending energy is included in the vibrational energy associated to the flexural ZA modes, as can be seen in their contribution to the specific heat (see Sec. V). Anharmonic couplings between in-plane and out-of-plane modes are also expected to show up, especially at high temperatures.Adamyan et al. (2016)

### iii.4 Enthalpy

Since we are working in the isothermal-isobaric ensemble, it is natural to consider the enthalpy as a relevant thermodynamic variable, in particular to calculate the specific heat of graphene at several applied stresses (see Sec. V). For a given stress , we define the enthalpy as . As a reference we will consider , such that converges in the classical limit to at .

In Fig. 5(a) we present the enthalpy of graphene, , as a function of temperature for a tensile stress eV Å. Results of our simulations are displayed as circles. Taking into account that the enthalpy can be written as

(6) |

we have plotted in Fig. 5(a) the temperature dependence of the three contributions in the r.h.s. of Eq. (6). The largest change with temperature appears for the vibrational energy (diamonds). The elastic energy, , as well as the term are also found to depend on temperature, due to the thermal expansion (or contraction) of the in-plane area , but their change is much smaller than that of . As a result, the enthalpy turns out to be smaller than the vibrational energy at low temperature, because . However, becomes larger than at K, due to the larger increase in for rising temperature.

To better appreciate the low-temperature region, in Fig. 5(b) we present the temperature dependence of the enthalpy, as derived from PIMD simulations up to K. Symbols indicate results of the simulations for = 0 (circles), (squares), and eV Å (diamonds) Dashed lines are fits to the expression in the temperature region from to 100 K. This polynomial form shows good agreement with the temperature dependence of the enthalpy obtained from the simulations in the fitted region. For K, the lines depart progressively from the data points. Note that a linear term in this expression for the enthalpy is not allowed for thermodynamic consistency, since the specific heat has to vanish in the limit . The coefficient changes from eV K for to eV K for eV Å. The coefficient varies from 3.5 to eV K in the same stress range. This is important for the low-temperature dependence of the specific heat discussed in Sec. V, as the contribution of the quadratic coefficient decreases for increasing stress, whereas is found to increase. This means that the specific heat in the lowest temperatures accessible to our simulation method will be given as a combination of a linear term and a quadratic one , with the latter becoming more important for larger tensile stresses.

## Iv Thermal expansion

In the low-temperature limit (), the real area and in-plane area derived from PIMD simulations converge to 2.6459 Å/atom and 2.6407 Å/atom, respectively. Comparing these values with the classical minimum Å, we find a zero-point expansion in and of about 0.02 Å/atom ( 1%), due to an increase in the mean C–C bond length (an anharmonic effect). There also appears a difference of 0.2% between real and in-plane areas at low temperature, caused by out-of-plane zero-point motion, so that the layer is not strictly planar. This is a genuine quantum effect, as in classical simulations for one finds a planar layer in which and coincide.Herrero and Ramírez (2016); Ramírez and Herrero (2017) The difference increases as temperature is raised, since is the projection of on the reference plane, and the real graphene surface becomes increasingly bent for rising temperature because of larger out-of-plane atomic displacements.

For , the area displays an almost constant value up to K, and increases at higher temperatures. However, decreases in the temperature range from to K, reaches a minimum and then increases at higher .Herrero and Ramírez (2016) These results for are qualitatively similar to those found from classical Monte Carlo and molecular dynamics simulations of graphene,Zakharchenko et al. (2009); Gao and Huang (2014) but in PIMD simulations the contraction of with respect to the zero-temperature value is significantly larger than for classical calculations.

To analyze the thermal expansion of graphene, and following our definitions of the areas and , we will consider two different thermal expansion coefficients. The first of them, , refers to changes in the real area:

(7) |

This coefficient takes mainly into account changes in the interatomic distances, and is rather insensitive to bending of the graphene layer. The second coefficient, , is a measure of variations in the in-plane area and has been widely used in the literature to analyze results of classical simulations and analytical calculations:

(8) |

In Fig. 6 we present both thermal expansion coefficients, as derived from our PIMD simulations for (circles), (squares), and eV Å (diamonds). Symbols are data points obtained from numerical derivatives of the areas (panel a) and (panel b). For these derivatives we took temperature intervals ranging from 10 K at low temperature to about 50 K at temperatures K. In general, the statistical uncertainty (error bars) in the values of obtained from numerical derivatives is larger than that found for , as a consequence of the larger fluctuations in . The behavior of shown in Fig. 6(a) is similar to that observed for most 3D materials, i.e., it goes to zero in the low-temperature limit and increases for rising temperature so that for .Ashcroft and Mermin (1976); Herrero and Ramírez (2000) This is related to the fact that the in-plane vibrational modes have positive Grüneisen parameters when they are calculated with respect to the real area . One observes in Fig. 6(a) that is higher for larger tensile stress. In fact, at 800 K the value obtained for eV Å is 40% larger than that corresponding to . Viewing the thermal expansion as an anharmonic effect, this can be interpreted as an increase in anharmonicity for larger tensile stress. This is associated with an increase in area and the corresponding decrease in the frequency of vibrational in-plane modes, which in turn causes a larger vibrational amplitude and consequently a larger anharmonicity.

The in-plane thermal expansion coefficient also converges to zero in the low-temperature limit, but contrary to it decreases for increasing temperature until reaching a minimum, as shown in Fig 6(b). The negative value of in the minimum approaches zero for rising tensile stress. In fact, it goes from K for to K for eV Å. At higher , approaches zero and eventually becomes positive at a temperature , which changes appreciably with the applied stress. As a result, we find , 590 K, and 360() K, for the three values of the stress presented in Fig. 6.

It is interesting to note that the difference , which vanishes at , rapidly increases for rising temperature, and for it takes a value K at temperatures higher than 1000 K. For larger tensile stress, this difference is smaller, since the amplitude of out-of-plane vibrations is reduced, and is closer to . Thus, for and eV Å, we find at high temperatures K and K, respectively.

The results presented here for at zero external stress display a temperature dependence similar to those found earlier using other theoretical and experimental techniques.Jiang et al. (2009); da Silva et al. (2014); Bao et al. (2009); Yoon et al. (2011) The solid line in Fig. 6(b) represents the thermal expansion coefficient obtained by Mounet and MarzariMounet and Marzari (2005) from DFT combined with a QHA for the vibrational modes. Similar curves were obtained from first-principles calculations in Refs. Sevik, 2014; Mann et al., 2017, converging to negative values at high temperature. This seems to be a drawback of this kind of calculations, which are optimal to obtain the total energy and electronic structure of the system, but the employed QHA may be not precise enough to capture the important coupling between in-plane and out-of-plane vibrational modes at relatively high temperatures. This mode coupling controls the in-plane thermal expansion.

Jiang et al.Jiang et al. (2009) studied free-standing graphene using a nonequilibrium Green’s function approach, and obtained a minimum for of about K, similar to our data for displayed in Fig. 6(b). Moreover, these authors found a crossover from negative to positive for a temperature K, lower than the results of our PIMD simulations. An even lower temperature K has been found by using other theoretical techniques.da Silva et al. (2014); Michel et al. (2015) Room-temperature Raman measurementsYoon et al. (2011) yielded K, whereas a value of K was derived from scanning electron microscopy.Bao et al. (2009) We note, however, that the agreement between measurements with different techniques is not so good for the temperature dependence of , since the change in at 300 K is faster in the former case Yoon et al. (2011) than in the latter.Bao et al. (2009)

The temperature dependence of can be qualitatively understood as a competition between two opposing factors. On one side, the real area increases as is raised in the whole temperature range considered here. On the other side, bending of the graphene surface causes a decrease in its projection, , onto the plane. For stress-free graphene at temperatures 1000 K, the decrease due to out-of-plane vibrations is larger than the thermal expansion of the real surface, and . For 1000 K, the thermal expansion of dominates over the contraction of associated to out-of-plane atomic motion, and thus . For graphene under tensile stress, the amplitude of out-of-plane vibrations is smaller, so the decrease in is also smaller and the region of negative is reduced.

Our results for the thermal expansion of graphene, derived from atomistic simulations, can be related with an analytical formulation of crystalline membranes in the continuum limit, for which the relation between and may be written asWaheed and Edholm (2009); Ramírez and Herrero (2017)

(9) |

Here is the height of the membrane surface, i.e. the distance to the reference plane. In this classical approach, the difference may be calculated by Fourier transformation of the r.h.s. of Eq. (9).Safran (1994); Chacón et al. (2015); Ramírez and Herrero (2017) In this procedure one needs to consider a dispersion relation for out-of-plane modes (ZA flexural band), where are 2D wavevectors. The frequency dispersion in this acoustic band is well approximated by the expression , consistent with an atomic description of grapheneRamírez et al. (2016) (; , surface mass density; , effective stress; , bending modulus). Then, for effective stress , which is the case at finite temperatures, even for zero external stress (), one findsRamírez et al. (2016); Ramírez and Herrero (2017)

(10) |

This expression was obtained in the classical limit, so it does not take into account atomic quantum delocalization. Nevertheless, it is expected to be a good approximation to our quantum calculations at relatively high temperature, , with 1000 K the Debye temperature corresponding to out-of-plane vibrations in graphene.Tewary and Yang (2009); Politano et al. (2011)

We note that both and change with temperature and stress, so that it is not straightforward to write down an analytical formula for from a temperature derivative of Eq. (10). For given temperature and stress we can write , being a parameter derived from Eq. (10), which is a good approximation for the difference . The effective stress appears at finite temperatures for zero external stress (), and vanishes for in the classical limit.Ramírez et al. (2016) Introducing the finite-temperature values of and given earlier,Ramírez et al. (2016); Ramírez and Herrero (2017) we find = K, K, and K for , , and eV Å, respectively. These values are close to those found for from our PIMD simulations (see above). However, the difference between both sets of results increases for rising tensile stress. Thus, for eV Å our simulations yield K, larger than the prediction based on Eq. (10), which is based on harmonic vibrations. We observe that the importance of anharmonic effects (including also quantum corrections) is manifested more clearly for large tensile stresses.

## V Specific heat

We have calculated the specific heat of graphene monolayers as a temperature derivative of the enthalpy derived from our PIMD simulations (see Sec. III.D): . One may ask if our simulations can yield reliable results for this thermodynamic variable, mainly because electronic contributions to are not taken into account in our numerical procedure. This is, however, not a problem for the actual precision reached in our calculations, as the electronic part is much less than the phonon contribution, actually considered in our method. The former was estimated in various works, and turns out to be three or four orders of magnitude less than the phonon part in the temperature range considered here.Benedict et al. (1996); Nihira and Iwata (2003); Fong et al. (2013)

In Fig. 7 we display the specific heat of graphene as a function of temperature for (circles) and eV Å (squares), as obtained from a numerical derivative of the enthalpy. For comparison, the classical Dulong-Petit specific heat is shown as a dashed-dotted line (). At room temperature the specific heat of graphene is about three times smaller than the classical limit, and is still appreciably lower than this limit at K. This is in line with the magnitude of the Debye temperature of graphene, which amounts to about 1000 K for out-of-plane modesTewary and Yang (2009); Politano et al. (2011) and K for in-plane vibrations.Tewary and Yang (2009); Pop et al. (2012) In Fig. 7 we also present results for at taken from earlier DFT calculations combined with a QHA: = 900 K from Ref. Mann et al., 2016 (open square), 850 K and 750 K from Ref. Ma et al., 2012 (open circles). These data agree well with the results of our PIMD simulations.

For a comparison with our numerical results for at zero and negative stress, we present a harmonic approximation (HA) for the lattice vibrations. This approximation is expected to be rather accurate at low temperatures, provided that a reliable description for the phonon frequencies is used. A comparison between results of the HA and PIMD simulations yields an estimate of anharmonic effects in the specific heat of graphene, given that both kinds of calculations employ the same interatomic potential.

The HA assumes constant frequencies for the graphene vibrational modes (calculated for the minimum-energy configuration), and does not take into account changes in the areas and with temperature. Thus, it will give us the constant-area specific heat per atom, , which for a cell with atoms is given by

(11) |

where the index ( = 1, …, 6) indicates the six phonon bands of graphene (ZA, ZO, LA, TA, LO, and TO),Mounet and Marzari (2005); Karssemeijer and Fasolino (2011); Wirtz and Rubio (2004) and the sum in runs over wavevectors in the hexagonal Brillouin zone, with discrete points spaced by and .Ramírez et al. (2016) For increasing system size , one has new long-wavelength modes with an effective cut-off , with , and the minimum wavevector is (i.e., ).

The dashed line in Fig. 7 was calculated with the HA using Eq. (11), with the frequencies ( = 1, …, 6) obtained from diagonalization of the dynamical matrix corresponding to the LCBOPII potential. Results of the simulations for follow closely the HA up to about 400 K, and they become progressively higher than the dashed line at higher temperatures, for which anharmonic effects become more important. At room temperature ( K) we obtained eV/K atom from the HA vs a value of eV/K atom derived from PIMD simulations. Both values are a little higher than that found by Ma et al.Ma et al. (2012) from DFT calculations ( J/K mol, i.e. eV/K atom). The difference between results obtained for from ab-initio calculations combined with a QHA and those presented here for a simpler HA are due to two main reasons. The first reason is that in the HA one does not take into account changes of frequencies with the temperature, and the second is the relative inaccuracy of the phonon bands derived from the considered effective potential far from the point, as compared with those obtained from DFT calculations. The HA provides, however, a consistency check for the results of the PIMD simulations at low temperature, as discussed below.

As shown in Sec. III, a part of the internal energy (and the enthalpy) corresponds to the elastic energy , i.e., to the cost of increasing the area of graphene by thermal expansion or an applied tensile stress. The contribution of this energy to the specific heat is given by , which is not included in the HA. This contribution can be calculated from the results of our PIMD simulations displayed in Fig. 3. For , it amounts to 2.3 and eV/K atom at 500 and 1000 K, which means a nonnegligible increase in the specific heat with respect to the pure HA. This increase is especially visible at K, capturing part of the anharmonicity of the system. In fact, at K it accounts for a 45% of the difference between the results derived from PIMD simulations and the HA presented in Fig. 7. The rest of that difference is associated to anharmonicity of the lattice vibrations.

In our present context, the most relevant effect of an applied stress in the phonon bands is a change in the low-frequency region of the ZA modes, for which

(12) |

The zero-stress band calculated for the minimum-energy structure (area ), follows for small : . Thus, for the small- region is dominated by the quadratic term (linear in ) in Eq. (12), and for 1 Å. The specific heat of stressed graphene in the HA has been calculated by taking the frequencies derived from Eq. (12).

To obtain deeper insight into the low-temperature dependence of the specific heat, is shown in Fig. 8 in a logarithmic plot. Symbols indicate results derived from PIMD simulations for (circles), (squares), and eV Å (diamonds). With our present method, we cannot give reliable values of for temperatures lower than 15 K, mainly due to the very large computation times required for Trotter numbers . Dashed lines show derived from the harmonic approximation using Eq. (11) for a large cell size ( atoms). Such a large cannot be reached in our quantum simulations, and gives us a useful reference where finite-size effects are minimized. The lines corresponding to the HA follow closely the data points derived from the simulations at low temperatures. For eV Å one observes that the PIMD results become progressively higher than the HA line at K. Such a difference is almost unobservable for eV Å, and disappears for . It is clear that the HA describes well the specific heat in the stress-free case at K, and becomes less accurate as tensile stress and/or temperature increase.

The low-temperature behavior of the heat capacity can be further analyzed by considering a continuous model for frequencies and wavevectors, as in the well-known Debye model for solids.Ashcroft and Mermin (1976) At low-temperatures, is mainly controlled by the contribution of acoustic modes with small . For graphene, these are TA and LA modes with , and ZA flexural modes with quadratic dispersion () for negligible at and low , whereas it is linear in for . In general, the low- contribution of a phonon branch with dispersion relation can be approximated by replacing the sum in Eq. (11) by an integral, which yields (see Ref. Herrero and Ramírez, 2018). We note that for -dimensional systems one finds an exponent .Hone (2001); Popov (2002)

Summarizing the above comments, the low-temperature contributions of the relevant phonon branches to the specific heat (those with for ), one has for graphene close to , with for and in the presence of a tensile stress (). This is in fact what we find for from the HA at low temperature in our calculations for a large graphene supercell. To see the convergence of to its low- value, we have calculated this exponent as a function of temperature from a logarithmic derivative: . Close to we obtain the values expected from our discussion above, but in the region displayed in Fig. 8, has not yet reached the zero-temperature value. In fact, for K we find for the values 1.05, 1.59, and 1.85, for , , and eV Å, respectively. The first of them is already close to its low-temperature limit (), but the exponents for graphene under stress are still somewhat lower than their zero- limit (). We note that the larger the tensile stress (larger ), the closer is to 2 at a relatively low , since the wavenumber region where the linear term dominates in the dispersion of the ZA band is larger (). For the results of our PIMD simulations we note that, although one can obtain an exponent from a linear fit in the logarithmic plot of Fig. 8 (in particular the fit is rather good for K), it is true that the actual slope is slowly changing in this temperature range.

Alofi and SrivastavaAlofi and Srivastava (2014) calculated the specific heat of few-layer graphene by using a semicontinuum model and analytical expressions for phonon dispersion relations. In particular, for single-layer graphene they found a temperature dependence for K, with an exponent that coincides with the mean value derived from our results for in the region from 10 to 50 K.

For comparison with our results for graphene, we also present in Fig. 8 experimental data for the specific heat of graphite, obtained by Desorbo and Tyler.Desorbo and Tyler (1953) For graphite, the dependence of on temperature has been analyzed in detail over the years.Komatsu and Nagamiya (1951); Krumhansl and Brooks (1953); Klemens (1953); Nicklow et al. (1972) For K, rises as (a temperature region not shown in Fig. 8), as in the Debye model, and for temperatures between 10 and 100 K, it increases as (). The major difference with stress-free graphene is that in graphite the dominant contribution to in this temperature range arises from phonons with a linear dispersion relation for small (). At room temperature the experimental specific heat of graphite is eV/K atom (8.59 J/K mol), somewhat less than the result for graphene derived from our PIMD simulations for : eV/K atom.

To compare with the results for of graphene derived from our PIMD simulations, we have also calculated the specific heat from constant- simulations. For each considered temperature, we took the equilibrium area obtained in the isothermal-isobaric simulations, and calculated as a numerical derivative of the internal energy

(13) |

from temperature increments (both positive and negative). One expects that at any temperature, but the difference between them turns out to be less than the statistical error bars of our results. A realistic estimation of this difference can be obtained from the thermodynamic relationLandau and Lifshitz (1980); Herrero and Ramírez (2018)

(14) |

where is the in-plane isothermal compressibility (see Sec. VI). Using Eq. (14), we find that is at least two orders of magnitude less than the values given above for the tensile stresses and temperatures considered here.

We finally note that a calculation of the low- specific heat of solids using this kind of path-integral simulations is in general not straightforward. The verification of the Debye law has been a challenge for path-integral simulations of 3D materials, because of the effective low-frequency cut-off corresponding to finite simulation cells.Noya et al. (1996); Ramírez et al. (2006) This kind of calculations at relatively low temperatures are in principle more reliable for 2D materials due to two different reasons. The first reason is that the length of the cell sides scales as , and the minimum wavevector available in the simulations is . For a simulation cell including atoms, is less for 2D than for 3D materials, so that the low-frequency region is represented better for , and therefore the low-temperature regime will be also better described. The second reason is that the internal energy (or enthalpy) for graphene rises as with or 3. At low , this increase is faster than the usual phonon contribution to the enthalpy in 3D materials (), and consequently it is more readily observable (less relative statistical noise) for 2D materials such as graphene.

## Vi Compressibility

The in-plane isothermal compressibility is defined as

(15) |

The variables in the r.h.s. of this equation refer to in-plane quantities, as the pressure in our isothermal-isobaric ensemble is a variable conjugate to the in-plane area . has been calculated here by using the fluctuation formulaLandau and Lifshitz (1980); Ramírez and Herrero (2017)

(16) |

where are the mean-square fluctuations of the area obtained in the simulations. This expression is more convenient for our purposes than obtaining , since a calculation of this derivative by numerical methods requires additional simulations at nonzero stresses. We have verified at some selected temperatures and pressures that both procedures give the same results for (taking into account the statistical error bars). Similarly, for the real area one can define a compressibility , which will be related to the fluctuations of the real area .

In Fig. 9 we present the stress dependence of the compressibilities (circles) and (squares) of graphene, as derived from our PIMD simulations at 300 K. At , one finds , as a consequence of the larger fluctuations in the in-plane area . In this figure we have included some points corresponding to (small) compressive stresses , to remark the very different behavior of and in this region. follows a regular dependence, in the sense that it displays a smooth change as is varied in the considered stress region. The in-plane compressibility , however, increases fast for , which is consistent with an eventual divergence at a critical stress . This divergence of shows the same fact as the vanishing of the in-plane bulk modulus (the inverse of ) discussed in Ref. Ramírez and Herrero, 2017 from classical calculations. Close to this critical compressive stress the planar morphology of graphene becomes unstable due to large fluctuations in the in-plane area, which is related to the onset of imaginary frequencies in the ZA phonon bands for . A detailed characterization of this wrinkling transition at is out of the scope of this paper and will be investigated further in the near future. It is not yet clear whether quantum effects may be relevant or not for a precise definition of . It is remarkable that the compressibility derived from the real area does not display any abrupt change in the region close to the critical stress, as both the area and its fluctuations are rather insensitive to the corrugation of the graphene surface imposed by compressive stresses. For large tensile stress, and converge one to the other, as shown in Fig. 9, since the amplitude of out-of-plane vibrations becomes smaller and the real graphene surface is closer to the reference plane.

For comparison with the results of our simulations, we also present in Fig. 9 data for the compressibility of graphene at derived from atomic force microscopy (AFM) indentation experiments. We have transformed the values for the Young modulus given in Refs. Lee et al., 2008; Lopez-Polin et al., 2015 to compressibility by using the expression , with a Poisson ratio .Ramírez and Herrero (2017) The resulting compressibilities are plotted as open symbols: a squareLopez-Polin et al. (2015) and a circleLee et al. (2008) (the error bar indicates one standard deviation for several measurements). We note that interferometric profilometry experimentsNicholl et al. (2015) have revealed that close to much smaller values of the Young modulus (much larger compressibilities) can be found for graphene. These authors have suggested that this apparent discrepancy can be associated to the difference between compressibilities of the real and projected areas, as discussed here. It seems that some experimental techniques can measure one of them, while other techniques may be sensitive to the other. This is an ongoing discussion that should be elucidated in the near future.Ramírez and Herrero (2017); Nicholl et al. (2017)

A thermodynamic parameter related to the thermal expansion and compressibility is the dimensionless Grüneisen parameter , defined as

(17) |

where is the contribution of mode to the specific heat, and are mode-dependent Grüneisen parameters.Ashcroft and Mermin (1976); Mounet and Marzari (2005) The overall can be related with our in-plane variables in graphene by using the thermodynamic relationAshcroft and Mermin (1976); Herrero and Ramírez (2018)

(18) |

For , the results of our PIMD simulations yield at 300 K, and at 1000 K, within the precision of the numerical results. At K, and are negative because the mode-dependent Grünesien parameter of the out-of-plane ZA modes is negative.Balandin (2011); Mounet and Marzari (2005); Mann et al. (2016) However, for K, this negative contribution is dominated by the positive sign of for in-plane modes, which are excited at these temperatures, so that and are positive.

At room temperature ( K) we find and , for and eV Å, respectively. Although these values of are still negative, they are very different from the zero-stress result (). Looking at Eq. (18), the main reason for this important change in the Grünesien parameter with tensile stress is the large variation in the in-plane thermal expansion coefficient , which takes values of , , and K for , , and eV Å, respectively (see Fig. 6(b)). From the point of view of the phonon contributions to , the small negative value found for eV Å indicates that the relative contribution of phonons with negative (out-of-plane ZA modes) is at room temperature less important than for . This is due to an increase in frequency of small- modes (), which causes a reduction in its corresponding at 300 K (see Eq. (17)).

## Vii Summary

In this paper we have presented and discussed results of PIMD simulations of graphene in a wide range of temperatures and tensile stresses. This technique has allowed us to quantify several structural and thermodynamic properties, with particular emphasis on the thermal expansion and specific heat. Nuclear quantum effects are clearly appreciable in these variables, even at higher than room temperature. Zero-point expansion of the graphene layer due to nuclear quantum motion is not negligible, and amounts to about 1% of the area . This zero-point effect decreases as tensile stress is increased.

The thermal contraction of graphene discussed in the literature turns out to be a reduction of the in-plane area (), caused by out-of-plane vibrations, but not a decrease in the real area . In fact, the difference rises as temperature increases, since the amplitude of those vibrations grows. On one side, the in-plane thermal expansion is negative at low temperature, and becomes positive for 1000 K. On the other side, the thermal expansion of the real area is positive for all temperatures and tensile stresses discussed here. The thermal contraction of is smaller as the tensile stress increases, due to a reduction in the amplitude of out-of-plane vibrations. Our PIMD simulations give at low for stresses so high as eV Å ( N/m). For this stress, becomes positive at K.

The anharmonicity of the vibrational modes is clearly noticeable in the behavior of and . The increase in real area for rising () is a consequence of anharmonicity of in-plane modes, similar to the thermal expansion in most 3D solids. Moreover, the peculiar dependence of (negative at low and positive at high ) is an indication of the coupling between in-plane and out-of-plane modes.

Other thermal properties of graphene can be well described by an HA at relatively low , once the frequencies of the vibrational modes are known for the classical equilibrium geometry at . This is the case of the specific heat, for which our PIMD results indicate that the effect of anharmonicity appears gradually for temperatures 400 K. Part of this anharmonicity is due to the elastic energy , associated to the expansion of the actual graphene sheet, which is not taken into account by the HA. At low the specific heat can be described as . At the lowest temperatures studied here ( 10 K) we find an exponent = 1.05 and 1.85, for and eV Å, respectively. These values are close to the corresponding low-temperature () values, i.e., = 1 for stress-free and = 2 for stressed graphene.

We note the consistency of the simulation results with the principles of thermodynamics, in particular with the third law. This means that thermal expansion coefficients have to vanish for , as found for both the in-plane area and the real area derived from our simulations for and . The same happens for the specific heat, i.e., for .

###### Acknowledgements.

The authors acknowledge the help 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.## References

- Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nature Mater. 6, 183 (2007).
- 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).
- 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).
- Balandin (2011) A. A. Balandin, Nature Mater. 10, 569 (2011).
- 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).
- Lee et al. (2008) C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008).
- Prasher (2010) R. Prasher, Science 328, 185 (2010).
- 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).
- 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).
- Fasolino et al. (2007) A. Fasolino, J. H. Los, and M. I. Katsnelson, Nature Mater. 6, 858 (2007).
- Herrero and Ramírez (2016) C. P. Herrero and R. Ramírez, J. Chem. Phys. 145, 224701 (2016).
- 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).
- 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).
- 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).
- Helfrich and Servuss (1984) W. Helfrich and R. M. Servuss, Nuovo Cimento D 3, 137 (1984).
- Fournier and Barbetta (2008) J.-B. Fournier and C. Barbetta, Phys. Rev. Lett. 100, 078103 (2008).
- Tarazona et al. (2013) P. Tarazona, E. Chacón, and F. Bresme, J. Chem. Phys. 139, 094902 (2013).
- Chacón et al. (2015) E. Chacón, P. Tarazona, and F. Bresme, J. Chem. Phys. 143, 034706 (2015).
- Pozzo et al. (2011) M. Pozzo, D. Alfè, P. Lacovig, P. Hofmann, S. Lizzit, and A. Baraldi, Phys. Rev. Lett. 106, 135501 (2011).
- 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).
- Ramírez et al. (2016) R. Ramírez, E. Chacón, and C. P. Herrero, Phys. Rev. B 93, 235419 (2016).
- Nelson et al. (2004) D. Nelson, T. Piran, and S. Weinberg, Statistical Mechanics of Membranes and Surfaces (World Scientific, London, 2004).
- Amorim et al. (2014) B. Amorim, R. Roldan, E. Cappelluti, A. Fasolino, F. Guinea, and M. I. Katsnelson, Phys. Rev. B 89, 224307 (2014).
- Pop et al. (2012) E. Pop, V. Varshney, and A. K. Roy, MRS Bull. 37, 1273 (2012).
- Wang et al. (2016) P. Wang, W. Gao, and R. Huang, J. Appl. Phys. 119, 074305 (2016).
- Fong et al. (2013) K. C. Fong, E. E. Wollman, H. Ravi, W. Chen, A. A. Clerk, M. D. Shaw, H. G. Leduc, and K. C. Schwab, Phys. Rev. X 3, 041008 (2013).
- Mann et al. (2016) S. Mann, P. Rani, R. Kumar, G. S. Dubey, and V. K. Jindal, RSC Adv. 6, 12158 (2016).
- Ma et al. (2012) F. Ma, H. B. Zheng, Y. J. Sun, D. Yang, K. W. Xu, and P. K. Chu, Appl. Phys. Lett. 101, 111904 (2012).
- Alofi and Srivastava (2013) A. Alofi and G. P. Srivastava, Phys. Rev. B 87, 115421 (2013).
- da Silva et al. (2014) A. L. C. da Silva, L. Candido, J. N. Teixeira Rabelo, G. Q. Hai, and F. M. Peeters, EPL 107, 56004 (2014).
- Burmistrov et al. (2016) I. S. Burmistrov, I. V. Gornyi, V. Y. Kachorovskii, M. I. Katsnelson, and A. D. Mirlin, Phys. Rev. B 94, 195430 (2016).
- Mounet and Marzari (2005) N. Mounet and N. Marzari, Phys. Rev. B 71, 205214 (2005).
- Sevik (2014) C. Sevik, Phys. Rev. B 89, 035422 (2014).
- Mann et al. (2017) S. Mann, R. Kumar, and V. K. Jindal, RSC Adv. 7, 22378 (2017).
- Shimojo et al. (2008) F. Shimojo, R. K. Kalia, A. Nakano, and P. Vashishta, Phys. Rev. B 77, 085103 (2008).
- de Andres et al. (2012) P. L. de Andres, F. Guinea, and M. I. Katsnelson, Phys. Rev. B 86, 245409 (2012).
- 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).
- 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).
- Feynman (1972) R. P. Feynman, Statistical Mechanics (Addison-Wesley, New York, 1972).
- Herrero and Ramírez (2014) C. P. Herrero and R. Ramírez, J. Phys.: Condens. Matter 26, 233201 (2014).
- Cazorla and Boronat (2017) C. Cazorla and J. Boronat, Rev. Mod. Phys. 89, 035003 (2017).
- 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. (2005) L. M. Ghiringhelli, J. H. Los, E. J. Meijer, A. Fasolino, and D. Frenkel, Phys. Rev. Lett. 94, 145701 (2005).
- Zakharchenko et al. (2009) K. V. Zakharchenko, M. I. Katsnelson, and A. Fasolino, Phys. Rev. Lett. 102, 046808 (2009).
- Politano et al. (2012) A. Politano, A. R. Marino, D. Campi, D. Farías, R. Miranda, and G. Chiarello, Carbon 50, 4903 (2012).
- Herrero and Ramírez (2018) C. P. Herrero and R. Ramírez, J. Chem. Phys. 148, 102302 (2018).
- Lambin (2014) P. Lambin, Appl. Sci. 4, 282 (2014).
- Shiba et al. (2016) H. Shiba, H. Noguchi, and J.-B. Fournier, Soft Matter 12, 2373 (2016).
- 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.
- Tuckerman et al. (1993) M. E. Tuckerman, B. J. Berne, G. J. Martyna, and M. L. Klein, J. Chem. Phys. 99, 2796 (1993).
- Tuckerman et al. (1992) M. E. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys. 97, 1990 (1992).
- Martyna et al. (1999) G. J. Martyna, A. Hughes, and M. E. Tuckerman, J. Chem. Phys. 110, 3275 (1999).
- Martyna et al. (1996) G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. Klein, Mol. Phys. 87, 1117 (1996).
- 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).
- Waheed and Edholm (2009) Q. Waheed and O. Edholm, Biophys. J. 97, 2754 (2009).
- Gao and Huang (2014) W. Gao and R. Huang, J. Mech. Phys. Solids 66, 42 (2014).
- Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders College, Philadelphia, 1976).
- 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).
- Adamyan et al. (2016) V. M. Adamyan, V. N. Bondarev, and V. V. Zavalniuk, Phys. Lett. A 380, 3732 (2016).
- Herrero and Ramírez (2000) C. P. Herrero and R. Ramírez, Phys. Rev. B 63, 024103 (2000).
- Jiang et al. (2009) J.-W. Jiang, J.-S. Wang, and B. Li, Phys. Rev. B 80, 205429 (2009).
- Bao et al. (2009) W. Bao, F. Miao, Z. Chen, H. Zhang, W. Jang, C. Dames, and C. N. Lau, Nature Nanotech. 4, 562 (2009).
- Yoon et al. (2011) D. Yoon, Y.-W. Son, and H. Cheong, Nano Lett. 11, 3227 (2011).
- Michel et al. (2015) K. H. Michel, S. Costamagna, and F. M. Peeters, Phys. Status Solidi B 252, 2433 (2015).
- Tewary and Yang (2009) V. K. Tewary and B. Yang, Phys. Rev. B 79, 125416 (2009).
- Politano et al. (2011) A. Politano, B. Borca, M. Minniti, J. J. Hinarejos, A. L. Vazquez de Parga, D. Farias, and R. Miranda, Phys. Rev. B 84, 035450 (2011).
- Benedict et al. (1996) L. X. Benedict, S. G. Louie, and M. L. Cohen, Solid State Commun. 100, 177 (1996).
- Nihira and Iwata (2003) T. Nihira and T. Iwata, Phys. Rev. B 68, 134305 (2003).
- Hone (2001) J. Hone, in Carbon nanotubes: Synthesis, structure, properties, and applications, edited by M. S. Dresselhaus, G. Dresselhaus, and P. H. Avouris (Springer, 2001), vol. 80 of Topics in Applied Physics, pp. 273–286.
- Popov (2002) V. N. Popov, Phys. Rev. B 66, 153408 (2002).
- Alofi and Srivastava (2014) A. Alofi and G. P. Srivastava, Appl. Phys. Lett. 104, 031903 (2014).
- Desorbo and Tyler (1953) W. Desorbo and W. W. Tyler, J. Chem. Phys. 21, 1660 (1953).
- Komatsu and Nagamiya (1951) K. Komatsu and T. Nagamiya, J. Phys. Soc. Japan 6, 438 (1951).
- Krumhansl and Brooks (1953) J. Krumhansl and H. Brooks, J. Chem. Phys. 21, 1663 (1953).
- Klemens (1953) P. G. Klemens, Austr. J. Phys. 6, 405 (1953).
- Nicklow et al. (1972) R. Nicklow, N. Wakabayashi, and H. G. Smith, Phys. Rev. B 5, 4951 (1972).
- Landau and Lifshitz (1980) L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, Oxford, 1980), 3rd ed.
- Noya et al. (1996) J. C. Noya, C. P. Herrero, and R. Ramírez, Phys. Rev. B 53, 9869 (1996).
- Ramírez et al. (2006) R. Ramírez, C. P. Herrero, and E. R. Hernández, Phys. Rev. B 73, 245202 (2006).
- Lopez-Polin et al. (2015) G. Lopez-Polin, C. Gomez-Navarro, V. Parente, F. Guinea, M. I. Katsnelson, F. Perez-Murano, and J. Gomez-Herrero, Nature Phys. 11, 26 (2015).