Introduction

Computational studies of material properties, particularly transport phenomena of heat, charge, and mass at micro- and nanoscales1, have recently provided accurate quantification of measured observables and fundamental understanding of complex physics, due in large part to increased computational power, theoretical advances, and availability of numerical tools2,3,4,5,6. First-principles calculations based on density functional theory (DFT) have demonstrated remarkable power in accurately predicting thermal and electrical properties of semiconductors7,8,9,10, for example, the ultrahigh thermal conductivity of cubic BAs and BN11,12,13,14, phonon hydrodynamics in graphitic materials15,16,17, and high mobility of BAs18. Many of these successful simulations have been based on materials with simple structures and compositions. More recently, researchers have been exploring more complex materials with correspondingly diverse and novel properties including large thermal resistance19,20, anomalous thermal Hall effect21,22,23, and topological superconductivity24,25,26, though this requires more extensive computational power due to much larger degrees of freedom.

Calculations and simulations of crystalline material properties are typically based on periodically arranged unit cells that contain the smallest number of degrees of freedom, i.e., the primitive unit cell. The most widely used primitive cells, for which basic properties can be described, are called Wigner-Seitz cells27 (WSC; real space) and their corresponding first Brillouin zones (FBZ; reciprocal space). These are typically used for calculations of lattice vibrational and transport properties because they have the cheapest computational cost. However, such calculations, particularly in complex crystals, often have awkward geometries with WSC and FBZ having irregular, non-intuitive shapes in Cartesian space. For example, bismuth telluride (Bi2Te3, space group R-3m) has a primitive cell with an angle of 24° between pairs of lattice vectors. This adds difficulty in analyzing properties along natural high-symmetry lines (in-plane and cross-plane) as these require projection from the oblique lattice vectors. On the other hand, there are alternative unit cells with different geometries that can fill space and respect natural symmetries. Experiments, in particular, tend to probe and analyze material properties from these more convenient conventional geometries where the symmetries of the lattice are obvious and straightforward. For example, the inelastic neutron scattering measurement studied the structural phase transition in a GeTe conventional cell28. This makes the comparison between primitive cell calculations and conventional cell measurements more challenging. One solution is to perform calculations based on conventional unit cells, however, the number of vibrational degrees of freedom of the conventional unit cell can be 2–3 times larger than those of the primitive unit cell, thus giving a much larger computational burden29. For example, the hexagonal conventional cell of Bi2Te3 has 15 atoms with 45 vibrational modes, which is three times larger compared to that of the primitive unit cell. As a result, calculations of material properties based on conventional geometries, especially quasiparticle dynamics and transport phenomena (to be discussed in this work), are more expensive because of the increased cost in calculating interactions. In particular, lowest order three-phonon scattering processes scale to the cubic power of the number of degrees of freedom, thus phonon transport calculations for the conventional cell of Bi2Te3 will cost around 27 times more than those of the primitive unit cell.

It is well known that the choice of unit cell is arbitrary and should underly the same physical behaviors; however, it is not well known how to use symmetry relations to make lattice dynamical calculations more efficient in larger unit cells. Taking advantage of both the convenience of conventional geometries and lower degrees of freedom of primitive geometries, we previously proposed and demonstrated an efficient dynamic method for calculating quasiparticle dispersions, phonons in particular30. There, we discussed that mapping of conventional cell dynamics to primitive cell dynamics using internal translational symmetries saves computational cost in dynamical matrix diagonalizations and leads to phase interference conditions demonstrated through comparison of calculated phonon dispersions with measured inelastic neutron scattering spectra.

In this work, we apply this primitive to conventional cell dynamic method based on primitive translational symmetry (PTS) to thermal transport calculations limited by anharmonic interactions. In conventional geometries, this PTS dynamic method significantly reduces the computational cost of thermal conductivity calculations by reducing the quasiparticle scattering phase space through a conservation rule that is hidden in typical conventional dynamics and reducing the number of summations in scattering matrix element calculations. We demonstrate the convenience of this PTS method by calculating phonon transport properties based on DFT for three materials from different space groups: GeTe with space group R3m, solid N2 with space group I213, and ferromagnetic CrCl3 with space group R-3. These are representative materials with different levels of complexity in conventional cells and have important applications in thermoelectrics (GeTe)31,32,33,34, energy storage (N2)35,36, and two-dimensional magnetism (CrCl3)37,38,39,40,41.

Results

PTS dynamics in conventional basis

Figure 1 depicts the primitive (a) and conventional (b) geometries for GeTe with space group R3m. In this figure, the translation vector \({\bf{S}}\) relates the primitive and conventional geometries and is fractional within the conventional cell. In general, \({\bf{S}}\) is (nxR1 + nyR2 + nzR3)/N where R1, R2, and R3 are conventional lattice basis vectors, N is the number of layers in the conventional unit cell, and nx, ny, and nz are integers related to a particular material’s space group (discussed below). This fractional translational symmetry can be used to reorganize the dynamical matrix of the conventional unit cell and more explicitly demonstrate the vibrational phases between layers:30

$${D}_{\alpha \beta }^{k{k}^{{\prime}}}({\bf{q}},l)=\frac{1}{\sqrt{{m}_{k}{m}_{{k}^{{\prime}}}}}{\sum}_{{h}^{{\prime}}{p}^{{\prime}}}{\varPhi }_{\alpha \beta }^{00k,{h}^{{\prime}}{p}^{{\prime}}{k}^{{\prime}}}{e}^{i{\bf{q}}{\cdot}{{\bf{R}}}_{p{\prime} }}{e}^{i{h}^{{\prime} }({\bf{q}}{\cdot}{\bf{S}}{{+}}2\pi l/N{{)}}}$$
(1)

where Greek subscripts are Cartesian directions, k loops over the atoms in a single layer of the conventional unit cell, q is a phonon wavevector, mk is the mass of atom k, \({{\bf{R}}}_{p^{{\prime} }}\) is a lattice vector locating the \(p^{{\prime} }\) conventional cell (integer multiples of the basis vectors above), N is the number of ‘primitive’ layers within the conventional cell, \({h}^{{\prime} }\) represents each of these layers ranging from 0 to N-1 (see labels in Fig. 1 for GeTe), \({\varPhi }_{\alpha \beta }^{00k,{h}^{{\prime} }{p}^{{\prime} }{k}^{{\prime} }}\) are harmonic interatomic force constants (IFCs) between atom \(k\) in the origin layer and atom \(k^{\prime}\) in the \(p^{\prime}\) conventional cell in layer \(h^{\prime}\), and \({e}^{{il}{h}^{{\prime} }2\pi /N}\) is the vibrational phase between layers. Derivation of Eq. 1 is provided in Supplementary Note 1. In this formalism, l is an integer representing the phase degree of freedom between the layers that also ranges from 0 to N-1. It acts as a conserved quantity in scattering processes (via phase interference conditions) within the conventional unit cell (discussed below) that limits the scattering phase space30. Equation 1 represents N \(3n\times 3n\) dynamical matrices for each q, which gives the same number of phonon branches as those from the single conventional \(3{Nn}\times 3{Nn}\) dynamical matrix, where \(n\) is the number of atoms in one layer of the conventional cell. Diagonalization of Eq. 1 for N \(3n\times 3n\) matrices is computationally more efficient than diagonalizing a single \(3{Nn}\times 3{Nn}\) matrix. Note that the primitive cell dynamics consists of a single \(3n\times 3n\) dynamical matrix.

Fig. 1: Crystal structure of GeTe (R3m space group).
figure 1

a Two-atom primitive unit cell where each pair of lattice vectors forms an angle of 57.8°. b Six-atom conventional unit cell with three layers related by translation vector S (green arrows). The black rectangle highlights that a single layer can be used to describe the dynamics in a conventional unit cell.

Full size image

Specifically, for GeTe (space group: R3m) with conventional unit cell shown in Fig. 1b, there are three layers with two atoms in each (\(N=3\) and \(n=2\)). The translation vector \({\bf{S}}=(2{{\bf{R}}}_{1}+{{\bf{R}}}_{2}+{{\bf{R}}}_{3})/3=[a/2,a/(2\sqrt{3}),c/3]\) where \(a\) and \(c\) are in-plane and out-of-plane lattice parameters, respectively. Instead of one conventional \(18 \times 18\) dynamical matrix, Eq. (1) with PTS dynamics gives three \(6\times 6\) dynamical matrices for \(l=\left(\mathrm{0,1,2}\right)\) for each q.

The PTS dynamics of Eq. 1 can be applied to 81 different space groups with symmetry operations x + nx/N, y + ny/N, z + nz/N, where N is either 2 or 3, and nα ranges from 0 to 2 with \({n}_{x}{n}_{y}{n}_{z}\,\ne\, 0\). Table 1 summarizes the allowed space groups corresponding to such operations. We demonstrate the application of PTS dynamics on representative materials: GeTe with space group R3m, solid N2 with space group I213, and ferromagnetic CrCl3 with space group R-3. We note that the solid N2 system considered here has not been synthesized, presumably due to predicted thermodynamic instability42, but provides a good test case of our method due to N2 having the I213 space group and a relatively simple structure. Complete results for GeTe and temperature dependent thermal conductivities for all materials are presented in the main text. Structures, phonon dispersions, and scattering rates for CrCl3 and solid N2 are described in a previous study30 and in Supplementary Figs. 2–4, respectively. Details of DFT calculations for all materials are provided in Methods.

Table 1 Space groups for which PTS dynamics can be applied to conventional unit cells.
Full size table

Phonon dispersion

We first examine the vibrational properties of GeTe based on its conventional unit cell (Fig. 1b). Diagonalizing Eq. (1) for each value of \(l\) = (0, 1, 2) gives the phonon dispersion depicted in Fig. 2 with six branches for each vibrational phase depicted by different colors. The phonon dispersions from the PTS dynamics overlap those from the conventional method (underlying black solid curves) as observed along high-symmetry lines and an arbitrary direction in reciprocal space. The colored dispersions provide insights into observables from scattering experiments in conventional geometries as only phonon branches with \(l=0\) give non-zero spectral intensities, which has been demonstrated for CrCl330, and other materials with internal twist symmetries43,44.

Fig. 2: Primitive and conventional phonon dispersions of GeTe.
figure 2

Underlying dashed black curves give the phonon dispersions of GeTe from conventional dynamics, while the colored solid curves give these from PTS dynamics, both along high-symmetry lines and an arbitrary direction from the Γ point to (\({q}_{x}=0.123\pi /a\), \({q}_{y}=0.456\pi /b\), \({q}_{z}=0.789\pi /c\)). Blue, green, and red curves correspond to integers l = 0, l = 1, and l = 2, respectively.

Full size image

The wavevector q and integer \(l\) are not limited to the FBZ and can be calculated for extended Brillouin zones. When comparing phonons from different zones, for instance mapping extended zones back to the FBZ, which is necessary for studying Umklapp phonon scattering, \(\Delta {\bf{q}}\) relating phonons between zones is given by a reciprocal lattice vector: \(\Delta {\bf{q}}={\bf{G}}=\left[{g}_{1}{{\bf{G}}}_{1},{g}_{2}{{\bf{G}}}_{2},{g}_{3}{{\bf{G}}}_{3}\right]\), an integer (\({g}_{i}\)) multiple of the conventional reciprocal lattice basis vectors (\({{\bf{G}}}_{i}\))29. In this mapping, q and \(l\) are related through a conserved phase relation obtained from Eq. (1): \({\mathbf{\Delta}}{\bf{q}}\,\cdot\,{\bf{S}}+2\pi \Delta l/N=0\) (see Methods for a detailed derivation), where the change of integer l is given by \(\Delta l={l}_{G}=-\left({n}_{x}{g}_{1}+{n}_{y}{g}_{2}+{n}_{z}{g}_{3}\right)\). Again, nα are given by the space group of the crystal (see Table 1). This relationship was demonstrated for bulk CrCl3 when comparing calculated dispersions to scattering measurements in varying Brillouin zones30 and is shown for GeTe in Supplementary Fig. 1 (see Methods for details).

Scattering conservation rule from phase relations

The PTS dynamic method also elucidates conservation conditions for intrinsic phonon-phonon interactions enforced by phase relations represented by the integer l that are hidden in the conventional dynamics used to build phonon linewidths30. This scattering is critically important when trying to determine and understand thermal transport in conventional geometries.

To demonstrate this, we start with anharmonic three-phonon transition rates from quantum perturbation theory (i.e., Fermi’s golden rule):45,46,47

$${\varGamma }_{\lambda {\lambda }^{{\prime} }{\lambda }^{\prime\prime}}^{\pm }=\frac{{{\hbar }}\pi }{4{N}_{{\rm{uc}}}}\frac{\left({f}_{\lambda }^{0}+1\right)\left({f}_{{\lambda }^{{\prime} }}^{0}+1/2\pm 1/2\right){f}_{{\lambda }^{\prime\prime}}^{0}}{{\omega }_{\lambda }{\omega }_{{\lambda }^{{\prime} }}{\omega }_{{\lambda }^{\prime\prime}}}{\left|{V}_{\lambda {,\pm \lambda }^{{\prime} }{,-\lambda }^{\prime\prime}}^{(3)}\right|}^{2}\delta \left({\omega }_{\lambda }\pm {\omega }_{{\lambda }^{{\prime} }}-{\omega }_{{\lambda }^{\prime\prime}}\right)$$
(2)

where \(\lambda\) represents a phonon mode with wavevector \({\bf{q}}\) and polarization \(J\), the Dirac delta function \(\delta\) ensures energy conservation, \({N}_{{\rm{uc}}}\) is the number of conventional unit cells, primes label different phonon modes, \({V}_{\lambda {,\lambda }^{{\prime} }{,\lambda }^{{\text{'}\text{'}}}}^{(3)}\) are three-phonon scattering matrix elements, and \({f}_{\lambda }^{0}\) is the Bose-Einstein distribution for mode \(\lambda\). The \(\pm\) designate two different types of three-phonon interactions, coalescence and decay48. The scattering matrix elements in conventional dynamics are

$${V}_{\lambda {,\lambda }^{{\prime} },{\lambda }^{\prime\prime}}^{(3)}={\sum }_{K}{\sum }_{{p}^{{\prime} }{K}^{{\prime} }}{\sum }_{{p}^{\prime\prime}{K}^{\prime\prime}}{\sum }_{\alpha \beta \gamma }{\varPhi }_{\alpha \beta \gamma }^{0K,{p}^{{\prime} }{K}^{{\prime} },{p}^{\prime\prime}{K}^{\prime\prime}}{e}^{i{{\bf{q}}}^{{\prime} }{{\cdot }}{{\bf{R}}}_{p^{\prime} }}{e}^{i{{\bf{q}}}^{\prime\prime}{{\cdot }}{{\bf{R}}}_{p^{\prime\prime}}}\frac{{\xi }_{\alpha K}^{\lambda }{\xi }_{\beta {K}^{{\prime} }}^{{\lambda }^{{\prime} }}{\xi }_{\gamma {K}^{\prime\prime}}^{{\lambda }^{\prime\prime}}}{\sqrt{{m}_{K}{m}_{{K}^{{\prime} }}{m}_{{K}^{\prime\prime}}}}$$
(3)

where \(K\) refers to atoms in the conventional cell, \({\xi }_{\alpha K}^{\lambda }\) is the \(\alpha\) component of an eigenvector in the conventional basis for phonon mode \(\lambda\), and \({\varPhi }_{\alpha \beta \gamma }^{0K,{p}^{{\prime} }{K}^{{\prime} },{p}^{{\prime\prime}}{K}^{{\prime\prime}}}\) are third-order anharmonic IFCs47. Translational invariance has already been enforced in Eq. (3), which leads to crystal momentum conservation for three-phonon interactions:49

$${\bf{q}}\pm {{\bf{q}}}^{{\prime} }-\left({{\bf{q}}}^{{\prime\prime}}+{\bf{G}}\right)=0$$
(4)

that limits the scattering phase space.

We now discuss an additional conservation rule of the integer l in the PTS dynamic method. An eigenvector in the conventional basis is related to that in PTS dynamics through30

$${\xi }_{\alpha K}^{\lambda }=\frac{1}{\sqrt{N}}{\varepsilon }_{\alpha k}^{\tilde{\lambda }}{e}^{{ih}\left({\bf{q}}{\cdot}{\bf{S}}+2\pi l/N\right)}$$
(5)

Where \({\varepsilon }_{\alpha k}^{\tilde{\lambda }}\) is the \(\alpha\) component of the PTS eigenvector for atom \(k\) in a single layer. \(\widetilde{\lambda }\) represents a phonon mode with wavevector q, phase integer \(l\), and polarization j (as determined by diagonalization of Eq. 1). Combining Eqs. (3) and (5) and replacing atom notation \(\left(K\right)\) with \(\left(h,k\right)\) and phonon state notation \(\lambda\) (\({\bf{q}},J\)) with \(\widetilde{\lambda }\) (\({\bf{q}},j,l\)), as with Eq. 1, we have

$$\begin{array}{ll}{V}_{{\tilde{\lambda }},{\pm {{\tilde{\lambda }}}}^{{\prime} },{-{{\tilde{\lambda }}}}^{\prime\prime}}^{(3)}={\left(\frac{1}{\sqrt{N}}\right)}^{3}{\sum }_{h0k}{\sum }_{{{h}^{{\prime} }p}^{{\prime} }{k}^{{\prime} }}{\sum }_{{h}^{\prime\prime}{p}^{\prime\prime}{k}^{\prime\prime}}\\ \qquad\qquad\quad\times\,{\sum }_{\alpha \beta \gamma }{\varPhi }_{\alpha \beta \gamma }^{h0k,{{h}^{{\prime} }p}^{{\prime} }{k}^{{\prime} },{{h}^{\prime\prime}p}^{\prime\prime}{k}^{\prime\prime}}{e}^{\pm i{{\bf{q}}}^{{\prime} }{\cdot}{{\bf{R}}}_{p{\prime} }}{e}^{-i{{\bf{q}}}^{\prime\prime}{\cdot}{{\bf{R}}}_{p^{\prime\prime}}}\\ \qquad\qquad\quad\times \,\frac{{\varepsilon }_{\alpha k}^{\tilde{\lambda }}{\varepsilon }_{\beta {k}^{{\prime} }}^{{\pm {{\tilde{\lambda }}}}^{{\prime} }}{\varepsilon }_{\gamma {k}^{\prime\prime}}^{{-\tilde{\lambda }}^{\prime\prime}}}{\sqrt{{m}_{k}{m}_{{k}^{{\prime} }}{m}_{{k}^{\prime\prime}}}}{e}^{{ih}\left({\bf{q}}{\cdot}{\bf{S}}+2\pi l/N\right)}{e}^{\pm i{h}^{\prime}\left({{\bf{q}}}^{\prime}{\cdot}{\bf{S}}+2\pi {l}^{{\prime} }/N\right)}{e}^{-i{h}^{\prime\prime}\left({{\bf{q}}}^{\prime\prime}{\cdot}{\bf{S}}+2\pi {l}^{\prime\prime}/N\right)}\end{array}$$
(6)

The interatomic potential (to all perturbative orders) is invariant with respect to translation by a lattice vector (in both primitive and conventional bases), thus anharmonic IFCs are invariant with respect to integer multiples of the layer translation vector S so that

$${\varPhi }_{\alpha \beta \gamma }^{h0k,{{h}^{{\prime} }p}^{{\prime} }{k}^{{\prime} },{{h}^{\prime\prime}p}^{\prime\prime}}{k}^{\prime\prime}={\varPhi }_{\alpha \beta \gamma }^{00k,{\left({h}^{{\prime} }-\,h\right)p}^{{\prime} }{k}^{{\prime} },{\left({h}^{\prime\prime}-\,h\right)p}^{\prime\prime}{k}^{\prime\prime}}$$
(7)

Note that the unit cell \({p}^{{\prime} }\) (\({p}^{{\prime}{\prime}}\)) will be shifted if \({h}^{{\prime} }\, <\, h\) (\({h}^{\prime\prime}\, < \,h\)). The layer phase relations in Eq. (6) are also unchanged if shifted by h:

$${e}^{{ih}\left({\bf{q}}{\cdot}{\bf{S}}+2\pi l/N\right)}{e}^{\pm i{h}^{{{{\prime}}}}\left({{\bf{q}}}^{{{{\prime}}}}{\cdot}{\bf{S}}{{+}}2\pi {l}^{{\prime}}/N\right)}{e}^{-i{h}^{{\prime}{\prime}}\left({{\bf{q}}}^{{{{\prime}{\prime}}}}{\cdot}{\bf{S}}{{+}}2\pi {l}^{{{\prime}{\prime}}}/N\right)}{{=}}{e}^{\pm i\left({h}^{{\prime}}-h\right)\left({{\bf{q}}}^{{{{\prime} }}}{\cdot}{\bf{S}}{{+}}2\pi {l}^{\prime}/N\right)}{e}^{-i\left({h}^{{{{\prime}{\prime}}}}-h\right)\left({{\bf{q}}}^{{{{\prime}{\prime}}}}{\cdot}{\bf{S}}{{+}}2\pi {l}^{{{{\prime}{\prime}}}}/N\right)}$$
(8)

Rearranging Eq. 8 and using Eq. 4 lead to the following constraint

$${e}^{{ih}{\bf{S}}{\cdot}{\bf{G}}}{e}^{i2\pi h\left(l\pm {l}^{{\prime} }-{l}^{{{{\prime}{\prime}}}}\right)/N}=1$$
(9)

Thus,

$$h\left[{-l}_{G}+\left(l\pm {l}^{{\prime} }-{l}^{{{{\prime}{\prime}}}}\right)\right]=0$$
(10)

where \({l}_{G}=-N{\bf{S}}{\cdot}{\bf{G}}/\left(2\pi \right)\) is an integer (see Methods for further discussion of the relationship of \({\bf{G}}\) and \({l}_{G}\)). Since \(h\) can be any integer, the expression in the bracket must be zero

$$l\pm {l}^{{\prime} }-\left({l}^{\prime\prime}+{l}_{G}\right)=0$$
(11)

Equation (11) reveals another conservation rule in terms of integer \(l\) that is hidden in conventional methods. Explicitly using this conservation condition reduces the scattering phase space by ~1/N and enables significantly more efficient transport calculations (described below). We note that the primitive and conventional cell calculations have the additional phase symmetry constraint (Eq. 11) built in, however, the primitive cell has awkward geometries and the conventional cell does not explicitly exploit it, simply giving no strength to scattering interactions that violate Eq. 11.

Scattering and thermal transport

We now continue to show how scattering matrix elements \({V}_{{\tilde{\lambda }},{\pm {\tilde{\lambda }}}^{{\prime}},{-\tilde{\lambda}}^{\prime\prime}}^{(3)}\) are calculated to obtain scattering rates (Eq. (2)) in PTS dynamics. Applying translational invariance to the third-order potential expansion (and relabeling the arbitrary layer indices in Eqs. 7 and 8), the scattering matrix elements in PTS dynamics become

$$\begin{array}{l}{V}_{{\tilde{\lambda }},{\pm {\tilde{\lambda }}}^{{\prime} },{-{\tilde{\lambda }}}^{\prime\prime}}^{(3)}={\left(\frac{1}{\sqrt{N}}\right)}^{3}{\sum }_{h0k}{\sum }_{{{h}^{{\prime} }p}^{{\prime} }{k}^{{\prime} }}{\sum }_{{h}^{\prime\prime}{p}^{\prime\prime}{k}^{\prime\prime}}\\ \qquad\qquad\quad\times\,{\sum }_{\alpha \beta \gamma }{\varPhi }_{\alpha \beta \gamma }^{00k,{{h}^{{\prime} }p}^{{\prime} }{k}^{{\prime} },{{h}^{\prime\prime}p}^{\prime\prime}{k}^{\prime\prime}}{e}^{i{{\bf{q}}}^{{\prime} }{\cdot}{{\bf{R}}}_{{{{\rm{p}}^{\prime}}}}}{e}^{i{{\bf{q}}}^{\prime\prime}{\cdot}{{\bf{R}}}_{{p}^{\prime\prime}}}\\ \qquad\qquad\quad\times \,\frac{{\varepsilon }_{\alpha k}^{\tilde{\lambda }}{\varepsilon }_{\beta {k}^{{\prime} }}^{\pm {\tilde{\lambda }}^{{\prime} }}{\varepsilon }_{\gamma {k}^{\prime\prime}}^{-{\tilde{\lambda }}^{\prime\prime}}}{\sqrt{{m}_{k}{m}_{{k}^{{\prime} }}{m}_{{k}^{\prime\prime}}}}{e}^{\pm i{h}^{\prime}\left({{\bf{q}}}^{\prime}{\cdot}{\bf{S}}+2\pi {l}^{{\prime} }/N\right)}{e}^{-i{h}^{\prime\prime}\left({{\bf{q}}}^{\prime\prime}{\cdot}{\bf{S}}+2\pi {l}^{\prime\prime}/N\right)}\end{array}$$
(12)

Since the terms in the summation are independent of h, the summation over h gives a factor N and Eq. (12) becomes

$$\begin{array}{l}{V}_{{\tilde{\lambda }},{\pm {\tilde{\lambda }}}^{{\prime} },{-{\tilde{\lambda }}}^{\prime\prime}}^{(3)}=\frac{1}{\sqrt{N}}{\sum }_{00k}{\sum }_{{{h}^{{\prime} }p}^{{\prime} }{k}^{{\prime} }}{\sum }_{{h}^{\prime\prime}{p}^{\prime\prime}{k}^{\prime\prime}}\\ \qquad\qquad\quad\times\,{\sum }_{\alpha \beta \gamma }{\varPhi }_{\alpha \beta \gamma }^{00k,{{h}^{{\prime} }p}^{{\prime} }{k}^{{\prime} },{{h}^{\prime\prime}p}^{\prime\prime}{k}^{\prime\prime}}{e}^{i{{\bf{q}}}^{{\prime} }{\cdot}{{\bf{R}}}_{{{{\rm{p}}^{\prime}}}}}{e}^{i{{\bf{q}}}^{\prime\prime}{\cdot}{{\bf{R}}}_{{{{\rm{p}}^{\prime\prime}}}}}\\ \qquad\qquad\quad\times\, \frac{{\varepsilon }_{\alpha k}^{\tilde{\lambda }}{\varepsilon }_{\beta {k}^{{\prime} }}^{\pm {\tilde{\lambda }}^{{\prime} }}{\varepsilon }_{\gamma {k}^{\prime\prime}}^{{-{\tilde{\lambda }}}^{\prime\prime}}}{\sqrt{{m}_{k}{m}_{{k}^{{\prime} }}{m}_{{k}^{\prime\prime}}}}{e}^{\pm i{h}^{{{{\prime} }}}\left({{\bf{q}}}^{\prime}{\cdot}{\bf{S}}+2\pi {l}^{{\prime} }/N\right)}{e}^{-i{h}^{\prime\prime}\left({{\bf{q}}}^{\prime\prime}{\cdot}{\bf{S}}+2\pi {l}^{\prime\prime}/N\right)}\end{array}$$
(13)

We show in Fig. 3a that the scattering rates of GeTe from PTS dynamics are identical (within numerical accuracy) to those from conventional dynamics. Figure 3b shows the transition rates49 for a representative phonon mode at the Γ point with (ω = 2.61 THz, j = 4) as a function of the frequency of one of the interacting phonons (\({\omega }_{{\lambda }^{{\prime} }}\)) calculated from conventional dynamics. Here, PTS dynamics was used to identify the integer \(l\) for each of the three phonon modes involved in the interactions to determine whether conservation of l occurred or not: \(\Delta l=l\pm {l}^{{\prime} }-({l}^{\prime\prime}+{l}_{G})=0\) or \(\Delta l\,\ne\, 0\). In Fig. 3b, the transition rates for \(\Delta l\,\ne\, 0\) are numerically zero compared to those for \(\Delta l=0\). This demonstrates that conservation of l, which is mandated by the internal translational symmetry, is hidden in the conventional formalism suggesting that this can be used to reduce computational cost in thermal transport calculations.

Fig. 3: Total scattering and individual transition rates for phonons in GeTe.
figure 3

a Scattering rates as a function of phonon frequency from separate conventional (blue) and PTS (red) dynamics. Inset: Direct comparison of scattering rates from the two methods. An underlying solid red line provides a guide for the eye showing equality. b Individual anharmonic phonon transition rates for a phonon mode at the Γ point with (ω = 2.61 THz, j = 4) as a function of the frequency of one of the other interacting modes in three-phonon interactions. Transition rates with \(\Delta l\,\ne \,0\), i.e., violating conservation of l, are numerically zero.

Full size image

The thermal conductivity matrix from the PTS dynamic method is calculated as

$${\mathbf{\kappa }}=\mathop{\sum}\nolimits_{\tilde{\lambda }}{C}_{\tilde{\lambda }}{{\bf{v}}}_{\tilde{\lambda }}^{2}{\tau }_{\tilde{\lambda }}$$
(14)

where \({C}_{\tilde{\lambda}}\) is the volumetric specific heat for phonon mode \(\widetilde{\lambda }\) with wavevector q, polarization j, and phase integer l, \({\bf{v}}\) is the phonon velocity vector, and \({\tau}_{\tilde{\lambda }}\) is the phonon lifetime due to intrinsic three-phonon scattering (related to the inverse of the sum of scattering rates that conserve energy, momentum, and l) and isotopic scattering (in natural samples). Note that phonon branch J in the conventional basis is broken down into j and l in PTS dynamics. This is the only difference from the widely used formula in the conventional method7,47,49,50,51. The PTS dynamic method gives identical thermal conductivities as the conventional method for both in-plane and cross-plane directions. Importantly, unlike phonon chirality previously described by a twist dynamical method43,44, here PTS is not limited to specific high-symmetry directions, and thus calculations that involve sums over the full Brillouin zone (i.e., including low-symmetry points in reciprocal space), are possible. This is an important feature for application of the PTS dynamics to transport phenomena of quasiparticles.

Figure 4 shows the thermal conductivities of bulk N2, GeTe, and CrCl3 as a function of temperature from 30 to 400 K. The crystalline N2 system considered here (a non-molecular solid that exists under high pressure35,36) has isotropic thermal transport, i.e., \({\kappa }_{{xx}}={\kappa }_{{yy}}={\kappa }_{{zz}}\) (off-diagonal terms are zero for all systems considered here), and has the highest \(\kappa\) of the three systems considered, with room temperature \(\kappa =150.35\) Wm-1K-1. We note that N2 has a reasonably high thermal conductivity, particularly at low temperature, due to its being comprised of covalently bonded light atoms in a relatively simple structure52, though it is significantly less thermally conductive than diamond53. We also note that naturally occurring N is nearly isotopically pure, thus phonon-isotope scattering does not provide significant thermal resistance in the temperature range considered here.

Fig. 4: Thermal conductivity of N2 (black curve, inset), GeTe (blue curve), and CrCl3 (red curve) from PTS dynamics as a function of temperature in bulk naturally occurring samples.
figure 4

Measurements for GeTe are represented by solid and open red circles54, squares55, and triangle56. In-plane thermal conductivity measurements for CrCl3 are represented by blue circles57.

Full size image

GeTe and CrCl3 have anisotropic thermal conductivity tensors with separate in-plane (\({\kappa }_{{xx}}={\kappa }_{{yy}}\)) and cross-plane (\({\kappa }_{{zz}}\)) values; furthermore, phonon-isotope scattering is relatively important in each (see Supplementary Fig. 5). In Fig. 4 we compare in-plane \(\kappa\) for CrCl3 and an effective thermal conductivity (\({\kappa }_{{\rm{eff}}}=\left(2{\kappa }_{{xx}}+{\kappa }_{{zz}}\right)/3\)) for GeTe with measured values54,55,56,57. Natural isotopes are included in the calculations for both cases. For GeTe, the crystal orientation of the thin film samples was not given (thus \({\kappa }_{{\rm{eff}}}\) was used here for comparison) and significant electronic contributions to \(\kappa\) are anticipated from measured resistance data. Thus, we compare our phonon calculations with both the total measured \(\kappa\) (filled red circles) for GeTe thin films and the expected lattice contribution (hollow red circles, subtracting estimated electronic \(\kappa\) via the Wiedemann-Franz law)54. For ferromagnetic CrCl3, spin-phonon interactions are expected to be important (not considered in this work), particularly at low temperatures. Here, we compare with measurements under the strongest applied in-plane magnetic field, which suppresses scattering by magnetic excitations so that nearly only intrinsic phonon scattering dominates the transport, particularly in the temperature range considered here57. Agreement with measurements is reasonable given that phonon-defect interactions (i.e., from point defects, grain boundaries, surfaces, and other extended point defects) have not been considered. This suggests that the samples are relatively crystalline and pure.

Computational efficiency

To demonstrate the improved efficiency of PTS dynamics over the conventional dynamics, we calculate the ratio of total cpu hours required for PTS and conventional thermal conductivity calculations for GeTe, N2, and CrCl3 (having different degrees of complexity) for different integration grid densities. The computational savings is generally independent of the density of the sampling q mesh (Supplementary Fig. 6), and the average (over different grid densities) computational cost ratios are 0.3942\(\pm 0.0151\) for N2 (N = 2, n = 4), 0.1334\(\pm 0.0051\) for GeTe (N = 3, n = 2), and 0.0925\(\pm 0.0054\) for CrCl3 (N = 3, n = 8). These results demonstrate that using PTS in the calculation of thermal transport in the conventional basis can give ~N2 numerical savings, particularly when the complexity of the material increases (as in going from N2 to GeTe to CrCl3). We also tested a very large toy system developed from GeTe and observed a similar speedup (see Supplementary Note 2 and Supplementary Fig. 7 for details). Overall, PTS dynamics in the conventional basis is significantly more efficient than conventional dynamics, which is especially beneficial for calculations of large complex conventional cells.

The efficiency of the PTS dynamic method for calculating phonon transport in a conventional basis comes in two-fold: reduced dynamical matrix sizes for harmonic calculations and reduced phase space and matrix sums for anharmonic calculations. In our particular transport framework3, the reciprocal space is discretized and phonon harmonic properties are calculated once for a given q mesh. Thus, the numerical savings for dynamical matrix diagonalizations are negligible compared with the much heavier anharmonic calculations. We find that the reduction of ~1/N2 in computational cost for GeTe and CrCl3 is primarily derived from two factors in anharmonic calculations. First, the conservation of integer \(l\) in Eq. (11) leads to reduced phase space that is ~1/N compared to that in conventional dynamics. Moreover, the first summation in Eq. (13) is only for atoms in a single layer rather than all atoms in the conventional cell, which contributes to another 1/N factor in the calculation of scattering rates, which is the primary numerical cost in thermal conductivity calculations.

Discussion

It is worth mentioning that our calculation framework, i.e., precalculated harmonic phonon properties on a fixed q mesh, is the same as that in widely used packages like ShengBTE3, almaBTE4, ALAMODE5, and phono3py6. Thus, the proposed PTS dynamic method can achieve a similar computational savings of ~1/N2 with appropriate modifications in these applications.

Also, this PTS dynamics can be particularly efficient when applied to higher-order anharmonic scattering, for example, four-phonon scattering58,59. In such cases, the scattering rates and scattering matrix elements are \({\varGamma }_{\lambda {\lambda }^{\prime}{\lambda }^{{\prime\prime}}{\lambda }^{{\prime\prime\prime}}}^{\pm }\) and \({V}_{\lambda {,\pm \lambda }^{{\prime} }{,\pm \lambda }^{{\prime\prime}}{,-\lambda }^{{\prime\prime\prime}}}^{(4)}\), respectively, where \({\lambda }^{{\prime\prime\prime}}\) represents a fourth phonon mode58. Replacing the conventional description with PTS dynamics and following the same derivation, the conservation of integer \(l\) in Eq. (11) becomes \(l\pm {l}^{{\prime} }\pm {l}^{\prime\prime}-({l}^{\prime\prime\prime}+{l}_{G})=0\) where \({l}^{\prime\prime\prime}\) corresponds to \({\lambda }^{{\prime\prime\prime}}\). A similar reduction of ~1/N2 in computational cost is expected, which is especially beneficial considering that the four-phonon scattering rate calculations are thousands of times more expensive than those of three-phonon calculations58.

In summary, we discussed lattice dynamics and phonon transport using an efficient dynamic method that uses primitive translational symmetry (PTS) within conventional cells. Such symmetry relations do not introduce new physical behaviors but can significantly simplify calculations in large complex unit cells. The underlying phase dynamics of quasiparticle interactions elucidate quantum phase interference conditions hidden in anharmonic phonon interactions within conventional unit cells, leading to conservation rules that can be exploited to make more efficient calculations of quasiparticle dynamics in larger, more convenient conventional unit cells, comparable to those of more awkward primitive unit cells. The proposed dynamics accurately describes transport phenomena and costs significantly less computationally compared to conventional dynamics, which is valuable for studying quasiparticle interactions in complex material systems.

Methods

Density functional theory

We performed density functional theory calculations using the projector augmented wave method (PAW)60, as implemented in the Vienna Ab-initio Simulation Package (VASP)61,62,63, for all the materials considered. The generalized gradient approximation, parameterized by Perdew, Burke, and Ernzerhof64 was used for exchange-correlations.

The plane wave energy cut-off is 500 eV (520 eV for CrCl3), and energy convergence criteria is 10-6 eV. Ionic relaxations were performed until Hellmann-Feynman forces converged to 10-5 eV/Å for GeTe and solid N2 (10-4 eV/Å for CrCl3). The structures were optimized with Γ-centered 9 × 9 × 3 (GeTe), 9 × 9 × 9 (solid N2), and 4 × 4 × 4 (CrCl3) k-meshes. Harmonic IFCs were calculated using the phonopy package65 with a 4 × 4 × 2 supercell and Γ-centered 3 × 3 × 1 k-mesh for GeTe, a 3 × 3 × 3 supercell and Γ-centered 3 × 3 × 3 k-mesh for solid N2, and a 2 × 2 × 1 supercell and Γ-centered 2 × 2 × 2 k-mesh for CrCl3.

The anharmonic IFCs were calculated by finite displacement method using the thirdorder.py package3. For GeTe, a 300-atom 5 × 5 × 2 supercell, Γ-only k-mesh, and a cutoff distance up to the 5th nearest neighbors (NN) were used. For solid N2, a 216-atom 3 × 3 × 3 supercell, Γ-centered 3 × 3 × 3 k-mesh, and a cutoff distance up to the 3rd NN were used. For CrCl3, the interlayer spacing is much larger than the in-plane nearest neighbor atomic distances, thus we used a cylindrical cutoff that is 0.43 nm in the plane and 0.6 nm across the plane. The cutoffs correspond to the 7th NN in the plane and include most of the atoms in adjacent layers across the plane. A 192-atom 2 × 2 × 2 supercell and Γ-only k-mesh were used. All the other calculation parameters are the same as those used for harmonic calculations.

Zone folding relation between q and l

Translational invariance of the harmonic interatomic potential with respect to PTS requires that

$${{e}^{i{\bf{q}}{\cdot}\left({{\bf{R}}}_{{p}^{{\prime} }}+{h}^{{\prime} }{\bf{S}}\right)}}{{e}^{\frac{{il}{h}^{{\prime} }2\pi }{N}}}={{e}^{i{\bf{q}}{\cdot}\left({{\bf{R}}}_{{p}^{{\prime} }}+\left({h}^{{\prime} }{\boldsymbol{-}}h\right){\bf{S}}\right)}}{e}^{{il}({h}^{{\prime} }-h)2\pi /N}={{e}^{i{\bf{q}}{\cdot}\left({{\bf{R}}}_{{p}^{{\prime} }}+{h}^{{\prime} }{\bf{S}}\right)}}{{e}^{i{\bf{q}}{\cdot}\left({\boldsymbol{-}}h{\bf{S}}\right)}}{{e}^{\frac{{il}{h}^{{\prime} }2\pi }{N}}}{{e}^{-\frac{{ilh}2\pi }{N}}}$$
(15)

Rearranging Eq. (15) and comparing terms leads to

$${e}^{i{\bf{q}}{\cdot}{\bf{S}}}{e}^{\frac{{il}2\pi }{N}}=1$$
(16)

and thus

$${\bf{q}}{\cdot}{\bf{S}}+\frac{2\pi }{N}l=0$$
(17)

where the translation vector \({\bf{S}}\) is \(\left[\frac{{n}_{x}}{N}{{\bf{R}}}_{1},\frac{{n}_{y}}{N}{{\bf{R}}}_{2},\frac{{n}_{z}}{N}{{\bf{R}}}_{3}\right]\). For extended zones, where \({\bf{q}}\) is changed by \(\Delta {\bf{q}}\,{\boldsymbol{=}}\,{\bf{G}}=[{g}_{1}\frac{2\pi }{{{\bf{R}}}_{1}},{g}_{2}\frac{2\pi }{{{\bf{R}}}_{2}},{g}_{3}\frac{2\pi }{{{\bf{R}}}_{3}}]\), the change of integer l is thus

$$\Delta l={l}_{G}=-\frac{N}{2\pi }{\bf{G}}{\cdot}{\bf{S}}{\boldsymbol{=}}{\boldsymbol{-}}\left({n}_{x}{g}_{1}+{n}_{y}{g}_{2}+{n}_{z}{g}_{3}\right)$$
(18)

This relationship is shown via the phonon dispersion across two BZs for GeTe in Supplementary Fig. 1.

Thermal conductivity calculations

The numerical cost of the thermal conductivity calculations was computed for a variety of reciprocal space integration mesh densities for each system (Supplementary Fig. 6). For the data presented in Fig. 4 the q-mesh samplings used are 29 × 29 × 29 (N2), 25 × 25 × 10 (GeTe), and 15 × 15 × 5 (CrCl3), for which thermal conductivities are converged.

For both conventional and PTS dynamic methods, thermal conductivities are calculated by solving the steady-state Peierls-Boltzmann equation under the relaxation time approximation. The Dirac delta functions for energy conservation in Eq. (2) are approximated by adaptive Gaussian functions as implemented in the ShengBTE package3.