Second order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: theory and stochastic implementation
Raffaello Bianco, Ion Errea, Lorenzo Paulatto, Matteo Calandra,, Francesco Mauri
TL;DR
This paper develops an analytic and stochastic method within the self-consistent harmonic approximation to study second-order phase transitions and anharmonic phonons in solids, enabling efficient first-principles calculations of free energy curvature and spectral properties.
Contribution
It derives an exact formula for the free energy curvature in the self-consistent harmonic approximation and introduces a stochastic implementation suitable for first-principles density functional theory.
Findings
Analytic expression for free energy second derivative derived.
Method effectively identifies phase transition instabilities.
Application demonstrated on a model mimicking ferroelectric transition.
Abstract
The self-consistent harmonic approximation is an effective harmonic theory to calculate the free energy of systems with strongly anharmonic atomic vibrations, and its stochastic implementation has proved to be an efficient method to study, from first-principles, the anharmonic properties of solids. The free energy as a function of average atomic positions (centroids) can be used to study quantum or thermal lattice instability. In particular the centroids are order parameters in second-order structural phase transitions such as, e.g., charge-density-waves or ferroelectric instabilities. According to Landau's theory, the knowledge of the second derivative of the free energy (i.e. the curvature) with respect to the centroids in a high-symmetry configuration allows the identification of the phase-transition and of the instability modes. In this work we derive the exact analytic formula for…
| Symbol | Meaning | First use |
|---|---|---|
| Atomic position (canonical variable) | Eq. (1) | |
| Potential energy | Eq. (1) | |
| Centroids position (parameter) | Eq. (6) | |
| SCHA free energy for the centroids | Eq. (14) | |
| Minimum point of | Eq. (30) | |
| Minimum point of | Eq. (15) | |
| Generic trial harmonic matrix (parameter) | Eq. (6) | |
| Probability distribution of for a given value of the parameters , | Eq. (6) | |
| that minimizes the SCHA free energy functional at a given | Eq. (16) | |
| Average of the -th derivative of with the probability | Eq. (20) | |
| Second derivative of in | Eq. (51) | |
| -th derivative of in | Eq. (52) | |
| Second derivative of in , divided by the square root of the masses | Eq. (29) | |
| Matrix divided by the square root of the masses | Eq. (32) | |
| Tensor divided by the square root of the masses | Eq. (40) | |
| Matrix divided by the square root of the masses | Eq. (30) | |
| Tensor divided by the square root of the masses | Eq. (60) | |
| Inverse of the matrix | Eq. (34) | |
| Green function associated to | Eq. (34) | |
| Green function associated to | Eq. (53) |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Second order structural phase transitions, free energy curvature, and temperature-dependent
anharmonic phonons in the self-consistent harmonic approximation: theory and stochastic implementation
Raffaello Bianco1,2
Ion Errea3,4
Lorenzo Paulatto1
Matteo Calandra1
Francesco Mauri2,5
1Institut de minéralogie, de physique des matériaux et de cosmochimie (IMPMC), Université Pierre et Marie Curie (Paris VI), CNRS UMR 7590, IRD UMR 206, Case 115, 4 place Jussieu, 75252 Paris Cedex 05, France
2 Dipartimento di Fisica, Università di Roma La Sapienza, Piazzale Aldo Moro 5, I-00185 Roma, Italy
3Fisika Aplikatua 1 Saila, Bilboko Ingeniaritza Eskola, University of the Basque Country (UPV/EHU), Rafael Moreno “Pitxitxi” Pasealekua 3, 48013 Bilbao, Basque Country, Spain
4Donostia International Physics Center (DIPC), Manuel de Lardizabal pasealekua 4, 20018 Donostia-San Sebastián, Basque Country, Spain
5 Graphene Labs, Fondazione Istituto Italiano di Tecnologia, Via Morego, I-16163 Genova, Italy
Abstract
The self-consistent harmonic approximation is an effective harmonic theory to calculate the free energy of systems with strongly anharmonic atomic vibrations, and its stochastic implementation has proved to be an efficient method to study, from first-principles, the anharmonic properties of solids. The free energy as a function of average atomic positions (centroids) can be used to study quantum or thermal lattice instability. In particular the centroids are order parameters in second-order structural phase transitions such as, e.g., charge-density-waves or ferroelectric instabilities. According to Landau’s theory, the knowledge of the second derivative of the free energy (i.e. the curvature) with respect to the centroids in a high-symmetry configuration allows the identification of the phase-transition and of the instability modes. In this work we derive the exact analytic formula for the second derivative of the free energy in the self-consistent harmonic approximation for a generic atomic configuration. The analytic derivative is expressed in terms of the atomic displacements and forces in a form that can be evaluated by a stochastic technique using importance sampling. Our approach is particularly suitable for applications based on first-principles density-functional-theory calculations, where the forces on atoms can be obtained with a negligible computational effort compared to total energy determination. Finally we propose a dynamical extension of the theory to calculate spectral properties of strongly anharmonic phonons, as probed by inelastic scattering processes. We illustrate our method with a numerical application on a toy model that mimics the ferroelectric transition in rock-salt crystals such as SnTe or GeTe.
I Introduction
Describing accurately atomic vibrations is crucial in many branches of physics and chemistry because thermodynamic, transport, and superconducting properties of materials and molecules as well as the spectra obtained in many spectroscopic techniques depend on how atoms vibrateBorn and Huang (1954). The standard harmonic approximation provides the simplest description of vibrations, which are also present at 0 K due to the quantum zero-point motion. The harmonic approximation is based on the expansion of the Born-Oppenheimer (BO) energy surface to the second order around the ionic equilibrium positions. It predicts well-defined non-interacting quasi-particles (phonons) with infinite lifetime and a temperature-independent energy spectrum. Within this approximation many physical effects cannot be described. For example, finite values of the thermal conductivity and temperature dependent effects, like the thermal expansion in solids, cannot be accounted for at the harmonic level. Therefore, it is of paramount importance to describe accurately the vibrations of atoms beyond the harmonic approximation.
Anharmonic effects, i.e. effects due to higher orders in the energy surface expansion, introduce interaction between phonons, thus finite scattering rates and finite lifetimes. Anharmonicity can be treated at different levels of theory. The basic approach is to consider higher order terms in the potential expansion as a small perturbation of the harmonic potentialMaradudin and Fein (1962). However, the perturbative approach can be used under quite restrictive conditions: the displacements of the atoms must be within the range in which the harmonic approximation is valid so that higher order terms are considerably smaller than the harmonic potential. Unfortunately, there are several cases in which a non-perturbative regime is reached. For example when light atoms such as hydrogen are present Errea et al. (2013, 2014, 2015, 2016); Borinaga et al. (2016); Rousseau and Bergara (2010), or when the system is close to a dynamical instability (a phase transition) as in ferroelectrics or materials undergoing a charge-density wave (CDW) instability Pawley et al. (1966); Delaire et al. (2011); Luspin et al. (1980); Weber et al. (2011); Leroux et al. (2015); Holt et al. (2001); Leroux et al. (2012); Ghosez et al. (1998); Calandra et al. (2009); Calandra and Mauri (2011); Haas (1965); Ong et al. (2001); Di Salvo et al. (1976); Bianco et al. (2015). In these cases, a non-perturbative approach is required in order to account for anharmonic effects Errea (2016).
Anharmonic effects at a non-perturbative level are commonly treated within molecular dynamics (MD) simulations or methods based on them de Koker (2009); Zhang et al. (2014); Wang et al. (1990); Magdau and Ackland (2013); Ljungberg and Íñiguez (2013); Zhang et al. (2014); Hellman et al. (2011, 2013); Hellman and Abrikosov (2013). However, these approaches are computationally expensive as long simulation times are needed to obtain converged renormalized phonon energies and have an intrinsic limitation because they are based on Newtonian dynamics, which limits their application to temperatures above the Debye temperature. This problem can be overcome by path-integral molecular dynamics Ceperley (1995), but the even greater computational cost that the method needs to incorporate the quantum character of atomic vibrations makes it challenging. To surmount these difficulties, several methods Errea et al. (2011); Tadano and Tsuneyuki (2015); Monserrat et al. (2013); Errea et al. (2014, 2013); Brown et al. (2013); Georgescu and Mandelshtam (2012); Patrick et al. (2015) have been developed, mainly inspired by the self-consistent harmonic approximation (SCHA) devised by Hooton Hooton (1955). The main idea of the SCHA is to use a variational principle, the Gibbs-Bogoliubov (GB) principle, in order to approximate the free energy of the true ionic Hamiltonian with the free energy calculated with a trial harmonic density matrix for the same system. In particular, in the stochastic self-consistent harmonic approximation (SSCHA) Errea et al. (2014, 2013), the free energy is explicitly minimized by using a conjugate-gradient algorithm with respect to the independent coefficients of the trial harmonic potential. In the SSCHA the free energy and its gradient are evaluated through averages computed with stochastic sampling of the configuration space and the importance sampling technique Errea et al. (2014, 2013). In that way the (approximated but non-perturbative) anharmonic free energy of the system is directly accessible. The stochastic approach is particularly suited to be used in conjunction with ab initio calculations, and it has been employed to study thermal anharmonic effects in several compounds such as hydrides and transition metal dichalcogenides Errea et al. (2013, 2015, 2014, 2016); Borinaga et al. (2016); Leroux et al. (2015); Borinaga et al. (2016).
Within the SCHA the free energy as a function of the average atomic positions, i.e. the centroids positions, can be estimated for any temperature. These can be used to study structural second order phase transitions like, for example, in some ferroelectric and CDW phase transitions Pawley et al. (1966); Delaire et al. (2011); Luspin et al. (1980); Weber et al. (2011); Leroux et al. (2015); Holt et al. (2001); Leroux et al. (2012); Ghosez et al. (1998); Calandra et al. (2009); Calandra and Mauri (2011); Haas (1965); Ong et al. (2001); Di Salvo et al. (1976); Bianco et al. (2015). In general, at any temperature the system is in equilibrium in the configuration which minimizes the free energy. According to Landau’s theory Landau and Lifshitz (1980), in a second order phase transition the high temperature free energy minimum is in a high-symmetry phase. As the temperature decreases, the minimum becomes less and less pronounced until it becomes a saddle point at the transition temperature , and, on lowering the temperature further, the equilibrium position moves continuously towards lower symmetry configurations, where the free energy is smaller (Cfr. Fig. 1). In this scheme the observable to be studied as a function of temperature is the second derivative of the free energy with respect to the centroids positions, i.e. the free energy curvature, in the high-symmetry phase. Indeed the Hessian in the high-symmetry phase is positive-definite at high-temperature, but lowering the temperature it develops first a null () and then negative () eigendirection, which indicates the instability distortion that lowers the free energy.
Using the SSCHA code Errea et al. (2013, 2014) it is possible to compute the free energy for several centroids positions and, therefore, to calculate numerically, by finite difference, the curvature in a point. This has been done, for example, to study the quantum H-bond symmetrization in the record superconductor H3S Errea et al. (2016). However, even if legitimate, this ‘brute force’ approach to compute the free energy curvature is computationally demanding. In fact, a careful calculation of the curvature by finite differences requires small statistical noise, implying a great deal of calls to the total-energy-force engine used. Moreover, it also requires SSCHA calculations in the low-symmetry distorted phase, which are always more statistically demanding because the number of free parameters in the trial harmonic Hamiltonian is larger due to the reduced symmetry.
Motivated by these considerations, in this paper we derive the general exact analytic expression of the SCHA free energy curvature for a generic atomic configuration. Our approach is similar to the one proposed by Götze and Michel in the context of elastic constants of anharmonic crystals Götze and Michel (1968). We also present an expression of the SCHA free energy curvature that only depends on atomic displacements and forces. The latter is suited for a stochastic implementation in conjunction with any total-energy-force engine. The method we are presenting here allows, thus, to compute the curvature of the SCHA free energy for a given structure straightforwardly once the GB functional has been minimized within the SSCHA. Since the method is much more efficient and precise than any finite-difference approach Errea et al. (2016), it is especially suited to be used in conjunction with first-principles calculations.
Besides the practical achievements, the curvature formula described here has also interesting conceptual consequences. Since the only physical observable given by the SCHA is the (approximated) free energy, the effective SCHA quadratic matrix that minimizes the GB functional must be understood just as an auxiliary quantity, even if its eigenvalues have been often used to calculate renormalized anharmonic phonon spectra Errea et al. (2013, 2015, 2014, 2016); Borinaga et al. (2016); Leroux et al. (2015); Borinaga et al. (2016); Errea et al. (2011). A significant parallelism can be traced with the Hartree-Fock approximation. In that case, the Rayleigh-Ritz functional of the total energy is minimized with trial Slater determinants describing non-interacting fermions. The energy obtained is an approximation of the true energy of the system, but, on the contrary, the corresponding trial non-interacting many-body wave function and related single particle spectrum do not have in general a physical meaning (except that in some aspects, e.g. see Koopmans’ theorem Koopmans (1934)). An analogous situation occurs with density-function theory (DFT) and the corresponding non-interacting electrons and energy bands Perdew et al. (1982). In the same way, the SCHA matrix (divided by the square root of the masses, in order to have the correct dimensions) cannot be considered a generalized dynamical matrix and, therefore, its temperature-dependent eigenstates do not represent phonons renormalized by anharmonicity. On the contrary, the free energy curvature (in the equilibrium position) divided by the square root of the masses defines a proper anharmonic temperature-dependent generalization of the harmonic dynamical matrix, whose eigenstates represent temperature-dependent anharmonic phonons. Indeed, at variance with the SCHA matrix, which is positive-definite by construction, the dynamical matrix based on the free energy curvature can have negative eigenvalues, and a softening in its spectrum corresponds to a system instability. The free-energy based dynamical matrix is a particularly important tool especially when we consider crystalline systems. Indeed in that case, exploiting the lattice periodicity and the Fourier interpolation technique, it allows to find structural instabilities with any modulations in real space by performing calculations only on a coarse grid of the Brillouin zone.
The theory based on the free energy curvature with respect to the centroids position is ‘static’ in the sense that there are no dynamical variables evolving with time. However, here we also propose a minimal ‘dynamic’ extension of the theory that resembles the work by Goldman et al. Goldman et al. (1970), which allows to study, in a full non-perturbative way, the spectral properties of anharmonic phonons, and thus allows to have finite scattering times and linewidths. This makes the theory able to interpret the results of inelastic scattering processes between anharmonic phonons and external incident particles (typically neutrons), as well as allowing the calculation of the thermal conductivity in strongly anharmonic solids where the harmonic approximation breaks down. Despite the proposed dynamic extension being based on an ansatz, it is reasonable because it yields good results in two limits: at the lowest perturbative level it reduces to the standard perturbation theory result and in the static limit it predicts the same instabilities found with the free energy curvature.
The paper is structured as follows. In section II, we present the fundamentals of the SCHA method, we define the SCHA free energy as a function of the centroids position, and we fix the notation used. In section III we show how to analyze structural second order phase transitions through the Hessian of the free energy with respect to the centroids position (i.e. the curvature), and in section IV we give the explicit expression of the free energy curvature. In section V, on the basis of the results obtained, we describe the temperature-dependent free-energy-based generalization of the harmonic dynamical matrix. In section VI, we express the theory developed so far in a diagrammatic way. In section VII, we show how to implement the curvature formula in a stochastic way and, in particular, how to take into account symmetries in order to speed up the statistical convergence. In section VIII, we find the lowest order perturbative limit of the results obtained. In section IX, inspired by the results obtained, we propose the ansatz to obtain a dynamical extension of the theory. Finally, in section X, we perform numerical tests on a toy model based on the ferroelectric transition in SnTe, with the double objective of demonstrating the correctness of our findings and showing the power of the method. In section XI, we summarize our results and draw some final conclusions. The paper is completed with several appendices including the proofs of all the equations given in the manuscript and the details of the toy model.
II Self-consistent harmonic approximation
We consider the quantum atomic free energy of crystal lattices and molecules. For notation clarity, we derive the main results by using a real space formalism for both cases. This means that in the case of periodic crystals, we are actually studying a supercell with periodic boundary conditions. Later, we will explicitly consider a crystalline case for the numerical example and we will take advantage of translational lattice symmetry. The dynamic of atomic degrees of freedom is determined by the quantum Hamiltonian
[TABLE]
where is the total number of atoms, and are the atom and Cartesian component indices, respectively, is the mass of the -th atom, and are the momentum and position operators, respectively, and is the Born-Oppenheimer potential, where the bold letter indicates in component-free notation. In what follows, we will use bold letters in component-free notation also for other observables and higher order tensors with respect to the index. Moreover, in order to simplify the notation, we will use a single composite index to indicate both atom and Cartesian indices. Notice that double index can be used also for the masses by defining .
Fixed the temperature , the free energy of the ionic Hamiltonian is given by the sum of the total energy and the entropic contribution
[TABLE]
where , ‘tr’ is the trace operation on the atom Hilbert space, and
[TABLE]
is the equilibrium density matrix. In systems comprising many interacting particles, calculating can represent a complicated task. Nevertheless, a quantum variational principle for the free energy can be established. By replacing the true density matrix
in Eq. (2) with a generic density matrix
, we can define the functional
[TABLE]
and the Gibbs-Bogoliubov (GB) variational principle Isihara (1968) states that
[TABLE]
Obviously the equality holds when \scalebox{1.1}{\tilde{\rho}}=\scalebox{1.1}{\rho}.
In particular, the SCHA is obtained by restricting the trial density matrix to equilibrium density matrices
[TABLE]
for the same temperature of a general trial harmonic Hamiltonian for the same particles Hooton (1955). The trial harmonic Hamiltonian is parametrized in terms of the vector of dimension and the square positive-definite matrix of order as
[TABLE]
In what follows we consider only trial harmonic potentials that respect the symmetries of the system.
With \Bigl{\langle}\mathds{O}\Bigr{\rangle}_{{\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}}} we indicate the average of an observable with respect to the density matrix \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}:
[TABLE]
In what follows it will be relevant to consider observables that are function of the position only. In that case, the average can be written as
[TABLE]
where \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}(\boldsymbol{R}) is the diagonal part of the density matrix \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}} in coordinate representation Errea et al. (2013, 2014):
[TABLE]
Here ‘det’ indicates the determinant and is the symmetric matrix associated to by
[TABLE]
where and are eigenvalues and corresponding eigenvectors of , and is the bosonic average occupation number associated to . Note that to have a normalizable distribution \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}(\boldsymbol{R}), and thus must be positive-definite matrices, as specified above.
From Eqs. (9)–(10) we see that the average positions of the atoms for the trial density matrix \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}, namely the ‘centroids’, coincide with :
[TABLE]
According to the GB variational principle, the best approximation of the free energy within the SCHA is , given by
[TABLE]
With we indicate the SCHA free energy for the centroids position ,
[TABLE]
and with we indicate the configuration that minimizes :
[TABLE]
Therefore, is the SCHA equilibrium configuration of the centroids at the considered temperature.
Given a configuration for the centroids, we define the corresponding SCHA square matrix as the matrix that minimizes the functional \mathcal{F}[\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}] with respect to . Therefore, from Eq. (14) we have
[TABLE]
The SCHA matrix satisfies the following self-consistent equation (see Eq. (100)):
[TABLE]
Notice that, for clarity, we are using two different symbols for the generic trial matrix (a ‘slanted’ phi), and for the SCHA matrix , the specific trial matrix that minimizes \mathcal{F}[\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}] with respect to for a given centroids position .
In the rest of this paper, we will consider exclusively the SCHA approximation for the free energy. Therefore, in order to simplify the notation, in what follows we can safely omit the superscript without ambiguity: will always refer to the SCHA free energy. Moreover, as guide and reference for the reader, we collect in table 1 some symbols used in the text with a concise description. Several symbols collected in the table will appear later in the course of the paper.
III Structural second order phase transition and curvature of the free energy
In second order phase transitions involving the position of the atoms, e.g. in ferroelectric and in charge-density wave phase transitions Pawley et al. (1966); Delaire et al. (2011); Luspin et al. (1980); Weber et al. (2011); Leroux et al. (2015); Holt et al. (2001); Leroux et al. (2012); Ghosez et al. (1998); Calandra et al. (2009); Calandra and Mauri (2011); Haas (1965); Ong et al. (2001); Di Salvo et al. (1976); Bianco et al. (2015), we can use the centroids to define the order parameter, which is the observable measured in diffraction experiments. The (temperature-dependent) function rules the phase transitions. At each temperature, the system is in equilibrium in the (temperature-dependent) configuration , where has a minimum. Therefore, in the first derivative of is zero, , and the Hessian matrix of (i.e. the curvature), , is positive-definite.
Landau’s theory of second order phase transitions Landau and Lifshitz (1980) shows that above a certain critical temperature the equilibrium configuration is in a high-symmetry phase . As decreases and approaches from above, the minimum of in becomes less and less pronounced. At , becomes a saddle point, i.e. the Hessian of in develops at least one null eigenvalue, which becomes negative by lowering further the temperature. At the same time, the minimum point , now depending on temperature, continuously deviates from to different configurations having lower symmetry. Since during the phase transition the equilibrium configuration remains a continuous function of temperature, these are also called “continuous phase transitions”. In Fig. 1 we show an example of a typical second order phase transition.
In conclusion, at any temperature it is and the phase transition is characterized by the change of character of the Hessian matrix : at it is positive-definite, whereas as it develops at least one negative eigendirection indicating the distortion which decreases the free energy. It follows that a method to estimate the transition temperature and the instability modes of a second order phase transition can be obtained by computing the Hessian of in the high-symmetry configuration and studying its evolution as a function of temperature. In the next section we will find explicit formulas for the first and second derivative of .
IV Derivatives of
From the definition of Eq. (16), we can calculate explicitly the derivatives of with respect to . Here we present only the results, while the derivation is given in appendix A. For the first derivative we have the intuitive result (see Eq. (101))
[TABLE]
The derivative of the free energy is the average of the potential derivative. In other words, the forces on the centroids are equal to the average of the mechanical forces on the atoms. From Eq. (18), for the equilibrium position defined in Eq. (15) it is
[TABLE]
For what follows it is convenient to define the -th order SCHA tensor, which generalizes Eq. (17) to higher orders:
[TABLE]
Notice that we did not use the superscript for the square SCHA matrix. The -th order SCHA tensor has the same properties of the -th order force constant (i.e. the -th derivative of the potential). Notably, it is invariant with respect to permutation of indices; it is invariant with respect to all the symmetry operations (including lattice translations in a crystal) associated to the configuration MARADUDIN and VOSKO (1968); and it satisfies the acoustic sum rule (ASR), i.e. the sum over any atom index vanishes (see Eq. (166)).
Deriving a second time Eq. (18) with respect to we obtain (see Eqs. (102)–(111) )
[TABLE]
where
[TABLE]
Here and are eigenvectors and eigenvalues of , respectively, and
[TABLE]
The tensor is the solution of the Dyson-like equation
[TABLE]
Notice that in all these equations the dependence of the quantities on is understood.
We have obtained for the second derivative a relation which is different from the one found for the first derivative. Indeed, as shown in Eq. (18), the first derivative of the SCHA free energy is equal to the average of the first derivative of the potential. On the contrary, the second derivative of the SCHA free energy is equal to the average of the second derivative of the potential, the SCHA matrix of Eq. (17), plus two terms depending on the third and fourth order SCHA tensors. In component-free notation, we can write Eq. (21) in compact form:
[TABLE]
where the contraction on the indices is understood. Moreover, it is convenient to introduce a ‘super-index’ . In this way, for example, , and are square symmetric ‘super-matrices’ of order , and the contraction of indices between them can be seen as a matrix product.
As explained in the previous sections, the curvature of the free energy in a high-symmetry phase as a function of temperature is essential in order to identify and characterize a second order phase transition. Diagonalizing the real symmetric matrix we obtain eigenvalues and eigenvectors as a function of temperature. In the presence of a second order phase transition, there is at least one eigenvalue that becomes negative at the transition temperature, and the corresponding eigenvector identifies the instability distortion pattern which reduces the free energy. By definition, the SCHA matrix is positive-definite (see comment after Eq. (11)). On the contrary, as shown in Eq. (150), is negative-definite thus is negative-semidefinite. It is this term, which for reasons that will be clear later we call ‘bubble’ (see Sec. VI), that allows the second derivative of the free energy to have negative eigenvalues. The formula obtained for also clarifies in this way the long-standing debate about the possibility of having second order phase transitions within the SCHA Pytte (1972); Gillis and Koehler (1972): the SCHA can describe a second order phase transition only if .
Using the interpretation of the 4th-rank tensors as super-matrices of order , is readily obtained by inverting Eq. (24) in matrix form:
[TABLE]
Substituting Eq. (26) into Eq. (25) we obtain the compact expression for the free energy Hessian:
[TABLE]
This is the equation that has been implemented . It is also interesting to write Eq. (27) in a more symmetric fashion:
[TABLE]
where is the adimensional real symmetric matrix that rules the convergence of the geometric series. For example, from this formula we clearly see that in the limiting case where the absolute values of ’s eigenvalues are much smaller than one and can be discarded with respect to the identity, the curvature is given by the SCHA matrix plus the bubble only.
V Phonons in the SCHA
From the results obtained, it is tempting to use the curvature of the free energy with respect to the centroids to define a phonon-like dispersion. To this purpose, for each temperature, we consider the free energy curvature in the corresponding equilibrium configuration , divided by the square root of the masses:
[TABLE]
This matrix can be considered as the temperature-dependent, free energy-based, generalization of the temperature-independent harmonic dynamical matrix
[TABLE]
Here is the temperature-independent configuration for which the potential has a minimum. We associate the ‘free energy dynamical matrix’ to ‘free energy phonons’, quasi-particles whose energies and polarization vectors are obtained by diagonalization as
[TABLE]
Since is positive-definite if and only if is positive-definite, an instability in the system corresponds to at least a frequency becoming imaginary. It is this fact that justifies the interpretation of as a temperature-dependent generalized dynamical matrix describing temperature-dependent anharmonic phonons. It is worthwhile to emphasize that the theory developed so far is ‘static’, in the sense that it is not based on time-dependent properties, but on the variation of the free energy with respect to a static variation of the centroids position. Moreover, it is important to observe that we cannot use
[TABLE]
to study system instabilities and defines phonon-like particles, even if in some cases it has given temperature-dependent anharmonic phonons in good agreement with experiments Leroux et al. (2015); Errea et al. (2013). Indeed, it is not given by the second derivative of the free energy. Moreover, by definition, is positive-definite, thus it is impossible to observe any softening in its eigenvalues.
The free energy dynamical matrix is a particularly important tool when we consider crystals. Indeed, in that case we can use the same techniques that are standard for the harmonic theory Baroni et al. (2001). Exploiting the translational lattice symmetry, we define the SCHA dynamical matrices in the unit cell as a function of the quasi-momentum . We can explicitly calculate on a coarse grid of the Brillouin zone (BZ) and later Fourier interpolate the result to obtain the matrix on an arbitrary finer grid or a path. Thus, diagonalizing , we obtain the spectrum and the polarization vectors on a path of the BZ. An imaginary phonon in a point indicates that the system is unstable for a distortion with modulation that reduces the lattice periodicity. This is, for example, what happens in charge-density wave instabilities. Therefore, with moderate workload, it is possible to have a complete picture of the crystal instabilities. In particular, with calculations on supercells of reasonable size it is possible, in principle, to study lattice instabilities which are periodic on very large supercells or even incommensurate.
VI Diagrammatic representation
In this section we give a perspicuous diagrammatic description of Eq. (27), in order to reformulate it in a language familiar to the field theorists. The diagrammatic description can also be useful as a basis for further developments of the theory, as we will see later in Sec. IX.
Fixed the temperature, with the corresponding we define the quadratic ‘SCHA Hamiltonian’
[TABLE]
and we consider the corresponding SCHA thermodynamic Green function for the displacements normalized by masses . Since is quadratic
[TABLE]
where indicates the inverse matrix of (similar notation for the inverse will be used later also in other formulas). We also consider , the SCHA ‘static’ loop, i.e. the loop with and total frequency equal to zero:
[TABLE]
where is the -th Matsubara frequency. With standard techniques for Matsubara frequency summation we obtain Mahan (2000); Maradudin and Fein (1962)
[TABLE]
with defined in Eq. (23). From Eq. (22), Eq. (35) and Eq. (36) we obtain a relation between the tensor and the static loop :
[TABLE]
Therefore, using Eq. (29) and Eq. (32), formula (27) divided by the square root of the masses gives
[TABLE]
where, as usual, we have used bold symbols in component-free notation and we have defined
[TABLE]
Here we have generalized the definition (32) to the -th order as
[TABLE]
Notice that we did not use the superscript for the second order tensor defined by Eq. (32). In terms of the SCHA Green function defined in Eq. (34), Eq. (38) is readily written as
[TABLE]
which is equivalent to the Dyson-like equation
[TABLE]
where the matrix product is understood. If the opportune diagram symmetry factors are taken into account, Eq. (42) with Eq. (39) have the Feynman diagrams representation shown in Fig. 2a and Fig. 2b. This is the diagrammatic representation of the curvature formula (27) (divided by the square root of the masses). Analogous diagrammatic series has been obtained by Götze and Michel in Ref. 42.
The first term of the series giving is the SCHA ‘bubble’ . It is given by the formula
[TABLE]
and corresponds to the diagram in Fig. 2c. The SCHA ‘bubble’ is the term of Eq. (25), divided by the square root of the masses. This explains the name ‘bubble’ given to that term.
Before concluding this section, it is worthwhile to remark that, in spite of the symbol used, at this level the defined in Eq. (39) is just an auxiliary quantity, without a specific physical meaning. However, the choice of the symbol is not casual because later we will interpreted it as a self-energy. This will give a deeper meaning to the results obtained.
VII Stochastic implementation
The stochastic implementation of the SCHA (SSCHA) has demonstrated to be an efficient method to analyze thermal properties of solids in situations where the harmonic approximation breaks down Errea et al. (2013, 2015, 2014, 2016); Borinaga et al. (2016); Leroux et al. (2015); Borinaga et al. (2016). The SSCHA is described in Ref. 4 and consists in minimizing with a conjugate-gradient (CG) method the functional \mathcal{F}[\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}] with respect to and . The functional and its gradient are expressed as averages taken with \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}} of observables \mathds{O}(\boldsymbol{R})=\mathds{O}\bigl{(}V(\boldsymbol{R}),\text{{f}}(\boldsymbol{R})\bigr{)} that are functions only of the potential and the forces . The method is ‘stochastic’ because these averages are evaluated with the importance sampling technique. Since the observables depend only on the position, Eqs. (9)–(10) apply. The space of configurations is statistically sampled with a (large) population of finite size , whose members are distributed according to the probability density \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}(\boldsymbol{R}). For each element , being the displacement from the centroids , the forces and the potential energy are calculated by any energy-force engine, i.e., making use of first-principles methods or empirical potentials. In that way the average integrals can be straightforwardly computed. However, at each step of the CG minimization algorithm the distribution probability \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}(\boldsymbol{R}) changes. Thus, in principle, at each minimization step a new population should be generated and for its members the energies and the forces should be calculated. In order to reduce the number of calls to the energy-force engine, in actual calculations a reweighting procedure is adopted Errea et al. (2014). Energy and forces are computed only once for the population elements that are distributed according to an initially fixed probability density \scalebox{1.1}{\tilde{\rho}}{}_{\scriptscriptstyle{\text{in}}}(\boldsymbol{R}). The approximate averages for a generic distribution probability \scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}(\boldsymbol{R}) are then computed as
[TABLE]
Obviously, the equality holds for .
We want to use the stochastic approach also to compute the free energy curvature through Eq. (27). Considering a configuration , after the SSCHA minimization of the functional \mathcal{F}[\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}] with respect to , the SCHA matrix for that configuration is available. Therefore we only need to express and in a form that is suited for the stochastic calculation (here and in what follows the dependence of the matrices on is understood). As demonstrated in Appendix C (see Eqs. (156), (171a) and Eqs. (162), (171b) with Eq. (172)), it can be shown making use of integration by parts that
[TABLE]
Here is the matrix obtained from through the definition (11), and
[TABLE]
The equations (45) express the third and fourth order SCHA tensors in terms of averages of forces and displacements only (in the definition (46) the term subtracted from the forces is computed analytically with negligible cost, since \Bigl{\langle}\partial V/\partial\boldsymbol{R}\Bigr{\rangle}_{{\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\Phi}}}} and are known). Therefore, they can be calculated through Eq. (44).
It is interesting to observe that, in the limit of an infinitely large population sampling, adding to a term odd in the displacements does not change the value of obtained from Eq. (45a). Therefore, the used in Eq. (45a) is actually defined only up to an additive factor that is odd in the displacements. Analogously, if we use an infinite sampling, the used in Eqs. (45b) is defined only up to an additive factor that is even in the displacements. However, depending on the actual used, we obtain different results when we use a finite sampling to compute the averages. The specific choice of Eq. (46), identical for both equations (45), guarantees that if the potential is quadratic, then the SSCHA tensors (i.e. the SCHA tensors calculated stochastically) and are correctly zero with any finite sampling used to compute the averages. Therefore, the definition (46) reduces the stochastic error and accelerates the convergence. Notice that if we compute the curvature of the free energy in a stationary point, since it is then from Eq. (18) the term \Bigl{\langle}\text{f}_{i}\Bigr{\rangle}_{{\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\Phi}}}} in Eq. (46) is zero. In particular, this is true when we evaluate the curvature in the equilibrium configuration , which is the relevant case when we study structural second order phase transitions.
In the limit of a fully converged stochastic calculation, the SSCHA tensors and satisfy both acoustic sum rule (ASR) and invariance with respect to permutations of indices and symmetry transformations. Actually, in Appendix C it is shown that the SSCHA -th order tensor satisfies the ASR with any finite population sampling, as long as the total force acting on the system is zero for any population element (as it must be), and the ASR is satisfied by . Therefore, it is not necessary to impose any extra-condition to make and satisfy the ASR.
For the invariance properties the situation is different. We can distinguish two kind of operators acting on a tensor : the operators , which permute the tensor indices according to the permutations , and the operators , whose action corresponds to the symmetry transformations (excluding lattice translations, if it is a crystal). If we are considering a crystal, the SSCHA calculation is performed on a supercell made of unit cells, with periodic boundary conditions. In that case we consider also the operators whose action corresponds to the translations by lattice vectors non commensurate with the supercell . The SSCHA tensors are invariant with respect to these operations only in the limit (for simplicity we consider the crystal case):
[TABLE]
For calculations performed with finite-size populations these conditions are not satisfied. We enforce them by applying the projectors , and to the result:
[TABLE]
For calculations with finite sampling the action of the projectors (48) has two benefits: we obtain SSCHA tensors with the correct properties and we reduce the statistical noise and improve, with negligible cost, the rapidity of the statistical convergence with respect to . Indeed, the necessity of imposing the property (47a) is due to the fact that Eqs. (45) are not symmetric with respect to permutation of indices. That is caused by the arbitrariness in the choice of the variables integrated by parts in the derivation of the formulas, shown in appendix C. As a consequence, an approximate evaluation of the averages causes spurious asymmetries, which are eliminated by applying the projector to the result. The necessity of imposing the properties (47b) and (47c) is instead due to the fact that, in general, the population generated to compute the averages is composed of elements whose distribution in configuration space does not respect the symmetries of the system. This leads to spurious fluctuations which spoil the symmetry properties of the result and which are eliminated by applying the projectors and . Applying these projectors to the result corresponds to computing the averages through Eq. (44) using a larger population of elements obtained by applying the symmetry operations on the members of the original population.
In conclusion, the formulas implemented in the SSCHA are:
[TABLE]
VIII Perturbative limit
In this section we analyze the lowest perturbative order of the SCHA and of the free energy dynamical matrix . First we set some definitions. Expanding the potential around its minimum , the Hamiltonian is written as
[TABLE]
where is the displacement with respect to the potential minimum,
[TABLE]
is the quadratic harmonic Hamiltonian, and
[TABLE]
is the -th order force constant tensor. Notice that for the second order force constant matrix we do not use the superscript . In order to avoid confusion, it is worthwhile to stress that the -th force constant is the -th derivative of the potential, evaluated at the potential minimum , whereas the -th SCHA tensor , defined in Eq. (20), is the -th derivative of the potential averaged with the distribution \displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\Phi}(\boldsymbol{\mathcal{R}})}.
The part of the Hamiltonian in Eq.(50) not included in defines the anharmonic part of the potential, which we treat as a (small) perturbation of . With and we indicate the Green function of and for the variable , respectively. The latter is given as
[TABLE]
where is the harmonic dynamical matrix, already defined in Eq. (30). The relation between the full and harmonic Green functions is given by the Dyson equation
[TABLE]
which is equivalent to
[TABLE]
where, in order to use a consistent notation, we have indicated with the harmonic self-energy, i.e. the self-energy obtained by taking as non-interacting unperturbed Hamiltonian. At the lowest perturbative order Maradudin and Fein (1962)
[TABLE]
where , and are the loop, tadpole and bubble harmonic self-energies, respectively, which have the following expressions:
[TABLE]
[TABLE]
[TABLE]
Here we have generalized the definition (30) of the harmonic dynamical matrix to the -th order:
[TABLE]
Notice that loop and tadpole self-energies do not depend on the value of the frequency . In fact they are real symmetric. On the contrary, the bubble is a complex symmetric matrix depending on . In Fig. 3 the diagrammatic representation of the harmonic perturbative result at the lowest order is shown.
From the SCHA equations, retaining only the lowest order corrections to the harmonic values and , using the SCHA matrix defined in Eq. (32) we obtain (see Eq. (182))
[TABLE]
Equivalently, using the SCHA propagator defined in Eq. (34) we can write
[TABLE]
that is
[TABLE]
At the lowest perturbative order we also have (see Eq. (185))
[TABLE]
where and are the quantities defined in Eq. (39) and Eq. (43), respectively. From Eqs. (61)–(63) we see that, at the lowest perturbative order, the SCHA and harmonic propagators are related through the harmonic loop and tadpole self-energies only Paulatto et al. (2015). However, from Eq. (38) and Eq. (64) we see that in order to obtain the SCHA dynamical matrix, defined in Eq. (29), we need the harmonic static bubble too:
[TABLE]
Notice that, in particular, this implies that the term in the curvature formula, Eq. (25), can be discarded at the lowest perturbative order. In terms of the harmonic propagator defined in Eq. (53), the formula (65) can be written as (Cfr. Eq. (42))
[TABLE]
Equations (63) and (66) are the main SCHA results at the lowest harmonic perturbative order. They are represented in diagrammatic form in Fig. 4a and in Fig. 4b, respectively.
It is interesting to observe that, at the lowest perturbative order, the free energy curvature takes into account only the static harmonic bubble, whereas in the full propagator the bubble actually depends on the frequency , as we can see from Eq. (56). This is consistent with the fact that we have developed only a ‘static’ theory (obviously, this fact does not have consequences for the tadpole and loop term, because they do not depend on the frequency). In the next section we will investigate possible dynamic extensions of the results found thus far.
IX Ansatz for a dynamic theory
In this section we propose a possible ‘dynamical’ extension of the ‘static’ results obtained above. This could be used to interpret the outcomes of inelastic scattering processes between phonons and external incident particles (typically neutrons) in the framework of the SCHA approximation. The extension that we are going to present is reasonable because it returns the expected results in two limits. In the static limit it gives results coherent with the ones already obtained for the free energy curvature and at the lowest perturbative order it gives the correct results already known in literature. Nevertheless, it is worthwhile to stress that, at variance with the ‘static’ results, the dynamical extension that we are going to propose is only an ansatz, reasonable but not based on a rigorous demonstration. For that reason it can be considered as a basis for a future rigorous extension of the static theory.
Fixed the temperature, and the relative , we consider the full Green function for and the Green function for in the variable . We consider a Dyson-type relation between them:
[TABLE]
which is equivalent to
[TABLE]
where is the SCHA self-energy. The aim of this section is to propose an expression for . The first assumption is that its static value, i.e. its value for , is given by Eq. (39). At that level, the symbol used did not have a physical meaning. Now we are explicitly interpreting it as the static SCHA self-energy. Comparing Eq. (67) to Eq. (41), this is equivalent to saying that
[TABLE]
This is the same kind of relation that exists between the harmonic static Green function and the harmonic dynamical matrix. Therefore, Eq. (69) gives a deeper meaning to the consideration in Sec.V that is the anharmonic generalization of the harmonic dynamical matrix. A real pole of the Green function corresponds to the energy of a phonon with zero linewidth, i.e. with infinite lifetime. Equation (69) means that we observe a phonon with zero energy, i.e. we see a phonon softening and therefore an instability, when has a null eigenvalue. This is exactly the result found in Sec. IV and Sec. V. Thus the interpretation of Eq. (39) as the static SCHA self-energy is consistent with the rigorous (static) results obtained for the free energy curvature.
The subsequent step is to give an expression for the SCHA self-energy at different from zero. As a second part of our hypothesis, we assume for the same structure of , given by Eq. (39) and illustrated by the diagrams in Fig. 2b, but readily generalized to any . Therefore it is
[TABLE]
with
[TABLE]
Using standard techniques for Matsubara frequencies summations Mahan (2000), we obtain an explicit expression for this term:
[TABLE]
where and are eigenvalues and corresponding eigenvectors of , respectively, and for
[TABLE]
The assumption expressed by Eqs. (70), (71) is reasonable because at the lowest perturbative limit it gives the correct result. Indeed, by using the same arguments of Sec. VIII, at the lowest perturbative order we readily generalize Eq. (64) to
[TABLE]
Thus, from Eq. (62) and Eq. (67) we obtain
[TABLE]
which is the correct perturbative result shown in Eqs. (54) and (56). In conclusion, according to our ansatz, the full Green function is (approximately) given by Eq. (68) with Eqs. (70), (71). In that way we obtain a minimal extension of the static theory which reproduces the correct instabilities and gives the correct results at the lowest perturbative level. By using this formula we can study anharmonic effects in a non perturbative way also for the dynamic case. In Fig 5 we give the diagrammatic expression for our ansatz, the self-energy being the one in Fig. 2b. An analogous diagrammatic series has been proposed in Ref. 45.
It is interesting to observe that, inspired by the perturbative result in Eq. (56), one could be tempted to naively obtain a dynamic SCHA theory simply by adding a dynamic bubble term on top of the standard SCHA results (which, as shown in Eq. (62), contain only tadpole and loop at the lowest perturbative level). This approach of adding a dynamic bubble has been taken, for example, in PbTe Li et al. (2014) and PdH Paulatto et al. (2013), where the strong anharmonicity induces satellite peaks in the spectral function. Now we can see that this essentially consists in adopting our ansatz, but discarding all the terms in described by the diagrams of Fig. 2b, except the non-perturbative SCHA dynamic bubble given in Fig. 2c:
[TABLE]
This, in general, is not justified. As long as we consider a non-perturbative situation, there is in principle no hierarchy that allows to discard the other terms. Therefore, the term given by Eq. (76) has to be considered an incomplete expression for and a better choice is to take into account the full expression of Eq. (70). Of course, there can be situations in which even if the regime is not perturbative, because the third order is not smaller than the harmonic term, nevertheless the superior orders are smaller. In that case it would be justified to use Eq. (76) to evaluate . However, this is a further assumption that, in order to be adopted, has to be justified case by case.
X Numerical test
In order to give a numerical demonstration of our findings, we apply the theory to a toy model based on the SnTe crystal (an analogous model could be used for GeTe). SnTe crystallizes at room temperature and ambient pressure in the NaCl-structure (Fm-3m), called -SnTe phase, where two fcc lattices of Sn and Te interpenetrate. At low temperature, around K, it undergoes a phase transition and stabilizes in a rhombohedral structure (R3m), called -SnTe. The phase transition can be described in terms of a two-step symmetry reduction: a fixed unit cell polar displacement, between the two fcc, along the [111] cubic direction, which eliminates the inversion center, and a strain of the unit cell along the cube diagonal Bauer Pereira et al. (2013). We concentrate on the first distortion. We define the interatomic potential of the toy-model as a function of the displacements from the equilibrium position of the rock-salt structure and we keep, beyond the quadratic part, only the anharmonic third and fourth order terms:
[TABLE]
The harmonic matrix has been obtained from first principle calculation for SnTe on a 2x2x2 grid of the Brillouin zone (BZ) (details in App. E). With the experimental lattice parameter we do not observe any instability in the total energy (i.e. the harmonic matrix is positive-definite). However, a lattice instability appears and increases at BZ as we increase the lattice parameter. Therefore in order to achieve, for explicative purposes, an increased instability at the harmonic level, we calculated ab initio the harmonic matrix with a higher lattice parameter: . Moreover, in order to keep the toy model as simple as possible and focus on the main purpose of the numerical test, we ignored the LO-TO splitting at , which is present in real undoped SnTe samples. In Fig. 6 we show the obtained (harmonic) phonon dispersion along a high-symmetry path of the fcc BZ. There are imaginary phonons in several points, the optical phonon in corresponding to the highest instability.
For the third and fourth order contributions, we follow the model described in Refs. 56; 57. We define short-range anharmonic terms by using reciprocal displacements of nearest-neighbor atoms (in the rock-salt structure each atom has 6 nearest-neighbors). In particular, as explained in App. E, in our model is proportional to a single parameter , and is a linear function of two parameters , . We take , and .
X.0.1 Free energy curvature
We consider the free energy profile obtained by displacing the atoms in the unit cell along the [111] cubic direction. In order to describe this distortion we write the atomic position as a function of a scalar, adimensional parameter :
[TABLE]
where is the configuration corresponding to the minimum of the potential energy along the distortion path. Therefore, is linear, and corresponding to the high symmetry phase (Fm-3m) and to the low-symmetry energy minimum (R3m), respectively. In Fig. 7 we show , the variation of the potential (per unit cell) along this distortion path. This curve depends on , , . The harmonic term is responsible for the initial decrease whereas the fourth order term gives the subsequent increase. On the contrary, due to the symmetry of the rock-salt structure, the value of is not relevant for the energy pattern (as a matter of fact, the value of does not affect the energy value of any unit cell configurations).
The free energy along the path has been calculated with the SSCHA on a 2x2x2 supercell. Fixed and the temperature, the GB functional \mathcal{F}[\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}] has been minimized with respect to as described in Ref. 3; 4. In Fig. 7 we show a complete variation path for the free energy at three temperatures. For reasons that will be clear in a while, we studied also the case without third order (). However, it is interesting to remark that at the SCHA result is independent from .
A first remarkable, somewhat counterintuitive, conclusion can be deduced from the results of Fig. 7. While the potential energy path is independent from , at given temperature the two free energy paths obtained with and are considerably different. This has important consequences. The presence or not of a second order phase transition and, when there is such a transition, the transition temperature and the low-symmetry equilibrium configuration for are properties which cannot be inferred from the potential energy profile.
From the values of the free energy computed near , the curvature in the origin has been evaluated by finite difference. The results at four temperatures are shown with dots in Fig. 8. We compare these values with the curvature in calculated by contracting with the formula for of Eq. (25):
[TABLE]
This formula is evaluated at . In order to be consistent with the finite differences result, all the ingredients have been calculated by using the SSCHA on a 2x2x2 supercell. Once the SSCHA minimization at has been completed and the converged value for has been obtained, and have been computed using Eq. (49). For each temperature, we used the converged value of to generate the population used to compute the averages. Therefore, in this case it is \scalebox{1.1}{\tilde{\rho}}{}_{\scriptscriptstyle{\text{in}}}(\boldsymbol{R})=\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\Phi}}(\boldsymbol{R}). Notice that, as explained in Sec. VII, since the calculation has been performed in a stationary point of the free energy, the term \Bigl{\langle}\text{f}_{i}\Bigr{\rangle}_{\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\Phi}}} on the right-hand side of Eq. (46) is zero.
For explicative purposes we have used Eq. (25) to express the curvature. Thus we have three terms in Eq. (80) and in Fig. 8 we plot three lines to show their different contributions. The term obtained from does not depend on the value of , whereas the other two terms depend quadratically on . As a consequence, gives the curvature in the high-symmetry phase when the third order is absent. This is confirmed, within the statistical error (), by comparing the red curve and the red dots in Fig. 8. For the other two terms and are necessary in order to obtain the curvature. This is confirmed by comparing the blue curve and the blue dots in Fig. 8. As explained in Sec. VIII, only at the lowest perturbative order it is possible to neglect the term . Indeed, in Fig. 8 we show with a yellow line the curvature computed with only the SCHA matrix and the bubble. In this case the difference with respect to the correct value increases with temperature and, even if small, it is already beyond the statistical error at K.
In this section we have numerically proved the correctness of Eq. (25). We conclude with a consideration. As already stressed in Sec. IV, the first term of Eq. (80) is always positive. Therefore, with it is possible to observe only a first-order phase transition within the SCHA approximation. With , the plot in Fig. 8 shows that the free energy curvature in changes sign for K. However, in Fig. 7 we see that at the free energy has already developed a lower minimum in . As a consequence, the toy model studied undergoes a first order phase transition even with different from zero.
X.0.2 Phonons
In this section we apply the concept of free energy dynamical matrix defined in Sec. V. To be precise, fixed the temperature, we compute the second derivative of the free energy in , divided by the square root of the masses. Notice that, properly speaking, this is only for , when is equal to , because at temperatures below the transition temperature departs from . Nevertheless, for explicative purposes and having this caveat in mind, we will use the same symbol even at .
The matrix is given by the matrix plus the static self-energy which, in turn, is made of the bubble term plus other factors, negligible at the lowest perturbative level (see Eqs. (38), (39), (43), (64)). Since we are considering a crystal, we exploit the lattice translational symmetry and we write the dynamical matrices in the unit cell as a function of the quasimomentum. In Fig. 9 we plot the spectrum of these matrices along a high-symmetry path of the BZ. We consider two temperatures. The matrix coincides with the free energy dynamical matrix when the third order is absent. Since is positive-definite, the spectrum is always positive. However, with the dynamical matrix is qualitatively different from . Below the transition temperature the phonon spectrum becomes imaginary (negative eigenvalue) in , and only in that point. The other instabilities that were present in the harmonic phonon spectrum, Fig. 6, have been washed out by the zero-point energy and anharmonicity. Notice that in this case the comparison between the harmonic and the free energy dynamical matrix is particularly meaningful because, for symmetry reasons, both are computed in the same point . In Fig. 9 we also show the spectrum obtained by adding only the bubble to . From the results shown in Fig. 8 we expected in a very small difference between the full formula and the one considering only the bubble. However, here we have a more complete picture. As we can see, in other points of the BZ the spectrum is more affected by the presence of terms beyond the bubble. For example, at K for the 5th mode in , the terms beyond the bubble change the spectrum around .
X.0.3 Convergence
Since for our test we used a toy model, i.e. an analytic potential, we could evaluate the averages using populations of very big size at small computational cost. However, in view of first principle applications for realistic materials, we carefully performed convergence tests of the curvature formula with respect to the population size .
First, we tested the convergence of at various temperatures. As said, for each temperature we calculated the curvature using the converged value of to generate the population used to compute the averages in Eq. (49). As shown in the upper left-hand panel of Fig. 10, the convergence can be considered reached with . However, it is worthwhile to say that, in general, fitting the values of the curvature versus temperature with a polynomial allows to wash out part of the stochastic noise and obtain good estimations for with smaller populations. In this case, for example, fitting with a th degree polynomial the results obtained with gives a value for which is only K smaller than the converged one.
As we have seen in Fig. 8, for the terms beyond the bubble, which depend on , have a limited relevance. For that reason, we performed an analogous convergence test for the frequency of the th mode in of . Indeed, as shown in Fig. 9, for that specific mode the terms beyond the bubble play a non negligible role in the determination of the spectrum. Therefore, this quantity is particularly significant to analyze the convergence of the different terms comprising the curvature formula. Here, as in the previous paragraph, with we are indicating the curvature of the free energy in divided by the square root of masses, even at temperatures below . As shown in the upper right-hand panel of Fig. 10, also in this case the convergence can be considered reached with . However, the absolute stochastic error is already smaller than with .
It is interesting to see how the two terms and affect the convergence, separately. To that end, we plot in the the other panels of Fig. 10 the curvature and the frequency of the chosen mode, versus temperature, obtained once with computed with different population sizes but with fixed to the converged value (obtained with a population of elements), and in the other case the inverse. The conclusion is that the total convergence is affected in a similar way from the two tensors and . This could be surprising since the 4th order tensor is obtained by averaging a quantity that depends three times on the displacements, whereas for the 3rd order tensors it is averaged a less complicated quantity which depends only two times on the displacements. However, it has to be considered that in the curvature formula is fully contracted (at variance with ). Indeed, the random fluctuations on the single components tend to cancel each other and, thus, the convergence of a contracted tensor is expected to be faster than the convergence of a single tensor component.
XI Conclusions
In this work we present an approach to study structural second order phase transitions in molecules and solids within the self-consistent harmonic approximation. The developed method allows to estimate transition temperature and instability modes. It is based on the analytic formula giving the second derivative of the SCHA free energy with respect to the average atomic positions. The Hessian of the SCHA free energy is also expressed in terms of thermal averages of forces and displacements. Therefore, the method is suitable for a stochastic implementation in conjunction with any energy-force engine. Considering a configuration, it allows to calculate directly the free energy curvature once the SSCHA calculation has been performed in that point. As a consequence, it permits to avoid the very computational demanding finite difference approach of computing the curvature through several SSCHA calculations for different configurations Errea et al. (2016). Moreover, the imposition of symmetries on the result reduces the statistical noise and speeds up statistical convergence with respect to the population size used to compute the averages with the importance sampling technique.
The efficiency of this method makes it ideal to be used in conjunction with first principle energy-force engines to study realistic materials, such as ferroelectrics or CDW materials. With the curvature formula it is possible to find the instabilities of a general condensed matter system. In particular, the method is especially convenient for crystals, since by exploiting the lattice translational symmetry and the Fourier interpolation technique it is possible to find distortions lowering the free energy with any modulation in space (e.g. periodic on large supercells or even incommensurate) with SSCHA calculations performed on supercells of moderate size. In order to demonstrate our findings, numerical tests have been performed on a toy model. The results confirm both correctness of the theory and numerical efficiency of the implemented method.
In addition to its practical utility, the developed theory sheds light on several fundamental aspects of the SCHA. In particular, the role of the auxiliary effective quadratic Hamiltonian is clarified. It is shown that the SCHA matrix is only a term of the free energy Hessian and, in general, it does not define an anharmonic dynamical matrix. On the contrary, an anharmonic temperature dependent, free energy based, dynamical matrix is obtained through the free energy curvature. It generalizes the temperature independent harmonic dynamical matrix and defines temperature dependent anharmonic phonons.
The theory developed for the SCHA free energy curvature is static, as it does not take into account any dynamical effects. Inspired by a perspicuous diagrammatic interpretation of the results, we propose a tentative minimal dynamic extension of the static theory in order to associate spectral functions with anharmonic phonons and interpret the results of scattering processes in a full non-perturbative way. Similarly, the dynamic theory allows to calculate phonon lifetimes in the non perturbative limit. At variance with the curvature formula, the suggested dynamic extension is not based on a rigorous demonstration. Nevertheless, it is expected to give good results, because it is correct in both static limit and lowest perturbative limit, and thus it opens the way to further theoretical developments and interesting applications.
ACKNOWLEDGMENTS
The authors acknowledge support from the Graphene Flagship. I.E. acknowledges financial support from the Spanish Ministry of Economy, Industry, and Competitiveness (Grant No. FIS2016-76617-P). M.C. acknowledges support from Agence Nationale de la Recherche under contract ANR-13-IS10-0003-01, from the Graphene Flagship, PRACE for awarding us access to resource on Marenostrum at BSC and the computer facilities provided by CINES, IDRIS, and CEA TGCC (Grant EDARI No. 2017091202).
Appendix A Proofs of the SCHA method
In this appendix we give an explicit demonstration of the SCHA self-consistent equation, Eq. (17), and of the expression for the first and second derivative of the SCHA free energy, Eq. (18) and Eq. (21), respectively. We find convenient to use a notation slightly different from the one used in the main text. Considering a trial harmonic matrix we define three matrices from it: the matrix ; the matrix , where and are eigenvalues and eigenvectors of , respectively, and , where is the bosonic average occupation number; and the matrix . Introducing the diagonal mass matrix , we can summarize these definitions in the compact form
[TABLE]
In this notation the matrix of Eq. (11) coincides with the inverse of , which we indicate with the symbol . Thus the thermal average of an observables that is function only of the position is given by Georgescu et al. (2013)
[TABLE]
For the subsequent derivation it is convenient to perform the change of variable
[TABLE]
where we have introduced the compact notation . With that change of variable, the average is written as a Gaussian integral:
[TABLE]
where
[TABLE]
We now demonstrate the following two relations:
[TABLE]
From the explicit expression of the average, the first identity is trivially obtained. In order to demonstrate the second one, we use integration by parts:
[TABLE]
the boundary terms at infinity being zero due to the exponential.
Denoting with the free energy of \displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}, we prove next the following relations:
[TABLE]
where is the trial potential. Indeed, since the trial Hamiltonian is quadratic, we have the standard result
[TABLE]
Therefore, since ,
[TABLE]
The matrices and have same eigenvectors but eigenvalues and , respectively. Thus,
[TABLE]
or, equivalently,
[TABLE]
Moreover,
[TABLE]
Indeed,
[TABLE]
From Eq. (94) and Eq. (95) the relation (88b) is readily obtained. The equation (88a) comes from the observation that neither nor \left\langle{\widetilde{V}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}}\right\rangle_{\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}} depend on .
The SCHA functional, i.e. the GB functional restricted to quadratic trial Hamiltonians, can be written as Errea et al. (2013, 2014)
[TABLE]
From this relation, and Eq. (86) and Eq. (88) we obtain:
[TABLE]
Fixed , with we indicate the matrix that minimizes \mathcal{F}[\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}] with respect to . This implies that \partial\mathcal{F}[\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\varPhi}}]/\partial\boldsymbol{\varPhi} is equal to zero in :
[TABLE]
Therefore, from Eq. (98b) we have the self-consistent relation
[TABLE]
Defining F(\boldsymbol{\mathcal{R}})=\mathcal{F}[\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\Phi}(\boldsymbol{\mathcal{R}})}], from Eq. (98a) and Eq. (99) we have
[TABLE]
Deriving one more time and using Eq. (86) and Eq. (100),
[TABLE]
where is defined as a generalization of Eq. (100) to higher orders (see Eq. (20)) and in the application of the chain rule we have derived only with respect to the independent components of the symmetric matrix . Moreover, using Eq. (100),
[TABLE]
In the next section we will give the explicit expression of the 4th order tensor satisfying the relation
[TABLE]
Using it we rewrite Eq. (103) and Eq. (104) in the following way:
[TABLE]
The second equation can be solved by iteration, in the hypothesis that the resulting series converges:
[TABLE]
Substituting this solution into Eq. (106) we obtain
[TABLE]
where is the 4th order tensor given by the series
[TABLE]
Thus, solves the Dyson-like equation
[TABLE]
already introduced in Eq. (24).
Appendix B The matrix
In this section we derive the explicit expression of the tensor used in Eq. (105). We break the derivation into several intermediate steps.
B.1 Derivatives of eigenvalues and eigenvectors with respect to matrix elements
Let us consider a real symmetric matrix with distinct eigenvalues and eigenvectors . From non degenerate perturbation theory, if depends on a parameter we have
[TABLE]
where . In particular, for ,
[TABLE]
Thus,
[TABLE]
and
[TABLE]
B.2 Calculation of
The matrix can be written as
[TABLE]
Given a regular function , is a real symmetric matrix having the same eigenvectors and eigenvalues :
[TABLE]
Then, we can write
[TABLE]
with
[TABLE]
and
[TABLE]
Therefore, from Eq. (115)
[TABLE]
and from Eq. (116)
[TABLE]
B.3 Contraction of with a symmetric matrix
We are interested in calculating the quantity
[TABLE]
with a symmetric matrix. First, we define the tensor
[TABLE]
with
[TABLE]
and
[TABLE]
In terms of the tensor
[TABLE]
Therefore,
[TABLE]
Moreover, notice that
[TABLE]
Therefore, given a symmetric tensor we have
[TABLE]
with defined by
[TABLE]
where
[TABLE]
Notice the following two symmetries:
[TABLE]
B.4 Degeneracy
Up to now, we have exclusively considered the non-degenerate case. However, if we write in the form
[TABLE]
Eq. (132) is well defined even in case of degeneracies, and gives the same tensor regardless of the gauge chosen, i.e. regardless of the basis set chosen in the degenerate spaces. In order to see that, let us consider two eigenvectors basis sets and related as
[TABLE]
where is an orthogonal matrix, i.e. it satisfies
[TABLE]
and it does not mix eigenvectors with different eigenvalues:
[TABLE]
Notice that an immediate consequence is that
[TABLE]
Then,
[TABLE]
Therefore, Eq. (132) and Eq. (136) are well defined even in the degenerate case, because they give the same result, regardless of the specific basis of eigenvectors chosen. This essentially proves that Eq. (131), is valid even if is degenerate. In fact, given a matrix degenerate, we can consider a real symmetric matrix function continuously depending on a parameter ,
[TABLE]
such that but with non degenerate for . Notice that, given , the arbitrariness in the choice of the function reflects in the arbitrariness of the eigenvectors for , i.e. in the choice of a specific gauge. We can apply Eq. (131), Eq. (132), and Eq. (133) with and consider the limit of the result for that goes to zero. That gives Eq. (131), Eq. (132), and Eq. (136) for the specific set . However, as observed, the formula is gauge invariant, i.e. the specific choice of in the degenerate subspaces is immaterial. Therefore, the procedure is well defined, as it gives a unique result.
B.5 The matrix
We apply the general results found to Eq. (105). We can write:
[TABLE]
Since and the matrix is symmetric in , applying Eq. (131) we obtain
[TABLE]
where
[TABLE]
and
[TABLE]
Here and are eigenvalues and eigenvectors of , respectively. Since and , we can write
[TABLE]
with
[TABLE]
Therefore, we obtain the final expression:
[TABLE]
For some derivations it is convenient to express the 4th-rank tensor as a square super-matrix , with and . From Eq.(135) we have that is real symmetric (we are using bold symbol in component free notation). Moreover, it is negative-definite. In fact, is monotonically decreasing. Thus from Eq. (148) , and, considering a vector we have from Eq. (145)
[TABLE]
Appendix C Stochastic calculation of SCHA matrices
Given an observable , it can be proved that
[TABLE]
where . In order to demonstrate this formula we use the change of variable in Eq. (83) and the inverse matrix of :
[TABLE]
The demonstration is obtained with integration by parts:
[TABLE]
We apply Eq. (151) to find expressions for and suited for a stochastic calculation. The goal is to express them as averages of functions of positions and forces only. In what follows the dependence of the matrices on is understood.
C.0.1 Stochastic formula for
Applying Eq. (151) we obtain
[TABLE]
where in the last line we used that \Bigl{\langle}u^{p}u^{q}\Bigr{\rangle}_{\displaystyle\scalebox{1.1}{\tilde{\rho}}_{\boldsymbol{\mathcal{R}},\boldsymbol{\Phi}}}=\varPsi^{pq}. In terms of the forces , we write
[TABLE]
Since the average of a function that is odd in the displacements is zero, we can write
[TABLE]
with
[TABLE]
where are generic odd functions. However, this is true only in the limit of an infinite population sampling. In actual calculations we compute the averages with populations of finite size. In that case, the value of obtained from Eq. (156) depends on the specific chosen, i.e. on the specific expression of . In order to reduce the statistical noie and speed up the convergence, we found convenient to utilize
[TABLE]
Indeed, this corresponds to , where is the difference between the true potential and the quadratic “reference” potential
[TABLE]
If is quadratic, then it coincides with , thus and . Therefore, considering Eq. (156) with Eq. (158) implies that if is quadratic then is identically zero, as it must be, with any finite sampling.
C.0.2 Stochastic formula for
With passages analogous to the ones used in the previous demonstration, we obtain:
[TABLE]
Therefore, we can write
[TABLE]
Again, since the average of a function that is odd in the displacements is zero, we can write
[TABLE]
with
[TABLE]
where are generic even functions.
In order to reduce the statistical noise and accelerate the convergence, in our simulations we calculated with Eq. (162) and we defined as in Eq. (158), since it is compatible with Eq. (163). the definition (158) for . In that way, if is quadratic, we obtain for any finite sampling, as it must be. Once has been calculated, another possibility, in principle better, would be to define
[TABLE]
With this choice it is sufficient that is cubic in order to obtain , as it must be, with any finite sampling. Indeed, Eq. (164) corresponds to , where is the difference between and the cubic “reference” potential
[TABLE]
If is cubic it coincides with this , thus calculated with Eqs. (162) and (164) is equal to zero with any finite sampling. However, for the simulations of this paper we did not use the definition (164).
C.0.3 Acoustic sum rule
The SCHA tensor satisfies the acoustic sum rule (ASR) if the sum over any atomic index vanishes. Considering that we are using a double (cartesian,atom) index , in our notation this means that
[TABLE]
with a global translation of the system by the 3D vector . The averages in Eq. (156) and Eq. (162) are evaluated stochastically with a finite-size population through Eq. (44). We demonstrate that if the matrix satisfies the sum rule and the total force on the center of mass of the system is zero for any population member (as it must be), then the approximate SCHA tensors given by Eqs. (156) and (162) with Eq. (158) satisfy the ASR with any finite population sampling. The two above mentioned conditions are expressed by the relations
[TABLE]
thus from Eq. (158) we also have
[TABLE]
This proves that and . The proof that the SCHA tensors given by Eq. (156) and Eq. (162) satisfy the ASR also on the other indices is obtained by considering that
[TABLE]
Indeed, Eq. (168) means that is null eigenvector for , which is equivalent to saying that is null eigenvector for . However, indicating with the inverse of defined in Eq. (81b), it is \overset{\scriptscriptstyle{-1}}{E}_{ab}=\xi^{\scriptscriptstyle{-2}}\bigl{(}D_{ab}\bigr{)} and . Therefore, is null eigenvector for , i.e. is null eigenvector of , which proves Eq. (170).
Notice that computed with Eq. (162) satisifes the ASR with any finite sampling even if we use the definition (164) for . In that case, in order to demonstrate the condition (169), we just have to consider that satisfies the ASR.
In conclusion, and using the matrix (see Eq. (11)) in order to ease the connection with the main text, the formulas that we used to compute the SSCHA 3rd and 4th order tensors are:
[TABLE]
where
[TABLE]
Appendix D perturbative limit of SCHA
In this section we explicitly calculate the lowest perturbative order of the SCHA results. The formalism used is the one introduced in Sec. VIII. We treat the anharmonic potential as a (small) perturbation of the harmonic Hamiltonian , with of order , being the adimensional perturbative expansion parameter (which can be estimated as the ratio between the mean harmonic displacement and the nearest-neighbor atomic distance Maradudin and Fein (1962)).
We consider the SCHA equilibrium position , Eq. (15), and the corresponding SCHA square matrix , Eq (17). From Eq. (19) we have
[TABLE]
Using the explicit expression Eq. (50) and Eq. (51) we obtain
[TABLE]
and
[TABLE]
where and is the matrix obtained from according to the definitions (81). At the lowest order
[TABLE]
and
[TABLE]
where is the matrix related to in the same way as is related to according to the definitions (81). Inverting Eq. (177) we obtain
[TABLE]
where is the inverse matrix of . By substituting Eq. (179) into Eq. (178), we obtain
[TABLE]
Denoting with and eigenvalues and eigenvectors of , respectively, with standard techniques for Matsubara frequency summation we obtain Mahan (2000)
[TABLE]
where is the harmonic Green function for the variable (see Eq. (53)). Therefore, dividing Eq. (180) by the square root of masses and considering the definition of , Eq. (32), and Eqs.(57)–(59) we obtain
[TABLE]
Moreover from Eq. (20) we readily have
[TABLE]
and
[TABLE]
Thus, from Eq. (39), Eq. (43), and Eq. (59), we have
[TABLE]
Therefore, from Eq. (38) and Eq. (182) we obtain
[TABLE]
Appendix E Toy model definition
In this section we define the toy model used in our numerical tests, which we rewrite here as:
[TABLE]
where , being the equilibrium configuration of the rock-salt structure. The harmonic matrix has been obtained with ab initio calculations for SnTe, on a 2x2x2 mesh of the BZ, performed within Density Functional Perturbation Theory (DFPT) Baroni et al. (2001) as implemented in Quantum ESPRESSO Giannozzi et al. (2009). In order to have increased harmonic instability, calculations have been performed with lattice parameter , which is higher than the experimental value . We used Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996), projector augmented wave (PAW) Blöchl (1994) potentials. We calculated Khon-Sham states, with a cutoff of 28 Ry and 280 Ry for the wave functions and the charge density, respectively. The BZ integration has been performed with a Monkhorst-Pack grid Monkhorst and Pack (1976) of 16x16x16 points. The self-consistent solution of the Kohn-Sham equations was obtained when the total energy changed by less than Ry.
In order to describe the anharmonic terms we follow the model described in Ref. 56; 57 and we define short-range anharmonic terms by using reciprocal displacements of nearest-neighbor atoms (in the rock-salt structure each atom has 6 nearest-neighbor). The third and fourth order terms are given by
[TABLE]
and
[TABLE]
with, for example,
[TABLE]
where and are the nearest-neighbor of the atom , along the cartesian direction and , respectively. Similar notation is used for the directions and . According to this definition, the third-order is proportional to the parameter and the fourth-order is linear function of the parameters and . In order to set reasonable values, and have been fixed by fitting the energy curve obtained ab initio for displacements with unit cell periodicity (the third order does not affect the value of the potential for configurations having unit cell periodicity): and . For the third order we set a value larger than the one reported in Ref. 56 for PbTe. This has been done to magnify the effect of the third order in the 2x2x2 supercell used for the SSCHA calculation.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Born and Huang (1954) M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Oxford University Press, Oxford, 1954).
- 2Maradudin and Fein (1962) A. A. Maradudin and A. E. Fein, Phys. Rev., 128 , 2589 (1962).
- 3Errea et al. (2013) I. Errea, M. Calandra, and F. Mauri, Phys. Rev. Lett., 111 , 177002 (2013).
- 4Errea et al. (2014) I. Errea, M. Calandra, and F. Mauri, Phys. Rev. B, 89 , 064302 (2014).
- 5Errea et al. (2015) I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Phys. Rev. Lett., 114 , 157004 (2015).
- 6Errea et al. (2016) I. Errea, M. Calandra, C. J. Pickard, J. R. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Nature, 532 , 81 (2016) , ISSN 0028-0836, letter. · doi ↗
- 7Borinaga et al. (2016) M. Borinaga, P. Riego, A. Leonardo, M. Calandra, F. Mauri, A. Bergara, and I. Errea, Journal of Physics: Condensed Matter, 28 , 494001 (2016 a) .
- 8Rousseau and Bergara (2010) B. Rousseau and A. Bergara, Phys. Rev. B, 82 , 104504 (2010).