General invariance and equilibrium conditions for lattice dynamics in 1D, 2D, and 3D materials

The long-wavelength behavior of vibrational modes plays a central role in carrier transport, phonon-assisted optical properties, superconductivity, and thermomechanical and thermoelectric properties of materials. Here, we present general invariance and equilibrium conditions of the lattice potential; these allow to recover the quadratic dispersions of flexural phonons in low-dimensional materials, in agreement with the phenomenological model for long-wavelength bending modes. We also prove that for any low-dimensional material the bending modes can have a purely out-of-plane polarization in the vacuum direction and a quadratic dispersion in the long-wavelength limit. In addition, we propose an effective approach to treat invariance conditions in crystals with non-vanishing Born effective charges where the long-range dipole-dipole interactions induce a contribution to the lattice potential and stress tensor. Our approach is successfully applied to the phonon dispersions of 158 two-dimensional materials, highlighting its critical relevance in the study of phonon-mediated properties of low-dimensional materials.


INTRODUCTION
The vibrations of atoms inside crystalline solids are fundamental to many aspects of condensed matter physics and materials science. The lattice vibrations play a crucial role in a wealth of phenomena including phase transitions 1 , phonon-mediated superconductivity 2,3 , and electronic 4 and thermal transport 5,6 . The lattice-dynamical theory of solids was established at the beginning of the twentieth century with seminal works from Einstein, Born, von Kármán and Debye [7][8][9][10] . With advances in density-functional theory (DFT) and density-functional perturbation theory (DFPT) 11,12 , the parameter-free calculations of phonon dispersions from first principles have become a routine task implemented in many software packages [13][14][15][16][17][18][19] . However, there still exist some controversies in the calculated phonon dispersions of low-dimensional (LD) materials, where the longwavelength behavior of flexural phonons are not fully understood. First-principles calculations performed by different groups showed either a linear or quadratic dispersion relation for the flexural acoustic (ZA) modes in two-dimensional (2D) materials and two bending modes in one-dimensional (1D) materials (see the detailed discussion in Ref. 20). It has been recently argued 19-21 that the lack of rotational invariance in the harmonic lattice Hamiltonian with non-vanishing external stress from DFT calculations is responsible for the linear dispersion of the ZA mode; it becomes quadratic after including rotational invariance and vanishing stress condition. Like translational invariance, there are various reasons * changpeng.lin@epfl.ch for the broken rotational invariance in the interatomic force constants (IFCs) obtained from DFT and DFPT; these include incomplete basis sets, insufficient Brillouin zone sampling or numerical stability of algorithms. In contrast, IFCs obtained from semi-empirical force fields or models do not typically suffer from such issues, since they can satisfy translational and rotational invariances by construction [22][23][24][25] .
Translational and rotational invariances drive the conservation of total momentum and angular momentum in crystals 26,27 (Noether's theorem 28 ); these are known as Born-Huang invariance conditions 10,29 . In the 1960s, several works from Keating 30,31 , Gazis and Wallis 32 questioned the scope of such conditions and derived additional conditions for a harmonic lattice to be rotationally invariant. From the requirement of unchanged potential energy under any arbitrary rotation operation, Gaizs and Wallis 32 derived additional constraints for the rotational invariance of a harmonic lattice. In a similar way, through the potential energy expressed in terms of the rotation parameter up to the second order, Keating 30,31 also formulated a new set of equations for harmonic IFCs of the simple cubic lattice. However, Sarkar and Sengupta 33 demonstrated that these are both equivalent to the Born-Huang invariance conditions which are complete. Finally, Pick, Cohen, and Martin 34 derived the connection between the microscopic theory and the phenomenological approach, and the microscopic acoustic sum rule for bulk crystals.
In this work, we investigate important aspects of the Born-Huang invariance conditions as well as the so-called equilibrium conditions (i.e. Huang conditions) in determining the phonon dispersions of materials and their consequence for LD materials. We introduce the polar Born-Huang invariances and polar Huang conditions for harmonic IFCs in infrared-active solids where the existing long-range forces from the dipole-dipole interactions need to be treated separately. Indeed, Huang had already demonstrated this in the case of bulk ionic crystals 35,36 , in which the long-range Coulomb interactions need a special treatment. However, this condition seems to have been omitted in recent investigations on the rotational invariance of 2D materials [19][20][21] . The detailed implementation of these conditions for correcting the IFCs either from real-space small displacements or reciprocal-space DFPT calculations are then presented. Last, based on the secular equation for the long-wavelength vibrations, we demonstrate that the bending modes in any LD material should have a purely out-of-plane polarization in the vacuum direction and a quadratic dispersion in the longwavelength limit. We apply the present approach to silicon, triclinic calcium triphosphide (CaP 3 ), graphene, 2D molybdenum disulfide (MoS 2 ), 1D single-wall (8,0) carbon nanotube (CNT) and (4,4) boron nitride nanotube (BNNT), and find that the fulfillment of the rotational invariance and vanishing stress condition is crucial in LD materials. We also verify that these conditions have negligible effects in 3D materials, as expected. Finally, these invariance conditions are applied to correct the phonon dispersions of 245 candidate 2D materials from the highthroughput database of Ref. 37. In that study, the authors used a simple sum rule, resulting in 187 phonon dispersions presenting soft modes or an incorrect linear dispersion of the ZA mode. With the rotational invariance and equilibrium conditions correctly applied, the phonon dispersions of 158 materials become real and display a quadratic ZA branch in the long-wavelength limit. The remaining 87 materials are either dynamically unstable (54) or require tighter numerical convergence (33).

RESULTS AND DISCUSSION
Born-Huang invariance conditions The vibrational properties of solids can be determined by the nuclear Hamiltonian constructed based on the second-order IFCs. Neglecting higher-order anharmonic contributions, a Taylor expansion of the Born-Oppenheimer potential energy E in the atomic displacements up to the second order gives: where E 0 is the total potential energy of the chosen equilibrium reference structure, and Φ κα = ∂E/∂u κα and Φ κα,κ β = ∂ 2 E/∂u κα ∂u κ β are the corresponding firstand second-order IFCs, respectively, with atomic indices κ and κ in the Born-von Kármán supercell and the Cartesian directions α and β; the atom index κ ≡ {l, κ} denotes the atomic index κ in the unit cell l, and u κα is the corresponding atomic displacement from the equilibrium position. Apart from the space group symmetry of the given crystal type and permutation symmetry, these IFCs are also required to satisfy both the global translational and rotational invariances due to the conservation of total crystal and angular momenta. The Born-Huang invariance conditions provide the acoustic sum rules for both the translational and rotational invariances of a harmonic lattice. These constraint equations on the first-and second-order IFCs read 27 : where τ κα is the equilibrium position of atom κ in the Cartesian direction α and δ αβ represents the Kronecker delta. In these expressions, Equations (2) and (3) are the acoustic sum rules for the translational and rotational invariances of the first-order IFCs, while Equations (4) and (5) are their counterparts for the second-order IFCs.
When atoms are at their equilibrium positions, this set of constraint equations can be derived by adding a global translation or an infinitesimal rotation to all atoms in the system, which results in an unchanged lattice potential and zero forces acting on atoms. These acoustic sum rules correspond to the conservation of total momentum and angular momentum in the crystal. In general, the rotational sum rules link the (n + 1)-th order to the n-th order of the Taylor-expanded potential 27 . For instance, Equation (5) links the second-order IFCs to the firstorder. Equilibrium conditions For an infinite crystal, the lattice sites are perfectly periodic throughout the whole bulk region, and this requirement has two important consequences. One is the total force on each atom must be zero, which means that the first-order IFCs in Equation (1) vanish at the equilibrium. This condition gives rise to the Born-Huang invariances imposed on the second-order IFCs as being rotationally invariant: The other consequence is the vanishing of the stress tensor that governs the equilibrium volume of the unit cell, eliminating any surface effect as an infinite lattice. The stress tensor for a finite lattice is given by 27,33 where Ω is the volume of a finite lattice, and the firstorder IFCs Φ κα account for the surface forces and vanish in the bulk region 27 . In order to derive the stress tensor for an infinite lattice, we multiply the rotational invariance, Equation (5), by the atomic position τ κδ and sum over κ. Then, using Equation (7) we arrive at where we have introduced Since T αδ,βγ might not be symmetric near the surface, we define With the help of the translational invariance, Equation (4), we can further express Equation (11) as where τ κκ γ ≡ τ κ γ − τ κγ is the distance between two atoms along the Cartesian direction γ. By construction, the symmetrized tensor has the following symmetry relations T sym αβ,γδ = T sym βα,γδ = T sym αβ,δγ .
To derive the stress tensor for an infinite crystal, we further subtract both sides of Equation (10) by the quantity Ωσ γδ δ αβ to obtain T sym αβ,γδ − Ωσ γδ δ αβ = 1 2 T αγ,βδ − Ωσ γδ δ αβ + T αδ,βγ − Ωσ δγ δ αβ , (14) where the symmetry of the stress tensor σ γδ = σ δγ is used. From Equations (8) and (9), it can be noted that the right-hand side of Equation (14) is invariant by interchanging α with δ and also β with γ. Therefore, the left-hand side of Equation (14) has the following symmetry which is the expression of the stress tensor for infinite crystals 27 . At equilibrium, the stress tensor vanishes, which gives the Huang conditions T sym αβ,γδ = T sym γδ,αβ : In general, there are 36 Huang conditions due to the fourth-rank tensor T sym αβ,γδ , which can be reduced to 15 by exploiting the symmetry relations in Equation (13). For the anisotropic components of stress tensor, they become for α = β. We emphasize here that the Huang conditions are not the constraints for a potential to be rotationally invariant but the ones for a vanishing stress tensor. The set of Equations (4), (6) and (16) constitute the complete description of invariance conditions for the lattice dynamics of crystalline solids at the equilibrium. Leaving the translational invariance aside, we focus on the rotational invariance, Equation (6), and the equilibrium conditions, Equation (16), which we collectively refers to as invariance conditions.
Invariance conditions to interatomic force constants In practice, the real-space IFCs can be obtained by fitting the DFT forces using the small displacement method 13 or a Fourier transformation within DFPT 11,12 . The IFCs are then transformed into reciprocal space on dense momentum grids, giving rise to the dynamical matrices. This can be achieved by a Fourier interpolation of the IFCs.
The small displacement method fits the IFCs in the Taylor-expanded interatomic forces: This can be recast in matrix form as a linear regression problem 38 where A is the displacement matrix, Φ is an IFC vector, and C is null space constructed from all symmetry constraints on the IFCs. The constraints contain those from space group symmetry, permutation symmetry and the invariance conditions discussed so far. Through the null space C, the IFCs Φ can be projected into the independent IFC parameters φ. Then, a least-square or compressive sensing technique 38 can be applied to determine the unknown parameters φ, and the IFCs Φ are further obtained via Φ = C · φ.
In DFPT calculations, the IFCs are corrected to make them fulfill the invariance conditions and then Fourier interpolated. Two approaches exist to correct the IFCs. The first one consists in projecting the IFCs onto a subspace, spanned by the vector basis of the invariance conditions and obtaining the solution by nearest distance minimization 39 . In order to get the projection Φ p of the IFCs in such a subspace, the orthogonal basis set {g i } of the subspace needs to be built from the constraint matrix of invariance conditions and yields which is the nearest distance between the original and the corrected IFCs. The corrected IFCs Φ c that satisfy the invariance conditions are obtained via Φ c = Φ − Φ p . We have implemented the invariance conditions through this optimal projection of IFCs in the Quantum ESPRESSO distribution 15,40 . The second approach adds a correction to the IFCs and then minimizes the correction to remain as close as possible to the original IFCs. We define the deviation from the invariance conditions as where I is the constraint matrix of the invariance conditions. The null space C introduced in Equation (22) ensures that any correction made will not break other symmetries, such as space group symmetry. The prime on null space indicates that C is evaluated excluding the invariance conditions that we want to correct. We apply a small correction ∆φ such that which means that I · C · ∆φ = −d. We can solve this problem efficiently by utilizing ridge regression 41 with the help of a penalty term µ ∆φ 2 2 where . . . 2 denotes the 2 -norm and µ is a hyperparameter with values usually in the range 10 −2 to 10 −6 , which ensures that the correction ∆φ on the IFCs is as small as possible: A similar implementation that applies ridge regression to impose the invariance conditions on the IFCs can be found in the hiPhive software 19 . In the case of crystals, the periodic-boundary conditions need to be taken into account in order to impose the invariance conditions properly, and this can be tackled easily by constructing the Wigner-Seitz supercell. Thanks to the periodicity of the lattice and translational invariance, all the interactions are considered between the atoms in the chosen central unit cell and the other atoms in the supercell. To construct the Wigner-Seitz supercell, their relative distances are recalculated as the nearest periodic image. In the case of equal distances among their periodic images, a weight inversely proportional to its cardinality is given. As a result, the sum of the weights is unity because all of the atom pairs in the supercell contribute uniquely to the potential energy. This approach is similar to the one implemented in the EPW code 42,43 to improve the accuracy of the interpolated dynamical matrix, which optimizes atomic pairs to ensure that the distance between each two of atomic po-sitions are within a cutoff radius.
Fourier interpolation of interatomic force constants The dynamical matrix D κα,κ β (q) is computed via a Fourier interpolation of the real-space IFCs as: where q is a phonon wavevector defined on an arbitrary grid in reciprocal space. Due to the transitional invariance, the IFCs Φ κα,κ β in crystalline solids depend only on the relative position between the two atoms κ and κ . Therefore, it is convenient to consider only the IFCs Φ κα,κ β (R) between the atom labelled by κ within the reference unit cell 0 and the atom labelled by κ within another unit cell R.
However, in the case of infrared-active materials, the atomic vibrations generate long-range Coulomb interactions that can deteriorate the quality of the Fourier interpolations 11,12,44 . A common strategy 12 to overcome this problem is to split the dynamical matrix into a shortrange (S) and long-range (L) component as: such that the short-range part is a smooth function of q. If we consider only dipole-dipole interactions 12,45 , the long-range part takes the following form: which is non-analytic at q = 0 and induces the directiondependent LO-TO splitting. To deal with the q = 0 term exactly, the Ewald summation technique 11,12 is adopted in practice with a sum over reciprocal G vectors. In this expression, e is the elementary charge unit, Ω 0 is the volume of unit cell, Z * κ is the Born effective charge tensor of the atom κ within the unit cell, and W c (q) is the screened Coulomb potential which depends on the dimensionality of the system as 44,46,47 where is the ion-clamped dielectric tensor of 3D system, r eff ≈ b/2 is the in-plane matrix of effective screening length, b is the thickness of the 2D material, 1D is the 1D dielectric constant, t is the effective screening radius, I n (x) and K n (x) are the nth-order modified cylindrical Bessel functions, and G mn pq is the Meijer G-function 48 . The 1D long-range formula, Equation (30), has been recently derived and more details can be found in Ref. 47. The long-range part is substracted from the dynamical matrix on the coarse grid, the remaining short-range part is Fourier transformed into real space to get the IFCs.
Inspired by this approach, we formulate the polar Born-Huang invariances and polar Huang conditions by separating the total IFCs into short-range and long-range IFCs, and Since the Φ L are analytically defined, we apply the corrections due to the invariance conditions only on Φ S . Due to the presence of the long-range terms, we have to solve an inhomogeneous set of linear equations which is here achieved by using a ridge regression method 41 . We then Fourier interpolate the resulting corrected shortrange IFCs at arbitrary momentum points and add back the analytic long-range part. This approach is an alternative to applying the Born-Huang invariances, Equation (6), and Huang conditions, Equation (16), on the total real-space IFCs. We have verified in Fig. 4 that both methods give the same interpolated 2D phonon band structures. We also reiterate that we have neglected dynamical quadrupoles in the long-range treatments which have been recently shown to be important in some cases 43,45,46,49 and have been left for future works. Finally, in the case of low-symmetry infrared-active solids, we have observed that the total dynamical matrices are not guaranteed to be Hermitian. As observed in Ref. 50, the dipole-dipole IFC ansatz (e.g. Equation (4) of Ref. 51 or Equation (70) of Ref. 12) which is used in the analytic form of the long-range dynamical matrix does not guarantee Hermiticity for a general crystal. Zhou et al. 50 addressed this issue by adding a real-space short-range correction to the IFC matrix. Here, for simplicity, we directly symmetrize the long-range part of the on-site dynamical matrix by its conjugate transpose and Fourier transform both the short-range and long-range parts into real space. We then apply the invariance conditions on the total IFCs and subsequently substract the symmetrized long-range in real space. Last, we Fourier interpolate the short-range IFCs to arbitrary momentum and add the long-range part symmetrized in the same way.
Bending modes in 2D materials The lattice dynamics of crystalline solids is governed by the secular equation 1,10 [ω ν (q)] 2 e ν,κα (q) = κ β where m κ is the mass of the κ-th atom within the unit cell. The eigenvector e ν,κα (q) is obtained by diagonalizing the mass-scaled dynamical matrix D κα,κ β (q)/ √ m κ m κ with the corresponding eigenfrequency ω ν (q) for the phonon mode ν at momentum q, and ν ranges from 1 to 3n with n atoms in the unit cell.
To analyze the dispersion relation of bending modes in the long-wavelength limit (q → 0), we perform a Taylor expansion of Equation (33) in terms of q. The resulting long-wavelength equation for lattice vibrations is 10,21 : where ω is an arbitrary polarization vector which is introduced via the solution for the zero-order equation of the Taylor expansion, e ν,κα (q) = √ m κ u να . The T sym αβ,γδ defined in Equation (12) can be rewritten in term of a sum over primitive cells as and T int αγ,βδ is the contribution associated with the relaxation of internal coordinates in response to external stress fields 10 : with the ν × ν matrix Γ κλ,κ µ denoting the inverse of the mass-scaled IFC matrix R Φ κλ,κ µ (R)/ √ m κ m κ . Unfortunately, the IFC matrix is singular because of the acoustic sum rules for translational invariance and the Γ κλ,κ µ is thus introduced as where the pseudoinverseΓ κλ,κ µ is calculated and satisfies the relation 10 Furthermore, the tensor T int αγ,βδ satisfies the following symmetry relations 10,36 : which result from the rotational invariance, Equation (6), and the fact that Γ κλ,κ µ = Γ κ µ,κλ . For the ZA bending modes of 2D materials in the longwavelength limit, if they were not purely polarized along the vacuum direction z and coupled with the in-plane vibrations, the leading linear term in the ZA dispersion would read: where the phonon band index ν ≡ ZA is omitted for clarity. By using the Huang conditions and the symmetry of T int αγ,βδ , the identity T sym holds. As demonstrated in Supplementary Note 1 52 , we have T sym γδ,zβ + T int γz,δβ = 0 even in case of non-purely outof-plane modes if the rotational invariance of the lattice potential is fullfilled and the system is free of stress, leading to the vanishing linear dispersion term. In the long-wavelength limit, the dispersion of ZA modes thus becomes quadratic with a purely out-of-plane polarization. Having proven this, we can therefore define the bending conditions as: which can be imposed in case of numerical noise in the calculations to ensure a quadratic dispersion of the flexural modes. More generally, these conditions for recovering quadratic flexural modes are also applicable to 1D materials, where there are two acoustic bending modes polarizing in two vacuum directions. In 1D materials with periodicity in the z direction, the bending conditions thus become T sym γδ,xβ = −T int γx,δβ and T sym γδ,yβ = −T int γy,δβ for polarizations in the x and y directions, respectively. We note that if LD materials are subjected to external stress fields, the linear term will then dominate the phonon dispersion in the long-wavelength limit. In bulk materials, no such bending conditions exist due to the periodicity in all three directions.
Application to bulk materials The influence of the rotational invariance on the lattice dynamics of bulk materials is almost unexplored. We show the phonon dispersion of silicon in Fig. 1 and find that both the Born- Influence of invariance conditions on the phonon dispersion of the infrared-active bulk CaP3. TI denotes the correction of translational invariance on the IFCs, while BH and H represent the corrections by Born-Huang rotational invariance and Huang conditions, respectively.
Huang invariances (red line) alone and the combination of Born-Huang invariances and Huang conditions (blue line) have little effect, compared with the phonon dispersion obtained from the IFCs without the aforementioned conditions imposed (black dashed line). This result is not surprising given that the IFCs of high-symmetry bulk crystals with space inversion automatically fulfill rotational invariance 53 .
As an example of infrared-active bulk materials, we calculate the phonon dispersion of triclinic CaP 3 shown in Fig. 2, with the LO-TO splitting included. The crystal symmetry of CaP 3 with the space group of P1 is low: only the identity and the spatial inversion operations. For a 2 × 2 × 2 supercell with 64 atoms, there are 160 unique atomic pairs which means 1440 IFC components. The total number of independent IFC parameters in CaP 3 is 1215 after considering the space group symmetry, the permutation symmetry and the translational invariance. When we add the Born-Huang invariances on the IFCs, the number of independent IFCs is reduced to 1179. The number of independent IFC parameters can further decrease to 1164 if the Huang conditions are included as symmetry constraints. Besides, there is a large deviation from these two invariance conditions reflected by the Frobenius norm of d in Equation (22); it vanishes with a magnitude of 10 −6 after the invariance conditions are enforced, which indicates that the Born-Huang invariances and Huang conditions are now satisfied. However, the phonon dispersions of CaP 3 in Fig. 2 are almost unchanged. There is only a small modification of the lowest-lying transverse acoustic (TA) branch. Hence, we conclude that the violation of the rotational invariance does not affect the phonon dispersion of bulk materials at equilibrium, and the corrections of the invariance conditions on the IFCs of bulk materials are not important, even for low-symmetry crystals. This can be rationalized from the fact that bulk materials are periodic along all directions and thus do not have any rigid rotation. Nonetheless, for the elastic tensor c ij (in Voigt notation) determined directly from the second-order IFCs, the broken symmetry relation (i.e. c ij = c ji ) of low-symmetry solids can be observed due to the violation of the rotational invariance and equilibrium conditions 53 .
Application to 2D materials As revealed by a few recent studies 19-21 , the violation of the rotational invariance and vanishing stress condition are expected to have a remarkable effect on the lattice-dynamical properties of 2D materials, which can lead to an unphysical dispersion relation of the out-of-plane ZA branch. The influence of such two invariance conditions on the phonon dispersions of infrared-inactive 2D materials is illustrated in Fig. 3, where we take graphene as a representative material. In general, the lack of these two invariance conditions can result in a ZA mode either displaying imaginary frequencies or a non-quadratic dispersion relation in the longwavelength limit. We confirm these findings by fitting the IFCs obtained from randomly displaced configurations. As can be seen in the inset of Fig. 3, the pristine phonon dispersion (black dashed line) of graphene without the Born-Huang invariances and Huang conditions exhibits a linear ZA branch close to the zone center. When we add the Born-Huang invariances to make the IFCs satisfy the global rotational invariance, the result (blue line) remains unchanged with respect to the pristine one. This happens because the crystal symmetry of graphene is high and the rotational invariance of its lattice potential is automatically satisfied. In contrast, only the red and yellow lines, where the additional Huang conditions are fulfilled, show the expected parabolic dispersion of the ZA branch around the Γ point. More generally, both the Born-Huang invariances and Huang conditions must be satisfied to observe the physically correct ZA phonons with quadratic dispersions in the long-wavelength limit. To illustrate this, we consider an MoS 2 monolayer (see Supplementary Fig. 1 52 ), which lacks the additional inversion symmetry. We see that fulfilling the vanishing stress condition alone is not sufficient for obtaining a parabolic ZA branch, and that the rotational invariance is required as well. We therefore deduce that residual stress remains after structural optimization and that this equilibrium condition enforces strict null stress to recover the quadratic dispersion relation of the ZA phonons. In particular, there is still a non-negligible stress in the vacuum direction even after a high-quality structural optimization with a Coulomb cutoff for 2D materials 44,54 , owing to the existing weak layer interactions from the periodic images. This observation is consistent with the dispersion model of long-wavelength bending modes in Equation (40) which indicates that any external stress will give rise to a non-vanishing linear dispersion term and hence suppresses the next-leading quadratic dispersion. Overall, both the Born-Huang rotational invariance and Huang conditions (vanishing stress condition) are of great importance for 2D materials to have the expected behavior of the ZA branch in the long-wavelength limit.
Most 2D insulators and semiconductors are infraredactive with non-vanishing Born effective charges and long-range interactions. In such cases, the long-range Coulomb interactions provide additional contributions to the IFCs 35,36 . The enforcement of the invariance conditions for infrared-active materials should therefore be performed with caution. We take the MoS 2 monolayer as an example to showcase this special treatment. In 2D materials, the LO-TO splitting vanishes at the zone center, and the LO dispersions become linear in q with a finite slope 44,54 . As shown in Fig. 4, the phonon dispersions of MoS 2 monolayer calculated directly from the pristine IFCs from DFPT (black dashed line) exhibit small imaginary frequencies in the ZA branch close to the zone center, as a consequence of the broken invariance conditions. If we add the corrections according to the invariance conditions from Equations (6) and (16) on the short-range IFCs only (blue line), there still exist some imaginary frequencies around the Γ point, belonging to the ZA phonons. This can be understood by realizing that the long-range IFCs indeed contribute to the stress tensor and the rotational invariance of the lattice potential. When the long-range interactions are added back analytically, the resulting total dynamical matrices in reciprocal space will not fulfill the Born-Huang invariances and Huang conditions. This is further confirmed by looking the yellow line in Fig. 4, where the non-analytic correction corresponding to the long-range interactions are removed, and the result presents a stable ZA branch in the parabolic shape near the Γ point (i.e. the shortrange IFCs satisfying the invariance conditions). To deal with the invariance conditions properly in infrared-active solids, we need to consider Equation (31) for the polar Born-Huang invariances and Equation (32) for the polar Huang conditions to make the total IFCs meet the invariance conditions (green line). The long-range IFCs contributed from the dipole-dipole interactions are considered when imposing the invariance conditions on the short-range IFCs only, and an analytic long-range correction is further introduced back in reciprocal space to recover the LO-TO splitting. As indicated by the red line in Fig. 4, an equivalent approach is to apply the normal invariance conditions on the total IFCs, substract the long-range part obtained by Fourier transformation, and Fourier interpolate the resulting short-range IFCs with an analytic long-range contribution added back. Both approaches yield the correct ZA branch and the 2D LO-TO splitting shown in the inset of Fig. 4. Apart from the effects on the LO-TO splitting of optical phonons, we observe that the long-range dipole-dipole interactions also modify the ZA modes by lowering their energies when comparing the results in yellow and red lines. In the rest of the manuscript, we always use the last approach (red line) in infrared-active materials.
Application to 1D materials Another important class of LD materials are 1D materials. The effects of the invariance conditions on their lattice-dynamical properties have not yet received much attention. We choose the single-wall (8,0) CNT and (4,4) BNNT as representative infrared-inactive and -active 1D cases, respectively. We illustrate in Fig. 5 the corrections from the invariance conditions on the phonon dispersion of (8,8) CNT, where we can see that the results with only the translational invariance applied yield only three vanishing acoustic frequencies at the Brillouin zone center. In 1D systems with more than one atom in the unit cell, four acoustic branches are expected and correspond to the three global translational degrees of freedom in the three cartesian directions and one rotational degree of freedom along the 1D axis. This is correctly recovered when applying the Born-Huang invariances as shown with a red line in Fig. 5. As shown previously 25,39,55,56 , the frequency of the twisting acoustic mode due to the global rotation along the nanotube axis can become finite at the Γ point (around 1.0 THz in our case) if the Born-Huang rotational invariance is not enforced. Also, the dispersion relation of the two acoustic branches corresponding to the bending motions becomes parabolic in the long-wavelength limit after adding the corrections from the Born-Huang rotational invariance. The Huang conditions for vanishing external stress seem to have no influence on the case of single-wall (8,8)  Corrections due to the invariance conditions on the phonon dispersion of an infrared-inactive 1D single-wall (8,0) carbon nanotube. TI denotes the correction of translational invariance on the IFCs, while BH and H represent the corrections by Born-Huang rotational invariance and Huang conditions, respectively. The inset shows a zoom-in result around the Γ point, with four acoustic branches labelled as longitudinal acoustic (LA), twisting acoustic (TWA) and flexural acoustic (ZA, doubly degenerate) modes, respectively. The blue line coincides with the red one, which indicates the Huang conditions have no effect and the external stress is already vanishing after structural optimization. tight structural optimization. In general, both the Born-Huang invariances and Huang conditions should be fulfilled in order to find the correct dispersions of two bending modes in 1D materials. It is the case for the infraredactive (4,4) BNNT shown in Fig. 6, where the two bending modes with quadratic dispersions are only recovered when both invariance conditions are imposed. Furthermore, there are long-range Coulomb interactions due to the non-vanishing Born effective charges in BNNT. The quadratic phonon dispersions of two bending modes and the emergence of twisting acoustic modes with vanishing frequency at the zone center can only be observed with the polar invariance conditions introduced in Equations (31) and (32). Overall, our proposed invariance scheme can be applied to any material ranging from bulk to 2D and to 1D systems to recover exact latticedynamical properties.
2D materials database A recent high-throughput screening by Mounet et al. 37 on experimentally known layered compounds has identified 1825 materials that can be easily or potentially exfoliated. The vibrational properties of 245 of these 2D materials were calculated and explored, with all the data available on the Materials Cloud 57 . However, a large proportion of the calculated phonon dispersions showed imaginary frequencies, which suggests that the predicted 2D materials would be dynamically unstable. Since most of the unstable phonon branches are ZA modes in the zone center region, we apply here our invariance scheme to see if these can be stabilized.
We select four representative 2D materials from the database, namely FeO 2 , PrOI, InSe and BiOCl, and report their phonon dispersions in Fig. 7. As can be seen, after the invariance conditions are enforced, all of the four phonon dispersions become stable with a quadratic ZA mode in the long-wavelength limit. Besides, FeO 2 and PrOI are metals, while InSe and BiOCl are inrared-active semiconductors with non-zero Born effective charges. The results in Fig. 7 highlight the effectiveness of the invariance conditions described in this work. Finally, our simulations shed a new light on the origin of imaginary frequencies existing in the 2D material database and shows that in most cases this was due to a violation of some invariance conditions. We therefore expect that most 2D materials reported are stable and will be synthesized in the future. ). There are mainly two reasons why the correction of the invariance conditions fails to recover a stable and quadratic ZA mode at the zone center. One of them is the numerical accuracy of the structure optimization and DFPT calculation, which might not be sufficient to produce accurate real-space IFCs. The other is more subtle: one might need to expand Equation (1) also in terms of strain-electric field 59 and include the treatment of higherorder multipolar interactions 43,45,46,49 , which we neglect here for simplicity. To reveal the first problem, we have selected a few of cases including HfS 2 , NbF 4 , CoO 2 , TaS 2 and C where the structure optimization and subsequent DFPT were repeated with more stringent calculation parameters. The corresponding phonon dispersions based on these new IFCs with the invariance conditions applied become stable and display a quadratic ZA mode in the long-wavelength limit (see the results with the legend this work in Supplementary Note 4 52 ). In a few cases involving magnetic 2D materials, the self-consistent solution found might not be the lowest-energy ground state. We did not investigate these cases further. For the second case, we observe that black phosphorus and arsenene have a non-vanishing yz components of the Born effective charge tensor, resulting in small soft modes due to the neglect of electric field in calculating IFCs and higherorder long-range electrostatic fields in the interpolation of IFCs, which could be important in the presence of piezoelectricity or flexoelectricity 45,59-61 .
In conclusion, we have presented a systematic study on the invariance and equilibrium conditions in the lattice dynamics of crystalline solids, as well as its implementation in the framework of modern first-principles simulations. The results highlight the importance of enforcing both the Born-Huang invariances and Huang conditions in order to recover the quadratic dispersion relation of the ZA phonons in LD materials around the Brillouin zone center. This is in agreement with the theoretical model discussed in Ref. 21 and this work for bending waves, where the first-order dispersion vanishes as the bending conditions regardless of the crystal type of LD solids. As expected, such two invariance conditions are found to have negligible influence on the lattice dynamics of bulk materials, even for low-symmetry crystals. Besides, we have introduced polar Born-Huang invariances and polar Huang conditions, which allow to separately treat the short-range and long-range IFCs for infrared-active solids. This scheme can work effectively in infrared-active 2D materials where the long-range dipole-dipole interactions contribute to the fulfillment of the invariance conditions. To demonstrate the effectiveness of the method developed, we apply this invariance scheme to over two hundred 2D material candidates from a high-throughput database, most of which originally displayed an unstable phonon dispersion with the imaginary ZA modes close to the Γ point. With four types of representative 2D materials from the database illustrated here, we demonstrate that most of their phonon dispersions become stable with a physically quadratic ZA branch in the long-wavelength limit, after the invariance conditions are imposed. We believe the proposed rotational invariance and equilibrium conditions will have wide applications in the field of LD materials, benefiting future lattice-dynamical studies of novel materials, and the calculations of their mechanical, thermodynamic, vibrational, and spectroscopic properties.

METHODS
The DFT and DFPT calculations are performed with Quantum ESPRESSO 15,40 and the v1.1 SSSP efficiency pseudopotentials' library 62 . To describe the exchange and correlation of electrons, we employ the Perdew, Burke, and Ernzerhof's parametrization for solids (PBEsol) 63 of the generalized gradient approximation (GGA). In the case of 2D materials, the Coulomb cutoff for 2D system implemented by Sohier et al. 54 is adopted to remove the long-range interactions in the vacuum direction, which is important to obtain the correct electronic screening and LO-TO splitting in 2D materials 44 . To compute the IFCs directly in real space from the small displacement method, we use an in-house code which is interfaced with Quantum ESPRESSO. The subsequent Fourier interpolation to reciprocal space is obtained with Phonopy 14 to compute the phonon dispersions.
Silicon The optimized lattice parameter is 5.43Å, with convergence thresholds for pressure and forces smaller than 10 −2 kbar and 10 −4 Ry/Bohr, respectively. The same criteria for other materials are used unless stated. To calculate the real-space IFCs, we use a planewave kinetic energy cutoff of 60 Ry with a Γ-point grid computed on a 4 × 4 × 4 supercell and perform a random displacement with the magnitude of 0.01Å for all atoms. A least-square fit of the Taylor-expanded forces to the DFT forces is carried out to find the best solution for the IFCs.
Triclinic CaP 3 The structure is relaxed with a plane-wave kinetic energy cutoff of 30 Ry and a 10 × 10 × 10 electronic grid. Similar to the case of silicon, we build a 2 × 2 × 2 supercell, and the real-space DFT calculations based on a 2 × 2 × 2 grid are carried out for 30 displaced structures to get the entire IFCs. In addition, a DFPT calculation at q = 0 is performed to get the dielectric tensor and Born effective charges.
Graphene To simulate graphene as a 2D material, a vacuum of 20Å is set in the out-of-plane direction to remove the interaction between periodic images. The structure is relaxed using a plane-wave kinetic energy cutoff of 80 Ry and a 20 × 20 × 20 electronic grid. The optimized lattice constant of graphene is 2.46Å. To calculate the phonon dispersions, the same real-space procedure is adopted. We create a 6 × 6 × 1 supercell and the corresponding DFT calculations for 10 displaced supercells, with an electronic sampling of 3 × 3 × 3.
2D MoS 2 We choose a vacuum of 25Å to avoid the spurious interactions in the out-of-plane direction. With a plane-wave kinetic energy cutoff of 55 Ry and a 14 × 14 × 1 electronic grid, we obtain the optimized lattice constant as 3.14Å based on tighter relaxation criteria for pressure (smaller than 10 −3 kbar) and forces (smaller than 10 −5 Ry/Bohr). For the lattice-dynamical calculations, we perform a DFPT calculation on a 6 × 6 × 1 phonon grid, and the IFCs are obtained through the inverse Fourier transform of the dynamical matrices.
1D materials For the single-wall (8,0) CNT, the periodic direction is chosen along the z direction with a vacuum of 20Å in both the x and y directions, and the lattice constant of the optimized structure is 4.26Å based on the relaxation criteria of pressure (smaller than 10 −2 kbar) and force (smaller than 10 −5 Ry/Bohr). We create a 1 × 1 × 4 supercell and generated 60 displaced configurations using the same procedure as before to obtain the real-space IFCs. The DFT calculations for these displaced supercells are performed with a plane-wave kinetic energy cutoff of 80 Ry and a Γ-point electronic grid. For the infrared-active single-wall (4,4) BNNT, the IFCs are obtained by DFPT calculations as in Ref. 47, with the appropriate Coulomb cutoff for 1D systems.

DATA AVAILABILITY
The material structures, pseudopotentials, interatomic force constants and related data in order to reproduce this study can be found on the Materials Cloud Archive 58 .

CODE AVAILABILITY
The code developed and used in this work is integrated into the q2r.x and matdyn.x executables of the Quan-tum ESPRESSO distribution, which will be released in the next version of Quantum ESPRESSO. force constants of the single-wall (4,4) boron nitride nanotube. We are also grateful to Alexander Croy for providing us with additional information on Ref.   Supplemental information: General invariance and equilibrium conditions for lattice dynamics in 1D, 2D, and 3D materials As discussed in the main text of this work, the first-order linear dispersion of the bending ZA mode in 2D materials reads In the second step of Equation (1), the Huang conditions and the symmetry of the tensor T int αγ,βδ (see Equations (15) and (38) in the main text) are used. We denote the z axis as the vacuum direction for 2D materials. As a result, the z component of distance between two atoms does not depend on the unit cell position R, and one can have the bending conditions with Φ R κγ,κ δ ≡ Φ κγ,κ δ (R) and τ R κκ β ≡ τ κκ β (R). The translational invariance and the properties of Γ κλ,κ µ matrix (see Equations (4) and (37) in the main text) are used in the derivation with the rearrangement of summation indices. Therefore, the linear dispersion term ω (1) ZA (q) for the ZA mode of 2D materials vanishes, and its dispersion relation becomes quadratic in the long-wavelength limit. These bending conditions T sym γδ,zβ = −T int γz,δβ also allow the ZA branch to have a purely out-of-plane polarization in the vacuum z direction, fully decoupling from the in-plane vibrational modes.

SUPPLEMENTARY NOTE 2: SUMMARY OF THE 2D MATERIALS DATABASE
The current 2D materials database available on the Materials Could archive 1 contains 245 entries for the phonon dispersions of 2D materials 2 . However, 54 entries are dynamically unstable; the phonon spectra display remarkable imaginary frequencies not located around the Brillouin zone center, as summarized in Supplementary Table 1. After excluding these cases, we have imposed the corrections of the Born-Huang rotational invariance and Huang conditions on the original IFCs from Materials Could archive 1 and recalculated the corresponding phonon dispersions. We find that the corrections of the two invariance conditions failed for 33 entries of the remaining 191 candidates in the 2D database, which are listed out in Supplementary Table 2. By analyzing these results, we identify two main reasons that result in the failure of the invariance conditions to yield stable quadratic dispersion for the ZA branch. One of them is inaccurate IFCs in real space as obtained from DFPT calculations, due to the choice of convergence parameters and not perfectly relaxed structures. To showcase this issue, we select a few entries (i.e. HfS 2 , NbF 4 , CoO 2 , TaS 2 and C) and perform further structural optimizations, and calculate again the IFCs by DFPT. The phonon dispersions from these new DFPT calculations are presented in Supplementary Note 4 with the legend this work. In addition, for some magnetic materials the self-consistent solution found might not be the lowest-energy ground state; the neglect of higher-order multipolar interactions in the interpolation of IFCs could also be important for piezoelectric and flexoelectric materials 3,4 .