Introduction

Phonons constitute an important topic in condensed matter physics and are required to explain phenomena ranging from thermal expansion to electrical resistivity to BCS superconductivity.1,2 Due to advancements in supercomputing and density functional theory (DFT), first-principles calculations of phonon properties are now routine.3 Within the harmonic approximation many phonon characteristics can already be calculated, including the phonon dispersion and density of states (DOS), specific heat, Helmholtz free energy, and thermal expansion.3 These calculations are even permissible at the DFT level for systems with hundreds of atoms, one example being single-layer black phosphorus or blue phosphorus adsorbed onto Au (111).4 However, some properties cannot be described within the quasi-harmonic approximation: Phonon anharmonicity results in phonon–phonon coupling and is required to describe frequency shift, lifetime, and thermal conductivity (all temperature-dependent). Anharmonic properties can be determined using many-body perturbation theory5,6,7,8,9,10,11,12,13 or density functional perturbation theory.14,15,16,17 For even the simplest materials the former method can require hundreds or thousands of accurate calculations on a large supercell using the finite displacement method. For example, with the phono3py computer code13 the four-atom primitive cell of single-layer black phosphorus18 requires 4741 separate DFT calculations on a 256 atom supercell without using a cutoff distance to describe the three phonon processes. Quartic order treatment may also be necessary to obtain quantitative agreement with experiment.9,12 Furthermore, at sufficiently high-temperature perturbation approaches break down.

On the other hand when molecular dynamics (MD) simulations are performed, phonon anharmonicity is treated exactly to all orders in principle19 and is only limited by how accurately the exchange correlation functional (quantum)20 or parameterized force-field21 reproduce the actual energy surface. Regardless of the level of theory used to calculate the forces, the dynamics are still classical, but the quantum effect of the zero point motion can be included by a simple rescaling of the MD temperature.22 The idea of obtaining vibrational properties from MD goes back to at least 1969: Temperature-dependent vibrational frequencies were extracted from the peak positions in spectra calculated from the Fourier transform of the velocity auto-correlation. It was also suggested early that the decay of the velocity auto-correlation of the normal mode coordinates is related to the phonon lifetime.19 Later the phonon lifetime was extracted from the inverse of the full-width at half-maximum (FWHM) of the spectra peaks fitted with Lorentzians.22,23,24 These methods have become popular and comprise fully dedicated analysis tools such as DynaPhoPy.25 An extension of the auto-correlation method was proposed recently and involves the ratio of correlations of conjugate variables.26 Furthermore, spectra obtained from the velocity auto-correlation without projections are related to the phonon DOS,27 while spectra from the polarizability tensor auto-correlation give the Raman scattering cross section.28 Another approach focuses on spectra derived from the kinetic energy and is termed the spectral energy density (SED) method.29 In the temperature-dependent effective potential (TDEP) method, the harmonic, cubic, and quartic force constants can be obtained directly by minimizing an error function based on all forces from the trajectory, but the phonon lifetimes are still evaluated using perturbation theory.30,31,32,33 On the other hand, self-consistent phonon methods offer non-perturbative results.34,35,36,37,38,39,40,41 In these last two methods the benefit is that deviation of the dynamics lineshapes from what is expected in theory is not an issue, although determining the convergence properties with the simulation time analytically could be challenging.

In this paper, we derive exact analytical expressions for spectra obtained from MD starting with a simple perturbed normal mode expression for the velocities. The effect of anharmonicity-induced frequency shift and lifetime as well as the simulation time are included. The derived fitting functions are applied to a simple model, graphene, hexagonal boron nitride (hBN), and silicon to study the convergence of the vibrational frequency and lifetime with simulation time. It is expected that the proposed method could have far reaching impact on the determination of phononic properties from MD, in areas ranging from strongly correlated systems to biological materials.

Results

Normal mode analysis

Consider a supercell subdivided into primitive cells each indexed by index ρ with corresponding atoms between primitive cells indexed by N. The atoms move about their equilibrium positions x0,ρN as time t evolves:

$${\boldsymbol{x}}_{\rho N}(t) = {\boldsymbol{x}}_{0,\rho N} + {\boldsymbol{d}}^{\rho N}(t).$$
(1)

Here xρN(t) are the time-dependent positions and dρN(t) are the time-dependent displacements from equilibrium. The potential energy V can be expanded in a Taylor series about the equilibrium configuration in the displacements of the atoms:

$$V = E_0 + \frac{1}{2}\mathop {\sum}\limits_{\rho ,\rho \prime } {\mathop {\sum}\limits_{N,N\prime } {\mathop {\sum}\limits_{\alpha ,\alpha \prime } {\left( {\frac{{\partial ^2V}}{{\partial x_{\rho N\alpha }\partial x_{\rho \prime N\prime \alpha \prime }}}} \right)} } _0} d^{\rho N\alpha }d^{\rho \prime N\prime \alpha \prime } + \cdots .$$
(2)

The ground state energy is E0, the Cartesian components are represented by the index α, and the second partial derivative is evaluated at equilibrium. In the harmonic approximation the higher order terms not shown in Eq. (2) are ignored and normal modes arise (indexed by n).42 The general solution for the time-dependent displacements can be written in terms of a superposition of the normal modes:42

$${\boldsymbol{d}}^{\rho N}(t) = \mathop {\sum}\limits_n {\mathop {\sum}\limits_{m = - 1}^1 {{\int}_{\!{BZ}} {\mathrm{d}} } } {\boldsymbol{q}}{\mathrm {e}}^{i{\boldsymbol{q}}\cdot {\boldsymbol{x}}_{0,\rho N}}a_{Nnm}({\boldsymbol{q}})e^{im\omega _n({\boldsymbol{q}})t}{\boldsymbol{u}}_{Nn}({\boldsymbol{q}}).$$
(3)

Here q is a phonon wave vector in the first Brillouin zone (BZ), the coefficients aNnm(q) are determined by initial conditions, i2 = −1, ωn(q) is the angular frequency of the nth normal mode, and uNn(q) is a normalized vector in the direction of the oscillation for the Nth atom and the nth normal mode. The sum over m is introduced for compactness and aNn0(q) = 0. It will be apparent in the discussion of the toy model that it is simpler to work with the velocities (bNnm(q) ≡ imωn(q)):

$${\boldsymbol{v}}^{\rho N}(t) = \frac{{\mathrm{d}}}{{{\mathrm{d}}t}}{\boldsymbol{d}}^{\rho N}(t) = \mathop {\sum}\limits_n {\mathop {\sum}\limits_m {{\int} {\mathrm{d}} } } {\boldsymbol{q}}{\mathrm {e}}^{i{\boldsymbol{q}}\cdot {\boldsymbol{x}}_{0,\rho N}}b_{Nnm}({\boldsymbol{q}}){\mathrm {e}}^{im\omega _n({\boldsymbol{q}})t}{\boldsymbol{u}}_{Nn}({\boldsymbol{q}}).$$
(4)

Equation (4) implies that the velocities calculated from MD will contain information about all phonon wave vectors and normal modes. Taking the Fourier transform of these velocities directly will lead to complicated spectra, so it is best to select a single phonon wave vector and normal mode at a time. We can project onto a specific phonon wave vector Q:19,43

$${\boldsymbol{v}}_{\boldsymbol{Q}}^N(t) \equiv \mathop {\sum}\limits_\rho {{\boldsymbol{v}}^{\rho N}} (t){\mathrm {e}}^{ - i{\boldsymbol{Q}}\cdot {\boldsymbol{x}}_{0,\rho N}}.$$
(5)

Further, we can also project onto a phonon wave vector Q and project onto normal mode η, provided the normal modes can be calculated in advance:44

$$v_{\eta {\boldsymbol{Q}}}^N(t) \equiv {\boldsymbol{v}}_{\boldsymbol{Q}}^N(t)\cdot {\boldsymbol{u}}_{N\eta }({\boldsymbol{Q}}).$$
(6)

If the number of atoms per primitive cell is NA, then we expect to find 3NA normal modes after the unfolding procedure. There are two general approaches to generating spectra: (i) projected velocities based on Eq. (6) and (ii) velocity auto-correlation based on Eqs. (5) or (6).

Spectral lineshapes for projected velocities

For the projected velocities in Eq. (6) the spectrum is defined by

$$S_{\eta {\boldsymbol{Q}}}(\omega ) \equiv \hat N\hat R\mathop {\sum}\limits_{N = 1}^{N_A} {M_N^{1/2}} {\int}_{\!\!t_i}^{t_f} {\mathrm{d}} t{\boldsymbol{v}}_{\eta {\boldsymbol{Q}}}^N(t){\mathrm {e}}^{i\omega t}.$$
(7)

The operator \(\hat N\) normalizes the spectrum so that the maximum deviation from zero is one and the operator \(\hat R\) takes the real part. The limits of the integral (ti and tf) are the initial and final simulation times, respectively. The reason for not using the complex modulus will be apparent later. Contributions from all atoms are taken into account, and the weighting by \(M_N^{1/2}\) puts atoms with different masses on equal footing since \(b_{Nnm}({\boldsymbol{q}}) \propto M_N^{ - 1/2}\).42 In the idealized case when ti → −∞ and tf → ∞, the resulting spectrum is

$$\begin{array}{c}S_{\eta {\boldsymbol{Q}}}(\omega ) = \hat N2\pi \hat R\mathop {\sum}\limits_N {M_N^{1/2}} \mathop {\sum}\limits_\rho {\mathop {\sum}\limits_n {\mathop {\sum}\limits_m {\displaystyle\int {\mathrm{d}} } } } {\boldsymbol{q}}{\mathrm {e}}^{i({\boldsymbol{q}} - {\boldsymbol{Q}})\cdot {\boldsymbol{x}}_{0,\rho N}} \cdots \\ \cdots b_{Nnm}({\boldsymbol{q}})({\boldsymbol{u}}_{Nn}({\boldsymbol{q}})\cdot {\boldsymbol{u}}_{N\eta }({\boldsymbol{Q}}))\delta (\omega + m\omega _n({\boldsymbol{q}})).\end{array}$$
(8)

From this expression we see that the spectrum will exhibit peaks at the normal mode frequencies. It is also apparent that there will be no direct comparison to Raman spectra, since the spectrum features all of the normal mode frequencies, i.e. there are no selection rules. For an actual simulation of duration T with ti = −T/2 and tf = T/2:

$$\begin{array}{c}S_{\eta {\boldsymbol{Q}}}(\omega ) = \hat NT\hat R\mathop {\sum}\limits_N {M_N^{1/2}} \mathop {\sum}\limits_\rho {\mathop {\sum}\limits_n {\mathop {\sum}\limits_m {{\displaystyle\int} {\mathrm{d}} } } } {\boldsymbol{q}}{\mathrm {e}}^{i({\boldsymbol{q}} - {\boldsymbol{Q}})\cdot {\boldsymbol{x}}_{0,\rho N}} \cdots \\ \cdots b_{Nnm}({\boldsymbol{q}})({\boldsymbol{u}}_{Nn}({\boldsymbol{q}})\cdot {\boldsymbol{u}}_{N\eta }({\boldsymbol{Q}})){\mathrm{sinc}}[(\omega + m\omega _n({\boldsymbol{q}}))T/2].\end{array}$$
(9)

Here sinc(x) = sin(x)/x. If we make the approximation that both projections are complete, the appropriate fitting function is:

$$S_{\eta {\boldsymbol{Q}}}(\omega ) = \mathop {\sum}\limits_m {C_{\eta m}} ({\boldsymbol{Q}}){\mathrm{sinc}}[(\omega + m\omega _\eta ({\boldsymbol{Q}}))T/2].$$
(10)

The contributions of a given normal mode at negative and positive frequencies need to be considered in particular when very few configurations are used and the two peaks overlap significantly. The fitting parameters are the normal mode frequency ωη(Q) and the Cηm(Q) amplitudes (all real). If the complex modulus were used instead of the real part, the amplitudes would be complex numbers and an artificial non-linearity between the peaks at negative and positive frequencies would be created. This would increase the computational cost of the fitting procedure. More importantly, it appears that the analytical manipulations would be more difficult when using the complex modulus. It makes no difference if what we have defined as the spectrum has negative values, since the precise fitting function is derived. For a sufficiently isolated peak, the full-width (FW) of the central peak as defined by the zeros of the sinc function is only dependent on the simulation time:

$$\Delta \omega _{{\mathrm{FW}}} = \frac{{4\pi }}{T}.$$
(11)

So far the results are based on the harmonic approximation, which predicts no temperature dependence for the vibrational properties. Provided that the anharmonicity in the potential energy surface accessed by the atoms at a particular temperature is sufficiently small, we can assume that the superposition of normal modes given in Eq. (4) is still valid with the harmonic vibrational frequencies ωn(q) replaced by their temperature-dependent renormalized values. This is analogous to perturbing eigenvalues on a fixed basis set in quantum theory. In the context of quantum field theory the phonon self-energy (PSEn(q)) describes how the properties of a given phonon are modified due to its anharmonicity-induced interaction with all other phonons. The real part of PSEn(q) is the frequency shift Δωn(q) and the imaginary part is the half-linewidth Γn(q) (both up to a factor of ħ).13 The requirement for a weakly anharmonic system is that the magnitudes of the phonon self-energies are much smaller than the harmonic vibrational frequencies.24 The lifetime τn(q) is related to the linewidth 2Γn(q) by

$$\tau _n({\boldsymbol{q}}) = (2\Gamma _n({\boldsymbol{q}}))^{ - 1}.$$
(12)

In the method developed here, Γn(q) becomes a fitting parameter in the spectral lineshape, which we achieve by making the following substitution in Eq. (4):45,46

$${\mathrm {e}}^{im\omega _n({\boldsymbol{q}})t} \to {\mathrm {e}}^{[im\omega _n({\boldsymbol{q}}) - \Gamma _n({\boldsymbol{q}})]t}.$$
(13)

A similar analysis results in the fitting function for spectra including the effect of phonon lifetime:

$$\begin{array}{r}S_{\eta {\boldsymbol{Q}}}(\omega ) = \frac{2}{T}\mathop {\sum}\limits_m {\frac{1}{{(\omega + m\omega _\eta ({\boldsymbol{Q}}))^2 + (\frac{{2\Gamma _\eta ({\boldsymbol{Q}})}}{2})^2}}} \cdots \\ \cdots \left\{ {[C_{\eta m}({\boldsymbol{Q}})(\omega + m\omega _\eta ({\boldsymbol{Q}})) + C_{\eta m}^\prime ({\boldsymbol{Q}})\Gamma _\eta ({\boldsymbol{Q}})]{\mathrm{cosh}}(\Gamma _\eta ({\boldsymbol{Q}})T/2){\mathrm{sin}}[(\omega + m\omega _\eta ({\boldsymbol{Q}}))T/2] \cdots } \right.\\ \left. { \cdots + [C_{\eta m}({\boldsymbol{Q}})\Gamma _\eta ({\boldsymbol{Q}}) - C_{\eta m}^\prime ({\boldsymbol{Q}})(\omega + m\omega _\eta ({\boldsymbol{Q}}))]{\mathrm{sinh}}(\Gamma _\eta ({\boldsymbol{Q}})T/2){\mathrm{cos}}[(\omega + m\omega _\eta ({\boldsymbol{Q}}))T/2]} \right\}.\end{array}$$
(14)

Γη(Q) and the amplitude \(C_{\eta m}^\prime ({\boldsymbol{Q}})\) are additional fitting parameters (both also real). In the limit that Γη(Q) goes to zero, Eq. (14) reduces to Eq. (10) as expected.

Spectral lineshapes for velocity auto-correlation

In the previous section, it was assumed that the eigenvectors for the normal modes were already known. In cases where the system of interest is too large to determine eigenvectors, it is necessary to generate a full spectrum which contains information about all the normal modes. First, note that a full spectrum cannot be obtained from the magnitude of the velocity squared (see SI section I). However, this can be accomplished by manipulating the velocity auto-correlation as follows (τ is the temporal offset):

$$A_{\boldsymbol{Q}}(\tau ) \equiv \mathop {\sum}\limits_N {M_N} {\int} {\mathrm{d}} t{\boldsymbol{v}}_{\mathbf{Q}}^N(t)\cdot {\boldsymbol{v}}_{\boldsymbol{Q}}^N(t + \tau ).$$
(15)

Note that the phonon DOS can also be obtained from the velocity auto-correlation.27 In that case the sum over the unit cells (ρ) occurs after the velocity auto-correlation is calculated, whereas in Eq. (15) the sum over the unit cells precedes the velocity auto-correlation calculation via Eq. (5). Similar to the previous section the spectrum is given by

$$S_{\boldsymbol{Q}}^A(\omega ) \equiv \hat N\hat R{\int} {\mathrm{d}} \tau A_{\boldsymbol{Q}}(\tau ){\mathrm {e}}^{i\omega \tau }.$$
(16)

To obtain analytical expressions for the spectra it is useful to write:

$$\begin{array}{l}{\boldsymbol{v}}_{\boldsymbol{Q}}^N(t) = \mathop {\sum}\limits_{\rho \prime } {\mathop {\sum}\limits_{n\prime } {\mathop {\sum}\limits_{m\prime } {{\displaystyle\int} {\mathrm{d}} } } } {\boldsymbol{q}}^{\prime} {\mathrm {e}}^{i({\boldsymbol{q}}^{\prime} - {\boldsymbol{Q}})\cdot {\boldsymbol{x}}_{0,\rho ^{\prime} N}}b_{Nn^{\prime} m^{\prime} }({\boldsymbol{q}}\prime ){\mathrm {e}}^{im^{\prime} \omega _{n^{\prime} }({\boldsymbol{q}}^{\prime} )t}{\boldsymbol{u}}_{Nn\prime }({\boldsymbol{q}}\prime ),\\ {\boldsymbol{v}}_{\boldsymbol{Q}}^N(t + \tau ) = \mathop {\sum}\limits_\rho {\mathop {\sum}\limits_n {\mathop {\sum}\limits_m {{\displaystyle\int} {\mathrm{d}} } } } {\boldsymbol{q}}{\mathrm {e}}^{i({\boldsymbol{q}} - {\boldsymbol{Q}})\cdot {\boldsymbol{x}}_{0,\rho N}}b_{Nnm}({\boldsymbol{q}}){\mathrm {e}}^{im\omega _n({\boldsymbol{q}})(t + \tau )}{\boldsymbol{u}}_{Nn}({\boldsymbol{q}}).\end{array}$$

The resulting spectrum is

$$\begin{array}{c}S_{\boldsymbol{Q}}^A(\omega ) = \hat N\hat R\mathop {\sum}\limits_N {M_N} \mathop {\sum}\limits_{\rho ,\rho ^{\prime} } {\mathop {\sum}\limits_{n,n^{\prime} } {\mathop {\sum}\limits_{m,m^{\prime} } {{\int} {\mathrm{d}} } } } {\boldsymbol{q}}{\int} {\mathrm{d}} {\boldsymbol{q}}^{\prime} {\mathrm {e}}^{i[({\boldsymbol{q}} - {\boldsymbol{Q}})\cdot {\boldsymbol{x}}_{0,\rho N} + ({\boldsymbol{q}}^{\prime} - {\boldsymbol{Q}})\cdot {\boldsymbol{x}}_{0,\rho ^{\prime} N}]} \cdots \\ \cdots b_{Nnm}({\boldsymbol{q}})b_{Nn^{\prime} m^{\prime} }({\boldsymbol{q}}^{\prime} )({\boldsymbol{u}}_{Nn}({\boldsymbol{q}})\cdot {\boldsymbol{u}}_{Nn{^{\prime}}}({\boldsymbol{q}}^{\prime} ))I.\end{array}$$
(17)

The integral I is defined as (the limits are not specified here):

$$I \equiv {\int} {\mathrm{d}} \tau {\mathrm {e}}^{i(\omega + m\omega _n({\boldsymbol{q}}))\tau }{\int} {\mathrm{d}} t{\mathrm {e}}^{i(m\omega _n({\boldsymbol{q}}) + m^{\prime} \omega _{n^{\prime} }({\boldsymbol{q}}^{\prime} ))t}.$$
(18)

In the idealized case, all of the integration limits are infinite:

$$\begin{array}{l}I = ({\int}_{ - \infty }^\infty {\mathrm{d}} t{\mathrm {e}}^{i(m\omega _n({\boldsymbol{q}}) + m\prime \omega _{n^{\prime} }({\boldsymbol{q}}\prime ))t})({\int}_{ - \infty }^\infty {\mathrm{d}} \tau {\mathrm {e}}^{i(\omega + m\omega _n({\boldsymbol{q}}))\tau }),\\ I = (2\pi )^2\delta (m\omega _n({\boldsymbol{q}}) + m^{\prime} \omega _{n^{\prime} }({\boldsymbol{q}}^{\prime} ))\delta (\omega + m\omega _n({\boldsymbol{q}})).\end{array}$$
(19)

The spectrum will then contain peaks at the normal mode frequencies, in similarity to the idealized result for the projected velocities (Eq. (8)). This is also similar to the result derived previously for the velocity auto-correlation.47 For a finite simulation, the integrals are no longer independent:

$$I = {\int}_{\!\!0}^T {\mathrm{d}} \tau {\mathrm {e}}^{i(\omega + m\omega _n({\boldsymbol{q}}))\tau }{\int}_0^{T - \tau } {\mathrm{d}} t{\mathrm {e}}^{i(m\omega _n({\boldsymbol{q}}) + m\prime \omega _{n\prime }({\boldsymbol{q}}\prime ))t}.$$
(20)

Note that special consideration is required in evaluating this integral when n(q) + mωn(q′) = 0. Assuming the Q projection is complete, the resulting fitting function is (\(C_{nm}^{(\lambda )}({\boldsymbol{Q}})\) for λ = 1, 2, 3, 4 are real fitting amplitudes):

$$\begin{array}{r}S_{\boldsymbol{Q}}^A(\omega ) = \mathop {\sum}\limits_n {\mathop {\sum}\limits_m {\left\{ {C_{nm}^{(1)}({\boldsymbol{Q}}){\mathrm{sinc}}[(\omega + m\omega _n({\boldsymbol{Q}}))T] \cdots } \right.} } \\ \cdots + C_{nm}^{(2)}({\boldsymbol{Q}}){\mathrm{sin}}[(\omega + m\omega _n({\boldsymbol{Q}}))T/2]{\mathrm{sinc}}[(\omega + m\omega _n({\boldsymbol{Q}}))T/2] \cdots \\ \cdots + C_{nm}^{(3)}({\boldsymbol{Q}}){\mathrm{sinc}}^2[(\omega + m\omega _n({\boldsymbol{Q}}))T/2] \cdots \\ \left. { \cdots + C_{nm}^{(4)}({\boldsymbol{Q}})\frac{{\partial \{ {\mathrm{sinc}}[(\omega + m\omega _n({\boldsymbol{Q}}))T]\} }}{{\partial [(\omega + m\omega _n({\boldsymbol{Q}}))T]}}} \right\}.\end{array}$$
(21)

Equation (20) can be extended to include the effect of lifetime based on Eq. (13):

$$I = {\int}_0^T {\mathrm{d}} \tau {\mathrm {e}}^{[i(\omega + m\omega _n({\boldsymbol{q}})) - \Gamma _n({\boldsymbol{q}})]\tau }{\int}_0^{T - \tau } {\mathrm{d}} t{\mathrm {e}}^{[i(m\omega _n({\boldsymbol{q}}) + m^{\prime} \omega _{n^{\prime} }({\boldsymbol{q}}^{\prime} )) - (\Gamma _n({\boldsymbol{q}}) + \Gamma _{n^{\prime}}({\boldsymbol{q}}^{\prime} ))]t}.$$
(22)

The appropriate fitting function can then be shown to be (Cnml(Q) and \(C_{nml}^\prime ({\boldsymbol{Q}})\) are real fitting amplitudes, and l = ±1):

$$\begin{array}{r}S_{\boldsymbol{Q}}^A(\omega ) = \frac{1}{T}\mathop {\sum}\limits_n {\mathop {\sum}\limits_m {\mathop {\sum}\limits_l {\frac{1}{{(\omega + m\omega _n({\boldsymbol{Q}}))^2 + (\frac{{2\Gamma _n({\boldsymbol{Q}})}}{2})^2}}} } } \cdots \\ \cdots \left\{ {[C_{nml}({\boldsymbol{Q}})(\omega + m\omega _n({\boldsymbol{Q}})) + C_{nml}^\prime ({\boldsymbol{Q}})l\Gamma _n({\mathbf{Q}})]{\mathrm {e}}^{l\Gamma _n({\boldsymbol{Q}})T}{\mathrm{sin}}[(\omega + m\omega _n({\boldsymbol{Q}}))T] \cdots } \right.\\ \left. { \cdots + [C_{nml}({\boldsymbol{Q}})l\Gamma _n({\boldsymbol{Q}}) - C_{nml}^\prime ({\boldsymbol{Q}})(\omega + m\omega _\eta ({\boldsymbol{Q}}))]\{ {\mathrm {e}}^{l\Gamma _n({\boldsymbol{Q}})T}{\mathrm{cos}}[(\omega + m\omega _n({\boldsymbol{Q}}))T] - 1\} } \right\}.\end{array}$$
(23)

Note that \(m\omega _n({\boldsymbol{q}}) + m\prime \omega _{n\prime }({\boldsymbol{q}}\prime ) = 0\) does not cause any problem when the effect of lifetime is included. When \({\mathrm{\Gamma }}_n({\boldsymbol{Q}}) = 0\) in Eq. (23), only the \(C_{nm}^{(1)}({\boldsymbol{Q}})\) and \(C_{nm}^{(2)}({\boldsymbol{Q}})\) terms in Eq. (21) will be recovered up to a constant factor. Equations (10), (14), (21) and (23) represent the key results of this paper. To summarize, the assumptions that are required to derive these lineshapes are quite simple: (i) Eq. (4) remains valid with the harmonic oscillation frequencies replaced by their renormalized values, (ii) the phonon lifetimes are included using a damped simple harmonic oscillator model (Eq. (13)), and (iii) the projections in Eqs. (5) and (6) perfectly select the desired phonon wave vector and normal mode. These assumptions are also required to obtain the standard fitting function used in the literature with the additional requirement that the simulation time be infinite (more discussion in the next section).

Comments on previous work

In the literature the standard procedure is to calculate spectra with the fast Fourier transform (FFT) and then fit with Lorentzian functions, however many configurations are required to obtain an accurate result.25 Furthermore, a reasonable question is: Why has the sinc-like lineshape not been resolved for spectra resulting from low-temperature (nearly harmonic) simulations? The reason is that the grid spacing for the FFT is too coarse:

$$\Delta \omega _{{\mathrm{FFT}}} = \frac{{2\pi }}{T} = \frac{1}{2}\Delta \omega _{{\mathrm{FW}}}.$$
(24)

That is, on average the FFT places two points on the central peak and one point on the side peaks of the sinc function. Another question is: Can the Lorentzian lineshape be rigorously justified? It is useful to see what happens when the effect of lifetime is included in the limit of an infinite simulation. This is most easily demonstrated for the projected velocities. To keep the real part of the exponential involving the half-linewidth finite, the time limits are changed to the asymmetric form ti = 0 and tf = T, and only the peak at positive frequency is considered. The integral to evaluate is thus:

$$\mathop {{\lim }}\limits_{T \to \infty } {\int}_{\!\!0}^T {\mathrm{d}} t{\mathrm {e}}^{[i(\omega - \omega _n({\boldsymbol{q}})) - \Gamma _n({\boldsymbol{q}})]t} = - \frac{1}{{i(\omega - \omega _n({\boldsymbol{q}})) - \Gamma _n({\boldsymbol{q}})}}.$$

If the complex modulus squared is used instead of the real part, then the spectral lineshape is

$$S_{\eta {\boldsymbol{Q}}}(\omega ) \propto \frac{1}{{(\omega - \omega _\eta ({\boldsymbol{Q}}))^2 + (\frac{{2\Gamma _\eta ({\boldsymbol{Q}})}}{2})^2}}.$$
(25)

This is indeed a Lorentzian peaked at ω = ωη(Q) with ΔωFWHM = 2Γη(Q). A similar result holds for the velocity auto-correlation approach. Therefore, the fitting function of the standard procedure is only valid for long simulation times. Also, part of the width of the Lorentzian will be due to the finite simulation time, which can be estimated by making a comparison between the Lorentzian and sinc functions using their limiting behavior to the Dirac delta function:

$$\delta (x) = \mathop {{\lim }}\limits_{\epsilon \to 0} \frac{1}{\pi }\frac{\epsilon }{{x^2 + \epsilon ^2}} = \mathop {{\lim }}\limits_{\epsilon \to 0} \frac{1}{{\pi \epsilon }}{\mathrm{sinc}}\left(\frac{x}{\epsilon }\right).$$
(26)

The FWHM of the Lorentzian (first term) is

$${\mathrm{FWHM}}_{{\mathrm{Lorentzian}}} = 2\epsilon .$$

The full-width of the central peak as defined by the zeros of sinc (second term) is

$${\mathrm{FW}}_{{\mathrm{sinc}}} = 2\pi \epsilon .$$

Combining these results implies:

$${\mathrm{FWHM}}_{{\mathrm{Lorentzian}}} = \frac{{{\mathrm{FW}}_{{\mathrm{sinc}}}}}{\pi }.$$
(27)

Next, using Eq. (11) we obtain an important result with a quantitative approximation for the contribution to the FWHM of a Lorentzian due to the finite simulation time:

$$\Delta \omega _{{\mathrm{FWHM}}} \approx \frac{4}{T}.$$
(28)

Illustration for a Toy model

We now illustrate the formalism developed in this paper (specifically, Eqs. (10), (14), (21) and (23)) on a toy model, see Fig. 1a, which consists of a classical particle in one dimension subject to a force given by

$$F(x,v) = - k_1x - k_2x^2 - k_3x^3 - cv.$$
(29)

Here x is the position, v is the velocity, the k’s are spring constants, and c is a damping constant. Without the damping term this corresponds to a potential:

$$V(x) = \frac{1}{2}k_1x^2 + \frac{1}{3}k_2x^3 + \frac{1}{4}k_3x^4.$$
(30)

Explicitly, the acceleration a and jerk j are (M is the mass of the particle):

$$a = - \frac{1}{M}(k_1x + k_2x^2 + k_3x^3 + cv),$$
(31)
$$j = \frac{{{\mathrm{d}}a}}{{{\mathrm{d}}t}} = - \frac{1}{M}[(k_1 + 2k_2x + 3k_3x^2)v + ca].$$
(32)

The system is propagated on the basis of a truncated Taylor expansion (Δt is the time-step and n = 1, 2, …):

$$x(t_{n + 1}) = x(t_n) + v(t_n)\Delta t + \frac{1}{2}a(t_n)(\Delta t)^2 + \frac{1}{6}j(t_n)(\Delta t)^3,$$
(33)
$$v(t_{n + 1}) = v(t_n) + a(t_n)\Delta t + \frac{1}{2}j(t_n)(\Delta t)^2.$$
(34)

The total energy of the system is monitored and compared to its initial value to validate the choice of time-step (n ≥ 2):

$$\begin{array}{c}E(t_n) = \frac{1}{2}M(v(t_n))^2 + V(x(t_n))\cdots\\ \cdots+ \sum\limits_{n^{\prime} = 2}^n c (\frac{{v(t_{n^{\prime} - 1}) + v(t_{n{^{\prime}}})}}{2})(x(t_{n^{\prime} }) - x(t_{n^{\prime} - 1})).\end{array}$$
(35)

In order, these terms correspond to the kinetic energy, potential energy, and the cumulative drag energy. The simulation parameters are reported in Table 1. The mass is equivalent to that of a hydrogen atom, the total initial energy is the thermal energy kBT at room temperature, and the time-step is the same as that used in our DFT calculations. For reference, for silicon and black phosphorus the vibrational frequencies are less than about 600 cm−1, and the frequency shifts and linewidths are no more than about 6 cm−1 at room temperature.25,48 The initial position is the amplitude for x < 0 and is found by solving V(x) = E0, with E0 = E(t1) as the specified initial energy. Therefore, the initial velocity is set to zero. Without the jerk term the percent deviation of the total energy from its initial value is as much as ~103%; however, the inclusion of the jerk term reduced this to ~10−3%.

Fig. 1
figure 1

Three systems considered as examples of the developed method. a Sketch of the toy model where a classical particle of mass M in one dimension is subjected to an external force. x, ki and c correspond to the position of the particle, spring constants, and damping constant, respectively. Top view of b graphene and c single-layer hexagonal boron nitride (hBN)59

Table 1 Simulation parameters for the toy model

Combinations of harmonic/anharmonic and no damping/damping are considered for the derived projection and auto-correlation lineshapes, as well as the standard Lorentzians: (i) harmonic, no damping; (ii) harmonic, damping; and (iii) anharmonic, no damping. The anharmonic with damping case is not considered, since no general analytical solution exists to the equation of motion. The Lorentzians are fit to the modulus squared of the FFT (|FFT|2), while the derived lineshapes are fit to the real part of the Fourier transform (Re{FT}) calculated on a fine non-uniform frequency grid. To construct the grid, the peak positions are estimated using the relative maxima of |FFT|2 and then a grid spacing of ΔωFW/50 is set for ±5ΔωFW about these estimated positions (see Eq. (11)). Outside of these regions the grid spacing is increased by a factor of five. The convergence of the fitted vibrational frequencies and Γ values (damping only) with the number of configurations is considered. The configuration grid consists of values ranging from 11 to 10,001 (increment of 10) and 10,001 to 20,001 (increment of 100). The only exception is for harmonic with damping (auto-correlation only) for which the values range from 11 to 10,001 (increment of 10), 10,001–100,001 (increment of 100), and 100,001–200,001 (increment of 500). Only odd values are used to enable symmetry about zero time when necessary. The parameters extracted from the fitting procedure are averaged over 30 consecutive configuration grid points to prevent outliers from defining the convergence. Therefore, on the horizontal axis of the plots to follow in the main text and SI maximum number of configurations refers to the largest of the configuration grid values in a given average. Vibrational frequencies and Γ values are considered converged when the current and all subsequent averaged values are within the exact numbers by at least 0.5 cm−1. Note that in the discussion and figures that follow frequency refers to the vibrational frequency ν, not the angular frequency ω.

Toy model: harmonic, no damping

Harmonic with no damping means that k2, k3, and c are zero. The exact analytical solution takes the form:

$$v(t) = v_ + {\mathrm {e}}^{i\omega _{\mathrm{H}}t} + v_ - {\mathrm {e}}^{ - i\omega _{\mathrm{H}}t}.$$
(36)

This is a special case of Eq. (4) with one mode along one dimension, which implies that the derived harmonic with no damping lineshapes will apply here. The frequency convergence and sample spectra for the projection lineshape are illustrated in Figs. S1 and S2, respectively. The corresponding information for the auto-correlation lineshape is shown in Figs. S3 and S4. Our analysis indicates that using the proposed method reduces the required number of configurations by factors of about 10 (projection) and 29 (auto-correlation) for the vibrational frequency.

Toy model: harmonic, damping

Harmonic with damping means that k2 and k3 are set to zero. The exact analytical solution takes the form:

$$v(t) = v_ + {\mathrm {e}}^{(i\omega _{\mathrm{D}} - \Gamma )t} + v_ - {\mathrm {e}}^{( - i\omega _{\mathrm{D}} - \Gamma )t}.$$
(37)

This is a special case of Eq. (4) with the inclusion of Eq. (13), which implies that the derived harmonic with damping lineshapes will apply here. The fitting procedure will extract a modified frequency:

$$\nu _{\mathrm{D}} = \omega _{\mathrm{D}}/(2\pi ) = \sqrt {\omega _{\mathrm{H}}^2 - \Gamma ^2} /(2\pi ) \approx 99.9997\,{\mathrm{cm}}^{ - 1}.$$

Since ωH Γ, νD ≈ νH. The frequency/damping constant convergence and sample spectra for the projection lineshape are illustrated in Figs. S5 and S6, respectively. The corresponding information for the auto-correlation lineshape is shown in Figs. S7 and S8. Using the new method reduces the required number of configurations by factors of about 10 (projection) and 26 (auto-correlation) for the vibrational frequency, and factors of about 12 (projection) and 99 (auto-correlation) for Γ.

Toy model: anharmonic, no damping

Anharmonic with no damping means that c is set to zero. The exact frequency for determining the convergence can be calculated starting from energy conservation:

$$\frac{1}{2}M\left(\frac{{{\mathrm{d}}x}}{{{\mathrm{d}}t}}\right)^2 + V(x) = E_0,$$
$$\nu _{\mathrm{A}} = \left(\sqrt {2M} {\int}_{{\!\!A}_ - }^{A_ + } {\frac{{{\mathrm{d}}x}}{{\sqrt {E_0 - V(x)} }}}\right )^{ - 1} \approx 86.8\,{\mathrm{cm}}^{ - 1}.$$
(38)

The integration limits A and A+ are the amplitudes at negative and positive positions, respectively, as determined by solving V(x) = E0. The anharmonicity causes a significant downshift in frequency (100 cm−1 → 86.8 cm−1). Since the motion remains periodic, the solution for the position can be written as a Fourier series (ωA = 2πνA):

$$x(t) = \mathop {\sum}\limits_{m = - \infty }^\infty {a_m} {\mathrm {e}}^{im\omega _{\mathrm{A}}t}.$$
(39)

The constant term for m = 0 in the Fourier series will not be zero here, since there is an odd term in the potential (k3 ≠ 0). Taking a Fourier transform will then place a peak at zero as well as the vibrational frequency (and its other integral multiples), which complicates the analysis. Taking the derivative to get the velocity removes this effect and suggests why the velocities should be used instead of the displacements from equilibrium. The frequency convergence and sample spectra for the projection lineshape are illustrated in Figs. S9 and S10, respectively. The corresponding information for the auto-correlation lineshape is shown in Fig. S11. No additional interesting information is obtained from the spectra. Using the new method reduces the required number of configurations by factors of about 8 (projection) and 10 (auto-correlation) for the vibrational frequency.

Comment on toy model

For a given set of toy model parameters, the projection and auto-correlation approaches for the developed formalism converge in a similar number of configurations. However, for the standard procedure based on |FFT|2 and a Lorentzian fit the auto-correlation technique takes much longer to converge than the projection method. A possible explanation is that while the Lorentzian is the exact lineshape for both methodologies in the infinite simulation time limit, it is a worse approximation to the exact auto-correlation lineshape for finite simulation time as compared to the corresponding projection case.

Application to real materials

With the developed theory illustrated for a toy model, we now show its application to three materials: graphene, single-layer hBN (see Fig. 1b, c), and silicon. Ab initio molecular dynamics (AIMD) simulations were carried out for a supercell of each material and more than 100,000 velocity configurations were obtained (see “Methods” section for computational details). The effective velocities are calculated using Eq. (5) with projection at Q = 0 and the velocity auto-correlation is evaluated with Eq. (15). For illustration purposes only the BZ center is analyzed, which is relevant to single-phonon Raman spectroscopy.9,49,50 For graphene and hBN only the in-plane components of the velocity auto-correlation are used to select the doubly degenerate in-plane optical modes at Q = 0. For silicon the analysis is simple since all optical modes are triply degenerate at Q = 0. The spectra are calculated using both the standard |FFT|2 approach and Eq. (16), and are correspondingly fit using Eqs. (25) and (23). To construct the non-uniform frequency grid for the new method, the peak positions are estimated using the relative maxima of |FFT|2 and then a grid spacing of ΔωFW/50 is set for ±5ΔωFW about these estimated positions (see Eq. (11)). The values of the configuration grid are from 11 to 10,001 (increment of 10) and 10,001 to 100,001 or 125,001 (increment of 100), and the convergence tolerance is still set to ±0.5 cm−1.

Graphene

The vibrational frequency and Γ convergence rates are illustrated in Fig. 2 and sample spectra are shown in Fig. 3 for single-layer graphene at a temperature of 50 K. Details are provided in the captions of the figures.

Fig. 2
figure 2

Convergence properties for graphene treated using ab initio molecular dynamics with a temperature of around 50 K and the velocity auto-correlation method. The E2g symmetry modes are isolated at the Γ-point by projections. The convergence is studied using two approaches, which are denoted by “|FFT|2 + Lorentzian” and “Re{FT} + Analytical”. In both cases, the effective velocities at the Γ-point were obtained using Eq. (5) and the velocity auto-correlation was subsequently calculated with Eq. (15). In the “|FFT|2 + Lorentzian” method, which is the standard procedure and assumes infinite simulation time, the spectrum is calculated as the modulus squared of the fast Fourier transform of the velocity auto-correlation and the resulting peaks are subsequently fit using Lorentzian functions (Eq. (25)). In the “Re{FT} + Analytical” method, which is the new approach and is valid at all simulation times based on our analytical derivations, the spectra are evaluated with Eq. (16) and the vibrational frequencies and half-linewidths are extracted by fitting with Eq. (23). a Convergence of the vibrational frequency with the number of configurations (=NC). Convergence is indicated by the vertical lines with NC = 14,501 for the Lorentzian fit and NC = 2,461 for the analytical model, with improvement by a factor of about 5.9. b Convergence of the half-linewidth Γ with the number of configurations. Convergence is indicated by the vertical lines with NC = 27,301 for the Lorentzian fit and NC = 12,301 for the analytical model, with improvement by a factor of about 2.2. The developed method results in a significant reduction in the required simulation time to obtain converged properties

Fig. 3
figure 3

Sample spectra obtained from ab initio molecular dynamics applied to graphene (based on Fig. 2). a and b 10,001 configurations. c and d 50,001 configurations. e and f 100,001 configurations. The derived fitting function only deviates from the calculated spectrum slightly even after 100,001 configurations

Classical dynamics based on parameterized force-fields can be used as an alternative to AIMD, especially for very large systems.21 In this case, it is easier to use larger supercells to include many more phonon wave vectors. To that end we carried out MD on a (20 × 20) supercell of graphene at around 50 K using four common potentials (see “Methods” section for computational details) and plotted the spectra after 50,001 time steps (see Fig. 4). The fitted vibrational frequencies are close to the values obtained in a previous comprehensive study using very long simulation times and a detailed statistical analysis.51 In general the analytical fitting function does not match the calculated spectra very well. If the potentials were producing accurate trajectories, the spectra in Fig. 4 would all match the spectrum in Fig. 3d obtained from AIMD for which the analytical fit exactly matches the calculated spectrum. This is clearly not the case, so we can conclude that the parameterized force-fields produce inaccurate trajectories for oscillations close to equilibrium despite their outstanding utility for other applications. Slightly inaccurate trajectories could also explain the asymmetries for the AIMD spectra in Fig. 3f, although to a significantly smaller extent. In general the spectra in Fig. 4 are much broader than expected, which the analytical fitting function superficially attempts to recreate by using large Γ values. This observation is consistent with previous work, where the half-linewidth values predicted for silicon with the Tersoff potential disagreed with those obtained from AIMD.25

Fig. 4
figure 4

Sample spectra obtained for graphene from parameterized force-field molecular dynamics at 50 K and the velocity auto-correlation method. The E2g symmetry modes are isolated at the Γ-point by projections. All spectra are calculated using 50,001 configurations. a REBO60, b AIREBO-M61, c LCBOP62, and d Tersoff63. Together, the deviation of the spectral lineshapes from that of the corresponding AIMD result (Fig. 3d) and the inability of the analytical model to fit the data suggest that the potentials create inaccurate trajectories for oscillations close to equilibrium

Hexagonal boron nitride

The vibrational frequency and Γ convergence rates are illustrated in Fig. 5. Details are provided in the caption of the figure. Sample spectra are not shown since they are very similar to those of graphene (Fig. 3).

Fig. 5
figure 5

Convergence properties for hexagonal boron nitride treated using ab initio molecular dynamics with a temperature of around 50 K and the velocity auto-correlation method. The E2g symmetry modes are isolated at the Γ-point by projections. a Convergence of the vibrational frequency with the number of configurations (=NC). Convergence is indicated by the vertical lines with NC = 16,801 for the Lorentzian fit and NC = 3561 for the analytical model, with improvement by a factor of about 4.7. b Convergence of the half-linewidth Γ with the number of configurations. Convergence is indicated by the vertical lines with NC = 52,301 for the Lorentzian fit and NC = 5851 for the analytical model, with improvement by a factor of about 8.9. The developed method results in a significant reduction in the required simulation time to obtain converged properties

Silicon

The vibrational frequency and Γ convergence are illustrated in Fig. S12 and sample spectra are shown in Fig. S13 for silicon at a temperature of 50 K. There is almost no difference between the two methods for the convergence of the vibrational frequency; however, the convergence of the half-linewidth is improved by a factor of about 1.4 by using the new approach. Corresponding results are shown in Figs. 6 and 7 for room temperature (300 K) where the anharmonicity is significantly stronger.

Fig. 6
figure 6

Convergence properties for silicon treated using ab initio molecular dynamics with a temperature of around 300 K and the velocity auto-correlation method. The triply degenerate optical modes are isolated at the Γ-point by projections. a Convergence of the vibrational frequency with the number of configurations (=NC), which is indicated by the vertical lines (almost no difference between the methods). b Convergence of the half-linewidth Γ with the number of configurations is not obtained due to the significant distortion of the spectral lineshapes from their theoretically predicted forms (see next figure)

Fig. 7
figure 7

Sample spectra obtained from ab initio molecular dynamics applied to silicon at 300 K (based on Fig. 6). a and b 10,001 configurations. c and d 50,001 configurations. e and f 100,001 configurations. At high-temperature the anharmonicity causes strong deviations from the theoretically predicted lineshape. Therefore, the extracted half-linewidth values obtained from either method cannot be expected to be accurate

Discussion

Up until this point we have demonstrated that the developed method improves the convergence of vibrational frequencies and lifetimes with respect to the total simulation time. There are three simple points that explain this observation: (i) the derived lineshapes (Eqs. (10), (14), (21) and (23)) explicitly include the simulation time, (ii) the Lorentzian function used in the standard fitting procedure is only exact in the infinite simulation limit, and (iii) the FFT poorly samples each peak, only placing a few points on the dominant part of the peak, whereas in the developed method the Fourier transform is calculated on a much finer grid. The first two points are more important for the lifetime, since the simulation time and the lifetime both contribute to the width of a peak. The two methods can be evaluated based on more than the number of configurations to reach convergence. In general, the new method produces smooth curves, while the standard procedure results in rapid oscillations (this is also true for the toy model, see SI).

The developed method performs better for the toy model than the real materials, which can be explained by several points. Since there is only one vibrational mode and one phonon wavevector in the toy model, the system is initially thermalized, i.e. all phonon modes follow the equilibrium distribution. However, in the case of graphene, hBN, and silicon there are multiple modes and phonon wavevectors included, so the random initialization of the velocities does not result in the equilibrium distribution immediately. Since the anharmonicity is weak, the transfer of energy between the phonon modes to reach the equilibrium distribution takes some critical time which is related to the phonon lifetimes. Furthermore, the distortion of the peaks in Figs. 3f, S13f, and 7f indicates at least a partial breakdown of the initial assumptions. These are to be contrasted with the corresponding toy model spectra based on the auto-correlation in the right column of Fig. S8, which are symmetric and have the fitting function matching the simulation data exactly. The anharmonicity-induced coupling between phonons partially invalidates the simplicity of the superposition in Eq. (4). Lastly, the inclusion of the phonon lifetime by a simple exponential decay (Eq. (13)) represents a continuous dissipation of energy from the mode during the simulation. This differs from the actual thermalization process, which is a continuous exchange of energy between different modes due to their interactions, with the total energy being constant. This simple assumption is therefore perhaps only valid over short timescales, which may be difficult to treat due to the simulation broadening.

Although silicon is often featured in elementary textbooks as a standard example, it is evident that the phonon anharmonicity is strong even at 50 K (Fig. S13) in comparing to the spectra from graphene (Fig. 3) and hBN, which may in part be due to the difference in dimensionality. More specifically, the 50 K silicon spectrum broadens more than the fitting lineshape anticipates as the simulation evolves and a satellite peak at higher frequency than the central peak emerges, both of which contribute to the fact that the two fitting methods converge similarly for the vibrational frequency. At 300 K the phonon–phonon coupling is much stronger, causing the fitting functions to strongly deviate from the raw spectrum data. Therefore, the lifetime cannot be extracted reliably (Fig. 7). It is important to emphasize that when it is not possible to obtain accurate results from the developed method, the standard procedure will not produce meaningful results either since the underlying physical assumptions are the same except that the latter further imposes that the simulation time be infinite. Future work can address these issues, although the improvement with our new method is already striking for graphene and hBN.

To summarize, a new approach to understand vibrational lineshapes was derived for spectra calculated from MD and include the effect of anharmonicity-induced frequency shift and lifetime, as well as the simulation broadening (Eqs. (10), (14), (21) and (23)). For a toy model we demonstrated at least an order of magnitude reduction in the number of simulation steps to obtain converged vibrational properties in nearly all cases considered as compared to the standard extraction procedure. The theory was also illustrated for graphene, hBN, and silicon with the dynamics carried out at the DFT level. For graphene the vibrational frequency and lifetime convergence were improved by factors of 5.9 and 2.2, and for hBN the corresponding factors were 4.7 and 8.9 (both materials at 50 K). For silicon the anharmonicity results in spectral lineshapes that deviate from the theoretically predicted ones, which is more extreme at higher-temperature. At 50 K the lifetime can be converged using both methods, with the new method reducing the required simulation time by a factor of about 1.4. However, at 300 K the deviation of the calculated spectra from the theoretically predicted lineshape is too extreme and claiming accurate values for the lifetimes from either method is not possible. Further analytical work will be required to accurately describe cases where the anharmonicity is moderate. Considering the simplicity of the initial assumptions, only resolving the spectrum on a fine grid and including the simulation time in the lineshapes has significantly reduced the computational effort to obtain temperature-dependent phonon properties from AIMD in the regime where the anharmonicity is sufficiently weak and results in well-defined renormalized phonon quasiparticles. In all cases, adoption of the developed method is encouraged, since, by design, it will never be outperformed by the conventional method. Our extension of signal analysis to material vibrations represents a state-of-the-art advance in calculating temperature-dependent phonon properties. Application of the proposed approach has the potential to be far reaching as the theory applies equally well to ab initio MD based on DFT, time-dependent DFT dynamics, and parameterized force-fields and is thus expected to have important impact on topics ranging from strongly correlated materials with sophisticated treatment of electron–electron interactions to biological systems (large systems precluding long simulation times).

Methods

VASP was used for plane-wave DFT calculations52,53 with the PBE exchange correlation functional.54 The projector augmented wave pseudopotentials55 were used with an energy cutoff of 500 eV. The tolerance for the electronic loop was set to 10−5 eV. Single-layer graphene and hBN (4 × 4) supercells with a 15 Å vacuum region and a silicon (3 × 3 × 3) supercell were each treated using a Γ-point k-point sampling and a Gaussian smearing of 0.01 eV. The lattice and ions were first relaxed to a force cutoff of 0.01 eV/Å. The ab initio MD56,57 were carried out with a time-step of 0.5 fs. No thermostat was used to prevent systematic deviations from the desired Newtonian trajectories implicit in the initial assumptions for the analytical derivations. The dynamics were initialized with all atoms starting at their equilibrium positions and with velocities randomly assigned so that the initial temperature was 100 K (600 K), which quickly stabilizes to an average temperature of around 50 K (300 K). No time-steps were discarded and more than 100,000 velocity configurations were obtained.

The MD with parameterized force-fields were carried out using LAMMPS.58 A (20 × 20) supercell of graphene was used, and the lattice and ions were first relaxed to a force cutoff of 0.01 eV/Å. The dynamics were carried out with a time-step of 0.05 fs only saving the trajectory every tenth step. The smaller time-step is necessary since VASP uses a fourth-order predictor-corrector method, while LAMMPS uses a velocity-Verlet integrator. No thermostat was used and the dynamics were initialized in the same manner as VASP so that an average temperature of 50 K was maintained (although the fluctuations were smaller due to the larger supercell). No time-steps were discarded and more than 50,000 velocity configurations were obtained.