Strong anharmonicity in pristine graphene

The thermodynamic coefficients of a free standing infinite graphene monolayer are calculated using the quasi-classical unsymmetrized self-consistent field method (USF). The basic nonlinear integral equations of this theory are solved numerically in the strong anharmonic approximation. The isothermal and adiabatic elastic bulk moduli, the isochoric and isobaric heat capacities, the thermal expansion, thermal pressure coefficient, and the macroscopic Gruneisen parameter are calculated in terms of the derivatives of a specifically chosen interatomic potential function for different values of stress and for temperatures ranging from below room temperature up to the point of loss of thermodynamic stability. The nearest-neighbor distances vary from ≳1.4 Å to ≲1.8 Å for zero stress. Under stress, these distances decrease. At room temperature the molar heat capacities are ∼5.0 Jmol−1K−1. The elasticity moduli vary from 15.0 eV Å − 2 up to zero at the temperature of loss of stability and are increased by stress. The thermal expansion coefficient has a strong dependence on the temperature and is negative for temperatures lower than ∼340 K. For high temperatures it monotonically increases and decreases with stress. The macroscopic Gruneisen parameter has a strong nonlinear dependence with temperature and is estimated in about 3.0 ÷ 3.7 ; at ∼340 K its value decreases to ∼1.0 K and for even lower temperature it shows a peak and deep structure similar to what has been earlier reported for fullerene C60.


Introduction
The discovery of graphene [1], a truly two-dimensional crystal structure formed of a single layer of carbon atoms, and the discussion about the possibility of the syntesis of other two-dimensional materials as is the case of silicene, germanene [2], and more recently the black phosphorus [3] and boron nitride [4], draw back the attention to the analytical and semi-analytical methods in the theory of crystalline solids. The graphene monolayer is the building block for other allotropes of carbon such as the fullerenes, the nanotubes, the nanowires, and the graphite [5]. It is argued that these two-dimensional crystals must have very unusual physical properties if compared to the three-dimensional crystals. For instance, some authors speculate that they must be strongly anharmonic [6][7][8][9]. For strongly anharmonic crystals, an efficient theory is the unsymmetrized selfconsistent field approximation (USF) [10][11][12][13][14] which is based on the localized particles approach for strongly anharmonic crystals [13]. This approach describes the real state of the particles.
In modern statistical physics, the methods of distribution functions have widespread applications mainly in their classical and quasiclassical variants. They are based on the BBGKY (Bogoliubov, Born, Green, Kirkwood, Yvon) hierarchy of equations [15][16][17] that follows from the Liouville equation [18]. Some simplifying physical assumptions are used which allow to express the two-particle functions in terms of the one-particles ones [10]. Basic microscopic properties can be expressed in terms of these functions. In particular they can be used to describe strong anharmonicity in crystals.
In contrast to another very well known theory of strong anharmonic crystals, the theory of selfconsistent phonons, or pseudoharmonic approximation [19][20][21], the USF deals with the states of the atoms themselves in the unit cells. Therefore, it is more suited to the investigation of the structural, dynamical, and thermal properties of the crystals and particulary to the investigation of defects and surfaces.
The thermodynamic properties of crystals are related to atomic vibrations around their equilibrium positions. The relative displacements are small in comparison with the interatomic distance. Even in quantum crystals, where the zero point vibrations are large, the mean square deviations from the sites are about 30% of interatomic distances while for most of the solids they do not exceed about 10% up to the melting temperatures. This allows for the expansion of the potential energy in a power series of the nuclei displacements in the Born-Oppenheimer approximation [22]. For crystals with strong anharmonicity, perturbative techniques based on harmonic or quasiharmonic approximations do not work. The pseudoharmonic approximation [19][20][21] reduces the study of strong anharmonicity to the analysis of weak interacting phonons and because of that it is more suited to the investigation of the phonon spectrum and related phenomena.
The USF has been applied over the years to several types of crystals but mostly with Bravais lattices, and for non Bravais ionic crystals [23,24]. It has been also applied to noncrystalline solids [25]. Recently [26], we have reported on a generalization of the USF to the case of layered structured crystals. In that work though, the thermodynamics of a free standing graphene monolayer was investigated in the weak anharmonic approximation. The anharmonicity has a close relation to the thermal properties of graphene and on-going debates about the contribution of different phonon modes to heat conduction. The intrinsic thermal conductivity of graphene near the room temperature has been estimated in about ∼2000-5000 Wm −1 K −1 [27][28][29][30]. In the weak anharmonic approximation, the self-consistent potentials are expanded in power series of the anharmonic terms of the potential energy. It is in effect an expansion on temperature and only the terms not higher than the second order are taken into account. The objective of the present work is then twofold. Firstly, we want to extend the method to higher temperatures and to achieve this goal we must put aside the weak anharmonic approximation. Secondly, we want to focus our attention on the behavior of the thermodynamic properties in the range of high temperatures. We choose here as a first application a strictly two-dimensional approach. Therefore, the third dimension is not considered at all. Such a simplified model could nonetheless describe some aspects of such real systems as for instance graphene encapsulated between hexagonal boron nitride (h-BN) [4] monolayers.
2. The USF for structured crystals 2.1. The classical approximation The USF is based on a set of nonlinear integrodifferential equations for the one-particle density matrices and the self-consistent potentials. Let us consider a structured crystal with central pairwise interactions, whose hamiltonian function is Here N c is the number of unit cells, σ the number of atoms in a cell; m μ is the mass and r m  and p m  are the position and momentum of the μ-th atom. The function r r is the energy of interaction between a pair of atoms. Assuming that the distribution function is not symmetrical with respect to the permutation of phase coordinates, it follows from the Liouville equation a hierarchical chain of equations for the partial distributions [10] that can be broken on the first link if the two-particle functions are taken as a product of two one-particles ones. The result is the following closed equations for the one-particle functions is the force exerted on the ith atom in the μth unit cell. The equilibrium solutions of (2) can be sought in the form [10] W r p m e w r , 2 , 4 is the Boltzmann constant and T the absolute temperature). Here, the spatial probability densities w r i i m m  ( )must obey a system of σ N c nonlinear integral equations. For a polyatomic, free of defects crystal, the spatial densities of all the atoms belonging to the same sub-lattice have the same form but are shifted at the lattice vectors w r w q are the displacements of the atom from its lattice site, An i  is the position of the unit cell, and R m  is the position of the lattice site inside that unit cell (Â is the lattice matrix [31] and n i  is a vector with any integer components (For a crystal with defects and surfaces these components assume noninteger values). Usually, the origin of coordinates is chosen at one of the lattice sites so that R 0 1 =


. Taking all this into account and also assuming that our crystal is formed of electrically neutral structural elements (atoms or molecules), we can write these equations in the following form of N c identical systems of σ nonlinear integral equations w q K q q w q dq ln 1 0. 5 Here, the constants of integration Λ μ are given by the normalization conditions The sum in the second term of the left-hand side of (5) is the one-atom self-consistent potential u q m  ( ). The kernels Here, d is the dimension of space. This form of the free energy allows for the calculation of quantum corrections [32,33].

The quasi-classical approximation
The diagonal elements of the one-particle equilibrium density matrix for a perfect structured crystal are given in the coordinate representation by Here, β=1/θ; q x m  ( ) is the first non-trivial term of the ÿ-expansion for the diagonal elements of the unnormalized one-particle density matrix. It contains quantum corrections both to the one-particle probability densities and to the classical self-consistent potentials. The bar over the variable denotes averaging with respect to the classical distribution function w q The functions ξ μ obey the equations [34] and also because of the the second term in the definition of the self-consistent potential u q m  ( ), the free energy recovers the usual formula (7) with a first quantum correction additional term [14] by

Solution of the basic equations
The expression for the self-consistent potential that obeys the Bogoliubov's variational principle [14] is There are at least three ways of solving the nonlinear equations (12). Either using an iterative schem, or a variational method, or else expanding the self-consistent potentials in a power series of atomic displacements from the lattice sites. We choose here the third way because it allows us to construct a zeroth approximation with the optimal inclusion of anharmonicity. Usually, it is enough to take into account terms up to the fourth order (but for crystals with defects and surfaces higher terms must be taken into account in the zeroth approximation). The above indicated expansion gives Here, the coefficients are given by The equations (13), (14), and (16) constitute a set of transcendental equations with respect to the moments (16) of the one-particle phase probability configurational density functions. To calculate the coefficients K m n a b mn , we express the atomic force constants m n F a b mn in terms of the total derivatives of the interatomic potential function Φ(R) [35]. Using the operator D R d dR we obtain the following formulae Here, for the case of a perfect crystal, The equation of state can be deduced from (7) with the aid of the solution of these equations. In the following, we shall apply this theory to an infinite two-dimensional graphene monolayer.

The thermodynamics of an infinite graphene monolayer
A monolayer of graphene has a two-dimensional honeycomb crystalline lattice with two carbon atoms per unit cell (σ=2). In figure 1 it is shown a fragment of this lattice. The components of the lattice vectors put on columns form the lattice matrix Â [31] of the lattice.
We restrict ourselves here to the case of nearest-neighbor interactions only. Then considering that for a homogeneous membrane r r F =F mn ( ) ( ), we obtain equal coefficients (14) for all the atoms and can supress the index μ K a a 3 , ) are the lattice vectors, and a is the nearestneighbors distance. Shown also the unit cell (in yellow) and the first (red) and second (blue) coordination circumferences.
From symmetrical considerations, we have q q q q .
x y x y n Now, let us assume that there is only planar displacements and use polar coordinates, q q q q cos , sin The role of the out-of-plane coordinates deserves an special attention and is beyond the scope of this work. Nevertheless, we notice here that the out-of-plane modes should give an essential contribution to the structural, dynamic and thermodynamic properties especially at higher temperatures). Then, owing to the symmetry, we can reduce the self-consistent potentials to the following form Thus, we reduce the set of transcendental equations to only two identical ones In [26], this equation was solved analytically in the weak anharmonic approximation. Here, though, we shall take into account strong anharmonicity and that will lead us to a numerical solution. For this two-dimensional crystal, the two identical equations (21) can be reduced to the following form [36] q K , 2 This function has the asymptotic behavior Here, D ν (Z) is the parabolic-cylinder function [37], and X an adimensional argument. Thus, we have for the specific free energy where N is the number of atoms of the membrane. The first and second terms, given by are respectively the zeroth approximation that carry the contribution of the main anharmonic terms, and the first quantum correction whose account expands the range of application of the method to the low temperature. Early estimates for ordinary three-dimensional crystals based on the agreement with experimental measurements have shown that this method works well down to one third of the Debye temperature. [ derivatives of the force coefficients. We also have to calculate the derivatives of the function β 2 (X) with respect to temperature and to the atomic separation, which are expressed in terms of the derivatives of the adimensional argument X It is convenient to express the total derivatives of X 2 b ( ) in terms of the derivatives of the function ζ 2 (X)=β 2 (X)/X. It follows that X X X a ; 4 2 3 We first calculate the equation of state, i.e. an equation for pressure, P P P Q 0 = + as a function of the atomic separation and of the temperature P=P(a, T). We call attentio to the reader that in the two-dimensional case, P stands for the isotropic stress in the membrane. We obtain the following expressions Here, P is the stress, v a 2 2 k = is the area of the unit cell, and 3 3 2 k = is the form factor for the graphene. The specific internal energy, e=e 0 +e Q , is obtained by means of the well known Gibbs-Helmholtz equation [40]. The result is e K a 1 2 The isothermal elastic modulus is calculated as We obtain for these quantities We calculate the thermal pressure coefficient as We find the following expressions The thermal expansion coefficient, the constant pressure specific heat, the adiabatic elastic modulus, and the macroscopic Gruneisen parameter are now respectively calculated using the formulae

Numerical results and discussion
In order to perform a numerical analysis of our results, we have to choose an interatomic potential function among the several types that have been proposed over the last years for different materials [41]. An empirical interatomic potential has been specially constructed for graphene [42] based on the well known Tersoff-Brenner potential [43,44]. This potential has the same form of the Tersoff-Brenner potential function but with a modification to impose the account of angular forces. It also introduces explicitly the dependence of the potential function on the z components of the atomic displacements and the non-nearest neighbor (long range) interactions. It is convenient to make the calculations in reduced units choosing the well depth ò and the minimum neighbor distance r 0 as the units of energy and distance respectively, and to express the quantum corrections through the de Boer parameter, r m 2 0   p L = [45]. We begin our analysis by numerically solving the equation of state P=P(a, T) for some given pressure. As a result, we obtain the temperature dependence of the nearest neighbor distance. In figure 2 we show some isobars for zero (top panel) and nonzero pressures (bottom panel) in the quasiharmonic (traced lines) and in the anharmonic (solid lines) approximations. The solutions have two branches. The lower branches of the curves correspond to thermodynamically stable states (positive isothermal compressibility), whereas the higher ones to the absolute unstable (negative compressibility); both branches coalesce at the point of loss of thermodynamic stability (nule compressibility). Notice that (top panel), in general, as one increases the temperature, the nearest neighbor distance increases as well, at least for temperatures higher than 500K. In the quasi harmonic approximation the crystal is stable up to about 6000K and the nearest neighbor distance reaches about 1.6 Å, whereas in the anharmonic approximation the temperature of the crystal is shifted to about 12000K and the lattice parameter to 1.75 Å. It means that the inclusion of anharmonicity stabilizes the crystal at high temperatures allowing the atoms to move more around reaching about 10% more than in the weak anharmonic regime in the thermodynamic stability threshold. The insets are meant to magnify the regime of low temperatures where the quantum effects are important enough. For temperatures below about 1250K the quantum corrections induces a different behavior of that of the classical case. The latter has a decreasing almost linear while the former decreases to a minimum around 350K beyond which it increases with increasing of temperature. It is shown also that our results for the anharmonic approximation are in good agreement with those from classical Monte Carlo Simulations (in solid red circles) [6] and in very good agreement also with Figure 2. Temperature dependence of the nearest neighbor distance for zero (top panel) and nonzero (bottom panel) pressures; the insets are meant to magnify the regime of low temperatures. Open blue squares and solid red circles for PIMC [46] and classical Monte Carlo [6]. Note that the x-axis starts at 200K simulations, respectively. recently published results of quantum path integral Monte Carlo (PIMC) simulations (in open blue squares) [46]. One can notice that a combination of quantum and anharmonicity effects changes the behavior of the isobar in the region of temperatures for which other theoretical and experimental approaches indicate a change of sign in the thermal expansion coefficient. It is quite intriguing that even choosing such a simplified model, we have been able to observe such an interesting physics. Of course, the effects related to thermal transverse fluctuations, the so called flexural phonons as well as those related to the formation of ripples [47], must be taken into consideration in order to get more close contact with the experiment. Nevertheless, the good agreement of our simple model results with other experimental and theoretical and simulation data makes us think that the combination of quantum effects with the anharmonicity of the interatomic interactions would somehow imitate the action of these out-of-plane contributions. We would like to reinforce that even neglecting out-of-plane vibrational modes, our results for the flat membrane are a good reference for first estimates because they resemble the physics of a free-standing graphene monolayer in agreement with experiments. Very recently, PIMC simulation results [48] in agreement with [49] have shown that the crossover between classical and quantum contribution for the out-of-plane vibrational modes should occur around 50K. For lower temperatures quantum fluctuations should leave the membrane flat on the large scale. For temperatures higher than that the membrane behaves classically (for the modes out-of-lane) and also would be nearly flat with small wavelength riples. Therefore, our present results are basically in the regime where the membrane is pratically flat and this should be the reason why our results agree quite well with the experiments around the room temperature. In the bottom panel, we analyse the effects of slightly changing the pressure at the regime of anharmonicity. At low temperatures, the effect of pressure is small but as one increases the temperature, the dependence of the nearest-neighbor distance with the pressure turns out to be more strong, and we observe a decreasing with the increasing of pressure. In the quasiharmonic approximation, the temperature of loss of stability is about 5500K. With the account of the high anharmonicity, this value shifts towards about 94% reaching a value around 10700K. One can expect that the account of dynamical correlations by means of classical thermodynamic perturbation theory would considerably reduce this shift. We notice here to avoid misconception that the dynamical instability lays above the melting temperature. Now we discuss the molar heat capacities. They are plotted as a function of temperature for different pressures in figure 3. For low temperatures, both isochoric and isobaric molar heat capacities grow equally sharp with temperature. As seen in the plot, the USF results for these two quantities tend to the low temperature limit given by the quantum generalization of the USF approximation. For high temperatures, the isochoric molar heat capacity remains below the Dulong-Petit value while decreasing with the increasing of temperature, whereas the isobaric molar heat capacity extrapolates the Dulong-Petit value and grows very fast with temperature. The growth of c P and decay of c V at high temperatures is are due strictly to the anharmonicity.
For lower temperatures, one can apply the quasiharmonic approximation. In this case, the quantum generalization for USF is the Einstein model [50], in which the atoms are taken as vibrating independently with the same frequency. In this model, the specific free energy has the following view For a nondeformed crystal all the frequencies are equal to The procedure of differentating the free energy to obtain the thermodynamic coefficients in this case is analogous to what has been done in the preceding section.
In order to show that our results for lower temperatures are reliable, we show in figure 4 the threedimensional case in its comparison with known results for graphite. We plot our three-dimensional low temperature quantum results together with results from experiment for graphite taken from [51], and theoretical calculations made by [52][53][54]. Our results for the 3d case are in good agreement with more recent results [30]. As regarding the low regime of temperature, here we only scratch the surface of the problem. Our main interest in the present work is the regime of high temperatures. Nevertheless, we find it useful to show that our results for high temperature strong anharmonic USF performed fairly well in reproducing the low temperature limit results given by the Einstein quantum generalization which in its turn, as seen in figure 4, are in good agreement with other experimental and theoretical results.
In figure 5 we plot our results for the temperature dependence of the elastic moduli for different values of pressure. As expected, these moduli are equal at sufficiently low temperatures when the anharmonicity is negligible. At higher temperatures, the isothermal modulus is smaller than the adiabatic which is due to the anharmonicity. Therefore, the anharmonicity increases the hardness of the crystal. These moduli tend to zero as we approach the thermodynamic stability threshold. In general, an increase of the volume of the crystal gives rise to a decrease of the moduli. At the vicinity of the point of loss of instability the heat capacity c P  ¥ and the modulus b 0 T  which means a violation of the thermodynamic criteria of stability. To calculate the thermal expansion coefficient, we use the first formula of (39). Thus, we have first to calculate the thermal pressure coefficient Γ. Its dependence on temperature and pressure is plotted in figure 6. Analogously to the case of the elastic moduli, here we can see that the strong anharmonic effects lead to the increase of its dependence with temperature.
The thermal expansion coefficient is plotted in figure 7. The obtained results are in agreement with the weak anharmonic approximation for low temperatures. It has been conjectured that thermal excitations could induce large-scale ripples perpendicular to the graphene layer [41] which could change the electronic properties of the graphene membrane. Due to that, the in-plane thermal expansion coefficient would turn negative at low temperatures. The microscopic mechanisms of the structural deformations associated with this and the temperature dependence of the thermodynamic properties are still not completely understood [41]. There have been several theoretical and experimental attempts to accurately calculate or measure the thermal expansion coefficient and other thermodynamic properties of crystalline graphene, but the results are still very discrepant. In this particular application of our method, the model is flat. We do not take into account the z coordinate and therefore the so called flexural phonons. As one can see though, here we again notice that the quantum    correction combined with the anharmonicity of interatomic interactions in our model is the one responsible for changing the sign of the thermal expansion coefficient. Thus, we are able to see this effect without any more sophisticated model like those studied in several other works [6,41,[55][56][57][58][59].
A quantitative estimate of the overall anharmonicity of the atomic vibrations is given by the macroscopic Gruneisen parameter, which represents an average of the individual Gruneisen parameters that give the dependence of the vibration frequencies with the crystal volume. This parameter establishes a link between the thermal and caloric properties. The assumption often used for solids that this parameter is independent of temperature [60] cannot be applied to the graphene membrane. As noted in the preceding section, we calculate here the macroscopic Gruneisen parameter using the last equation (14). In figure 8, γ is shown as a function of temperature for different values of pressure. Analogously to what happens to fluids [60], the γ decreases with increasing temperature until the region of loss of stability is reached from which it increases sharply.
The strong dependence of the Gruneisen parameter on the temperature and pressure is a clear evidence of the strong anharmonicity of the graphene membrane. As shown in the inset of figure 8, around the temperature for which the thermal expansion coefficient changes its sign the Gruneisen parameter shows a sharp (peak and dip) structure very similar to what has been earlier reported for the Buckminsterfullerene, C 60 [61], at the region of a glassy phase transition. For other temperatures, the Gruneisen paramenter of C 60 was found [61] to be close to that for graphite along the c-axis. We find that in wide range of temperature γ≈3.4 and this estimate is an evidence of very high degree of anharmonicity in the graphene membrane. This estimate is higher than values prior used for graphite [62].

Conclusion
In conclusion, we have extended the USF theory to take into account strong anharmonicity in layered structured crystals. The theory was used here to investigate the thermodynamic properties of a free standing perfect infinite graphene monolayer. The main thermodynamic coefficients have been calculated for several values of pressure in a wide range of temperature up to the temperature of loss of thermodynamic stability which is about 10700K. Although we have considered here an infinite nondefective membrane, we have observed negative thermal expansion coefficient for low temperatures. Our calculations show that quantum effects are important in the thermodynamics of graphene for T<1500 K. The combined effect of anharmonicity and of the quantum correction is important for temperatures at the vicinity of the point where experiment [55] predicts the change in the sign of the thermal expansion coefficient. Here though, we have extended our calculations to the regime of high temperatures up to the spinodal point. It is important that we can include any degree of anharmonicity in the zeroth approximation without perturbative schemes that if calculated would yield the account of dynamical interatomic correlations. One of the effects that one could expect of the account of dynamical correlations is the estimate of a lower value for the temperature of loss of thermodynamic stability. The present method of calculation allows for a free choice of the interatomic potential function. The developed theory in this work can be used to investigate multilayer graphene as well as other layered crystals with in principle any dimension.