Primitive to conventional geometry projection for efficient phonon transport calculations

The primitive Wigner-Seitz cell and corresponding first Brillouin zone (FBZ) are typically used in calculations of lattice vibrational and transport properties as they contain the smallest number of degrees of freedom and thus have the cheapest computational cost. However, in complex materials, the FBZ can take on irregular shapes where lattice symmetries are not apparent. Thus, conventional cells (with more atoms and regular shapes) are often used to describe materials, though dynamical and transport calculations are more expensive. Here we discuss an efficient anharmonic lattice dynamic method that maps conventional cell dynamics to primitive cell dynamics based on translational symmetries. This leads to phase interference conditions that act like conserved quantum numbers and a conservation rule for phonon scattering that is hidden in conventional dynamics which significantly reduces computational cost. We demonstrate this method for phonon transport in a variety of materials with inputs from first-principles calculations and attribute its efficiency to reduced scattering phase space and fewer summations in scattering matrix element calculations .


INTRODUCTION
Computational studies of material properties, particularly transport phenomena of heat, charge, and mass at micro-and nanoscales 1 , 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 tools [2][3][4][5][6] .Firstprinciples calculations based on density functional theory (DFT) have demonstrated remarkable power in accurately predicting thermal and electrical properties of semiconductors [7][8][9][10] , for example, the ultrahigh thermal conductivity of cubic BAs and BN [11][12][13][14] , phonon hydrodynamics in graphitic materials [15][16][17] , and high mobility of BAs 18 .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 resistance 19,20 , anomalous thermal Hall effect [21][22][23] , and topological superconductivity [24][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 cells 27 (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 crossplane) 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.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 burden 28 .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.
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 particular 29 .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 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) [30][31][32][33] , energy storage (N2) 34,35 , and two-dimensional magnetism (CrCl3) [36][37][38][39][40] .

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  relates the primitive and conventional geometries.In general,  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 translational symmetry can be used to reorganize the dynamical matrix of the conventional unit cell and more explicitly demonstrate the vibrational phases between layers 29 : 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 ,  ′ is a lattice vector locating the ′ conventional cell (integer multiples of the basis vectors above), N is the number of 'primitive' layers within the conventional cell, ℎ ′ represents each of these layers ranging from 0 to N-1 (see labels in Fig. 1 for GeTe),   00,ℎ ′  ′  ′ are harmonic interatomic force constants (IFCs) between atom  in the origin layer and atom ′ in the ′ conventional cell in layer ℎ′, and  ℎ ′ 2  ⁄ is the vibrational phase between layers.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 is a conserved quantity in scattering processes (discussed below), like phonon energy and momentum 29 .Equation 1represents N 3 × 3 dynamical matrices for each q, which gives the same number of phonon branches as those from the single conventional 3 × 3 dynamical matrix, where  is the number of atoms in one layer of the conventional cell.Though Eq. 1 has N times more matrices, their smaller size makes diagonalization more efficient than that of the conventional dynamics.
Note that the primitive cell dynamics consists of a single 3 × 3 dynamical matrix.Specifically, for GeTe (space group: R3m) with conventional unit cell shown in Fig. 1b, there are three layers with two atoms in each (  = 3 and  = 2 ).The translation vector  = ⁄ ] where  and  are in-plane and out-of-plane lattice parameters, respectively.Instead of one conventional 18 × 18 dynamical matrix, Eq. ( 1) with PTS dynamics gives three 6 × 6 dynamical matrices for  = (0,1,2) for each q.

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  = (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  = 0 give non-zero spectral intensities, which has been demonstrated for CrCl3 29 , and other materials with internal twist symmetries 42,43 .The wavevector q and integer  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,  relating phonons between zones is given by a reciprocal lattice vector: an integer (  ) multiple of the conventional reciprocal lattice basis vectors (  ) 28 .In this mapping, q and  are related through a conserved phase relation obtained from Eq. ( 1) (see Methods for a detailed derivation), where the change of integer l is given by ∆ =   = −(   1 +    2 +    3 ).Again, nα are given by the space group of the crystal (see above).This relationship was demonstrated for bulk CrCl3 when comparing calculated dispersions to scattering measurements in varying Brillouin zones 29 and is shown for GeTe in Fig. S1 (see Methods for details) .

Scattering conservation rule from phase relations
The PTS dynamic method also elucidates conservation conditions for intrinsic phononphonon interactions enforced by phase relations represented by the integer l that are hidden in the conventional dynamics used to build phonon linewidths 29 .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) [44][45][46] : where  represents a phonon mode with wavevector  and polarization  (3)     where  refers to atoms in the conventional cell,    is the  component of an eigenvector in the conventional basis for phonon mode  , and   0, ′  ′ , ′′  ′′ are third-order anharmonic IFCs 46 .
Translational invariance has already been enforced in Eq. ( 3), which leads to crystal momentum conservation for three-phonon interactions 48 : 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 through 29 Where    ~ is the  component of the PTS eigenvector for atom  in a single layer. ~ represents a phonon mode with wavevector q, phase integer  , and polarization j (as determined by diagonalization of Eq. 1).Combining Eqs. ( 3) and ( 5) and replacing atom notation () with (ℎ, ) and phonon state notation  (, ) with  ~ (, , ), as with Eq. 1, we have 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 Note that the unit cell  ′ ( ′′ ) will be shifted if ℎ ′ < ℎ (ℎ ′′ < ℎ).The layer phase relations in Eq. ( 6) are also unchanged if shifted by h: Rearranging Eq. 8 and using Eq. 4 lead to the following constraint Thus, where ⁄ is an integer (see Methods for further discussion of the relationship of  and   ).Since ℎ can be any integer, the expression in the bracket must be zero Equation ( 11) reveals another conservation rule in terms of integer  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   ~,± ~′,− ~′′ (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 Since the terms in the summation are independent of h, the summation over h gives a factor N and Eq. ( 12) becomes 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 rates 48 for a representative phonon mode at the Γ point with (ω=2.61THz, j=4) as a function of the frequency of one of the interacting phonons (  ′ ) calculated from conventional dynamics.Here, PTS dynamics was used to identify the integer  for each of the three phonon modes involved in the interactions to determine whether conservation of l occurred or not: ∆ =  ±  ′ − ( ′′ +   ) = 0 or ∆ ≠ 0. In Fig. 3b, the transition rates for ∆ ≠ 0 are numerically zero compared to those for ∆ = 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.The thermal conductivity matrix from the PTS dynamic method is calculated as where   ~ is the volumetric specific heat for phonon mode  ~ with wavevector q, polarization j, and phase integer l,  is the phonon velocity vector, and  is the phonon lifetime due to intrinsic threephonon 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 method 7,46,[48][49][50] .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 method 42,43 , 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 pressure 34,35 ) has isotropic thermal transport, i.e.,   =   =   (offdiagonal terms are zero for all systems considered here), and has the highest  of the three systems considered, with room temperature  = 112.39W/m-K.The large calculated  values for N2 are expected given its strong covalent bonding and light atomic masses, with density and mass similar to those of diamond, the prototypical ultrahigh thermal conductivity material 51 .We 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.
GeTe and CrCl3 have anisotropic thermal conductivity tensors with separate in-plane (  =   ) and cross-plane (  ) values; furthermore, phonon-isotope scattering is relatively important in each (see Fig. S5 in Supplemental Materials).In Fig. 4 we compare in-plane  for CrCl3 and an effective thermal conductivity (  = (2  +   ) 3 ⁄ ) for GeTe with measured values [52][53][54][55] .Natural isotopes are included in the calculations for both cases.For GeTe, the crystal orientation of the thin film samples was not given (thus   was used here for comparison) and significant electronic contributions to  are anticipated from measured resistance data.Thus, we compare our phonon calculations with both the total measured  (filled red circles) for GeTe thin films and the expected lattice contribution (hollow red circles, subtracting estimated electronic  via the Wiedemann-Franz law) 52 .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 here 55 .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., squares 53 , and triangle 54 .In-plane thermal conductivity measurements for CrCl3 are represented by blue circles 55 .

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 (Fig. S6 in Supplemental Material 41 ), and the average (over different grid densities) computational cost ratios are 0.3942±0.0151for N2 (N=2, n=4), 0.1334±0.0051for GeTe (N=3, n=2), and 0.0925±0.0054for 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).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 framework 3 , 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/N 2 in computational cost for GeTe and CrCl3 is primarily derived from two factors in anharmonic calculations.First, the conservation of integer  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.
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 ShengBTE 3 , almaBTE 4 , ALAMODE 5 , and phono3py 6 .Thus, the proposed PTS dynamic method can achieve a similar computational savings of ~1/N 2 with appropriate modifications in these applications.This PTS dynamics can be particularly efficient when applied to higher-order anharmonic scattering 56,57 .

DISCUSSION
We discussed lattice dynamics and phonon transport using an efficient dynamic method that uses primitive translational symmetry (PTS) within conventional 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.

Density functional theory
We performed density functional theory calculations using the projector augmented wave method (PAW) 58 , as implemented in the Vienna Ab-initio Simulation Package (VASP) [59][60][61] , for all the materials considered.The generalized gradient approximation, parameterized by Perdew, Burke, and Ernzerhof 62 was used for exchange-correlations.
The anharmonic IFCs were calculated by finite displacement method using the thirdorder.pypackage 3 .For GeTe, a 300-atom 5x5x2 supercell, Γ-only k-mesh, and a cutoff distance up to the 5 th nearest neighbors (NN) were used.For solid N2, a 216-atom 3x3x3 supercell, Γ-centered 3x3x3 k-mesh, and a cutoff distance up to the 3 rd 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 7 th NN in the plane and include most of the atoms in adjacent layers across the plane.A 192atom 2x2x2 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 Rearranging Eq. ( 15) and comparing terms leads to Such relationship is shown through phonon dispersion across two BZs for GeTe in Fig. S1.

Thermal conductivity calculations
Numerical cost of the thermal conductivity calculations was computed for a variety of reciprocal space integration mesh densities for each system (Fig. S6 in Supplemental Material).
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 package 3 .
equality.b Individual anharmonic phonon transition rates for a phonon mode at the Γ point with (ω=21.54 THz, j=7) as a function of the frequency of one of the other interacting modes in three-phonon interactions.Transition rates with ∆ ≠ 0, i.e., violating conservation of l, are numerically zero.where N is 2 and 3.The computational saving approaches the 1/N 2 as the material system becomes more complex.

Fig. 1 .
Fig. 1.Crystal structure of GeTe (R3m space group).a Two-atom primitive unit cell.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.

Fig. 2 .
Fig. 2. Comparison of phonon dispersions of GeTe from conventional dynamics (underlying black solid curves) and PTS dynamics (colored dashed-dotted curves) along high-symmetry lines and an arbitrary

Fig. 3 .
Fig. 3. Total scattering and individual transition rates for phonons in GeTe. 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.61THz, j=4) as a function of the frequency of one of the other interacting modes in three-phonon interactions.Transition rates with ∆ ≠ 0, i.e., violating conservation of l, are numerically zero.

Fig. S5 .
Fig. S5.Anisotropic thermal conductivity of a GeTe and b CrCl3.Thermal conductivities for naturally occurring samples are smaller than those for pure samples due to phonon-isotope scattering.

Fig. S6 .
Fig. S6.Ratio of computational cost for thermal conductivity calculations of N2 (black curve), GeTe (red curve), and CrCl3 (blue curve) as a function of q mesh density.The dotted lines correspond to 1/N 2 values

TABLE I .
Space groups for which PTS dynamics can be applied to conventional unit cells.