si.pdf

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

Changpeng Lin<EMAIL_ADDRESS>Theory and Simulation of Materials
(THEOS), École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne,
Switzerland  National Centre for Computational Design and Discovery of Novel
Materials (MARVEL), École Polytechnique Fédérale de Lausanne, CH-1015
Lausanne, Switzerland  Samuel Poncé  European Theoretical Spectroscopy
Facility, Institute of Condensed Matter and Nanosciences, Université
catholique de Louvain, Chemin des Étoiles 8, B-1348 Louvain-la-Neuve, Belgium
Theory and Simulation of Materials (THEOS), École Polytechnique Fédérale de
Lausanne, CH-1015 Lausanne, Switzerland  Nicola Marzari  Theory and Simulation
of Materials (THEOS), École Polytechnique Fédérale de Lausanne, CH-1015
Lausanne, Switzerland  National Centre for Computational Design and Discovery
of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, CH-1015
Lausanne, Switzerland

###### Abstract

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 long-
wavelength 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 20, 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 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 20, 19, 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 long-wavelength limit. We apply the present
approach to silicon, triclinic calcium triphosphide (CaP3), graphene, 2D
molybdenum disulfide (MoS2), 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 high-
throughput 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:

$E=E_{0}+\sum_{\varkappa\alpha}\Phi_{\varkappa\alpha}u_{\varkappa\alpha}+\frac{1}{2}\\!\\!\sum_{\varkappa\alpha,\varkappa^{\prime}\beta}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}u_{\varkappa\alpha}u_{\varkappa^{\prime}\beta},$
(1)

where $E_{0}$ is the total potential energy of the chosen equilibrium
reference structure, and $\Phi_{\varkappa\alpha}=\partial E/\partial
u_{\varkappa\alpha}$ and
$\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}=\partial^{2}E/\partial
u_{\varkappa\alpha}\partial u_{\varkappa^{\prime}\beta}$ are the corresponding
first- and second-order IFCs, respectively, with atomic indices $\varkappa$
and $\varkappa^{\prime}$ in the Born-von Kármán supercell and the Cartesian
directions $\alpha$ and $\beta$; the atom index
$\varkappa\equiv\\{l,\kappa\\}$ denotes the atomic index $\kappa$ in the unit
cell $l$, and $u_{\varkappa\alpha}$ 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:

$\displaystyle\sum_{\varkappa}\Phi_{\varkappa\alpha}$ $\displaystyle=0,$ (2)
$\displaystyle\sum_{\varkappa}\Phi_{\varkappa\alpha}\tau_{\varkappa\beta}$
$\displaystyle=\sum_{\varkappa}\Phi_{\varkappa\beta}\tau_{\varkappa\alpha},$
(3)
$\displaystyle\sum_{\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}$
$\displaystyle=0,$ (4)
$\displaystyle\sum_{\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}\tau_{\varkappa^{\prime}\gamma}+\delta_{\alpha\gamma}\Phi_{\varkappa\beta}$
$\displaystyle=\sum_{\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\gamma}\tau_{\varkappa^{\prime}\beta}+\delta_{\alpha\beta}\Phi_{\varkappa\gamma},$
(5)

where $\tau_{\varkappa\alpha}$ is the equilibrium position of atom $\varkappa$
in the Cartesian direction $\alpha$ and $\delta_{\alpha\beta}$ 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 first-order.

### 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:

$\sum_{\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}\tau_{\varkappa^{\prime}\gamma}=\sum_{\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\gamma}\tau_{\varkappa^{\prime}\beta}.$
(6)

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

$\sigma_{\alpha\beta}=\frac{1}{\Omega}\sum_{\varkappa}\Phi_{\varkappa\alpha}\tau_{\varkappa\beta},$
(7)

where $\Omega$ is the volume of a finite lattice, and the first-order IFCs
$\Phi_{\varkappa\alpha}$ 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
$\tau_{\varkappa\delta}$ and sum over $\varkappa$. Then, using Equation (7) we
arrive at

$T_{\alpha\delta,\beta\gamma}-\Omega\sigma_{\gamma\delta}\delta_{\alpha\beta}=T_{\alpha\delta,\gamma\beta}-\Omega\sigma_{\beta\delta}\delta_{\alpha\gamma},$
(8)

where we have introduced

$T_{\alpha\delta,\beta\gamma}\equiv\sum_{\varkappa\varkappa^{\prime}}\Phi_{\varkappa\alpha\varkappa^{\prime}\beta}\tau_{\varkappa\delta}\tau_{\varkappa^{\prime}\gamma}=T_{\beta\gamma,\alpha\delta}.$
(9)

Since $T_{\alpha\delta,\beta\gamma}$ might not be symmetric near the surface,
we define

$\displaystyle T^{\rm sym}_{\alpha\beta,\gamma\delta}$
$\displaystyle\equiv\frac{T_{\alpha\gamma,\beta\delta}+T_{\alpha\delta,\beta\gamma}}{2}$
(10)
$\displaystyle=\frac{1}{2}\sum_{\varkappa\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}(\tau_{\varkappa\gamma}\tau_{\varkappa^{\prime}\delta}+\tau_{\varkappa\delta}\tau_{\varkappa^{\prime}\gamma}).$
(11)

With the help of the translational invariance, Equation (4), we can further
express Equation (11) as

$\displaystyle T^{\rm sym}_{\alpha\beta,\gamma\delta}=$
$\displaystyle-\frac{1}{2}\sum_{\varkappa\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}(\tau_{\varkappa\gamma}-\tau_{\varkappa^{\prime}\gamma})(\tau_{\varkappa\delta}-\tau_{\varkappa^{\prime}\delta})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{\varkappa\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}\tau_{\varkappa\varkappa^{\prime}\gamma}\tau_{\varkappa\varkappa^{\prime}\delta},$
(12)

where
$\tau_{\varkappa\varkappa^{\prime}\gamma}\equiv\tau_{\varkappa^{\prime}\gamma}-\tau_{\varkappa\gamma}$
is the distance between two atoms along the Cartesian direction $\gamma$. By
construction, the symmetrized tensor has the following symmetry relations

$T^{\rm sym}_{\alpha\beta,\gamma\delta}=T^{\rm
sym}_{\beta\alpha,\gamma\delta}=T^{\rm sym}_{\alpha\beta,\delta\gamma}.$ (13)

To derive the stress tensor for an infinite crystal, we further subtract both
sides of Equation (10) by the quantity
$\Omega\sigma_{\gamma\delta}\delta_{\alpha\beta}$ to obtain

$T^{\rm
sym}_{\alpha\beta,\gamma\delta}-\Omega\sigma_{\gamma\delta}\delta_{\alpha\beta}=\\\
\frac{1}{2}\Big{[}T_{\alpha\gamma,\beta\delta}-\Omega\sigma_{\gamma\delta}\delta_{\alpha\beta}+T_{\alpha\delta,\beta\gamma}-\Omega\sigma_{\delta\gamma}\delta_{\alpha\beta}\Big{]},$
(14)

where the symmetry of the stress tensor
$\sigma_{\gamma\delta}=\sigma_{\delta\gamma}$ is used. From Equations (8) and
(9), it can be noted that the right-hand side of Equation (14) is invariant by
interchanging $\alpha$ with $\delta$ and also $\beta$ with $\gamma$.
Therefore, the left-hand side of Equation (14) has the following symmetry

$T^{\rm
sym}_{\alpha\beta,\gamma\delta}-\Omega\sigma_{\gamma\delta}\delta_{\alpha\beta}=T^{\rm
sym}_{\gamma\delta,\alpha\beta}-\Omega\sigma_{\alpha\beta}\delta_{\gamma\delta},$
(15)

which is the expression of the stress tensor for infinite crystals 27. At
equilibrium, the stress tensor vanishes, which gives the Huang conditions
$T^{\rm sym}_{\alpha\beta,\gamma\delta}=T^{\rm
sym}_{\gamma\delta,\alpha\beta}$:

$\sum_{\varkappa\varkappa^{\prime}}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}\tau_{\varkappa\varkappa^{\prime}\gamma}\tau_{\varkappa\varkappa^{\prime}\delta}=\sum_{\varkappa\varkappa^{\prime}}\Phi_{\varkappa\gamma,\varkappa^{\prime}\delta}\tau_{\varkappa\varkappa^{\prime}\alpha}\tau_{\varkappa\varkappa^{\prime}\beta}.$
(16)

In general, there are 36 Huang conditions due to the fourth-rank tensor
$T^{\rm sym}_{\alpha\beta,\gamma\delta}$, which can be reduced to 15 by
exploiting the symmetry relations in Equation (13). For the anisotropic
components of stress tensor, they become

$\displaystyle\sigma_{\alpha\beta}=$ $\displaystyle\frac{T^{\rm
sym}_{\beta\beta,\beta\alpha}-T^{\rm sym}_{\beta\alpha,\beta\beta}}{\Omega}$
(17) $\displaystyle\sigma_{\alpha\alpha}-\sigma_{\beta\beta}=$
$\displaystyle\frac{T^{\rm sym}_{\beta\beta,\alpha\alpha}-T^{\rm
sym}_{\alpha\alpha,\beta\beta}}{\Omega},$ (18)

for $\alpha\neq\beta$. 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
12, 11. 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:

$F_{\varkappa\alpha}=-\sum_{\varkappa^{\prime}\beta}\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}u_{\varkappa^{\prime}\beta}+\cdots.$
(19)

This can be recast in matrix form as a linear regression problem 38

$\mathbb{F}=\mathbb{A}\cdot\mathbb{\Phi}=\mathbb{A}\cdot\mathbb{C}\cdot\boldmath{\phi},$
(20)

where $\mathbb{A}$ is the displacement matrix, $\mathbb{\Phi}$ is an IFC
vector, and $\mathbb{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 $\mathbb{C}$, the IFCs $\mathbb{\Phi}$ can be projected
into the independent IFC parameters $\boldsymbol{\phi}$. Then, a least-square
or compressive sensing technique 38 can be applied to determine the unknown
parameters $\boldsymbol{\phi}$, and the IFCs $\mathbb{\Phi}$ are further
obtained via $\mathbb{\Phi}=\mathbb{C}\cdot\boldsymbol{\phi}$.

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 $\mathbb{\Phi}^{\rm 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

$\mathbb{\Phi}^{\rm p}=\sum_{i}(\mathbb{\Phi}\cdot g_{i})g_{i},$ (21)

which is the nearest distance between the original and the corrected IFCs. The
corrected IFCs $\mathbb{\Phi}^{\rm c}$ that satisfy the invariance conditions
are obtained via $\mathbb{\Phi}^{\rm c}=\mathbb{\Phi}-\mathbb{\Phi}^{\rm p}$.
We have implemented the invariance conditions through this optimal projection
of IFCs in the Quantum ESPRESSO distribution 40, 15.

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

$\mathbb{d}\equiv\mathbb{I}\cdot\mathbb{C^{\prime}}\cdot\boldmath{\phi},$ (22)

where $\mathbb{I}$ is the constraint matrix of the invariance conditions. The
null space $\mathbb{C^{\prime}}$ 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 $\mathbb{C^{\prime}}$ is evaluated
excluding the invariance conditions that we want to correct. We apply a small
correction $\Delta\boldsymbol{\phi}$ such that

$\mathbb{I}\cdot\mathbb{C^{\prime}}\cdot(\boldsymbol{\phi}+\Delta\boldsymbol{\phi})=\mathbb{0},$
(23)

which means that
$\mathbb{I}\cdot\mathbb{C^{\prime}}\cdot\Delta\boldsymbol{\phi}=\mathbb{-d}$.
We can solve this problem efficiently by utilizing ridge regression 41 with
the help of a penalty term $\mu\|\Delta\boldsymbol{\phi}\|^{2}_{2}$ where
$\|\dots\|_{2}$ denotes the $\ell_{2}$-norm and $\mu$ is a hyperparameter with
values usually in the range 10-2 to 10-6, which ensures that the correction
$\Delta\boldsymbol{\phi}$ on the IFCs is as small as possible:

$\Delta\boldmath{\phi}=\arg\min_{\Delta\boldmath{\phi}}\|\mathbb{d}+\mathbb{I}\cdot\mathbb{C^{\prime}}\cdot\Delta\boldmath{\phi}\|^{2}_{2}+\mu\|\Delta\boldmath{\phi}\|^{2}_{2}.$
(24)

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 positions are within a cutoff radius.

### Fourier interpolation of interatomic force constants

The dynamical matrix $D_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{q})$ is
computed via a Fourier interpolation of the real-space IFCs as:

$D_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{q})=\sum_{\mathbf{R}}\Phi_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{0},\mathbf{R}){\rm
e}^{-{\rm i}\mathbf{q}\cdot\mathbf{R}},$ (25)

where $\mathbb{q}$ is a phonon wavevector defined on an arbitrary grid in
reciprocal space. Due to the transitional invariance, the IFCs
$\Phi_{\varkappa\alpha,\varkappa^{\prime}\beta}$ in crystalline solids depend
only on the relative position between the two atoms $\varkappa$ and
$\varkappa^{\prime}$. Therefore, it is convenient to consider only the IFCs
$\Phi_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{R})$ between the atom
labelled by $\kappa$ within the reference unit cell $\mathbf{0}$ and the atom
labelled by $\kappa^{\prime}$ within another unit cell $\mathbf{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 short-range ($\mathcal{S}$)
and long-range ($\mathcal{L}$) component as:

$D_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{q})=D_{\kappa\alpha,\kappa^{\prime}\beta}^{\mathcal{S}}(\mathbf{q})+D_{\kappa\alpha,\kappa^{\prime}\beta}^{\mathcal{L}}(\mathbf{q}),$
(26)

such that the short-range part is a smooth function of $\mathbf{q}$. If we
consider only dipole-dipole interactions 12, 45, the long-range part takes the
following form:

$D_{\kappa\alpha,\kappa^{\prime}\beta}^{\mathcal{L}}(\mathbf{q\rightarrow
0})=\frac{W_{\rm
c}(\mathbb{q})e^{2}}{\Omega_{0}}\Big{[}\big{(}\mathbb{q}\cdot\mathbb{Z}^{*}_{\kappa}\big{)}_{\alpha}\big{(}\mathbb{q}\cdot\mathbb{Z}^{*}_{\kappa^{\prime}}\big{)}_{\beta}\Big{]},$
(27)

which is non-analytic at $\mathbf{q}=\mathbf{0}$ and induces the direction-
dependent LO-TO splitting. To deal with the $\mathbf{q}=\mathbf{0}$ term
exactly, the Ewald summation technique 12, 11 is adopted in practice with a
sum over reciprocal $\mathbf{G}$ vectors. In this expression, $e$ is the
elementary charge unit, $\Omega_{0}$ is the volume of unit cell,
$\mathbb{Z}^{*}_{\kappa}$ is the Born effective charge tensor of the atom
$\kappa$ within the unit cell, and $W_{\rm c}(\mathbb{q})$ is the screened
Coulomb potential which depends on the dimensionality of the system as 44, 46,
47

$\displaystyle W_{\rm c}^{\rm
3D}(\mathbb{q})=\frac{4\pi}{\mathbb{q}\cdot\boldsymbol{\epsilon}\cdot\mathbb{q}},$
(28) $\displaystyle W_{\rm c}^{\rm
2D}(\mathbb{q})=\frac{2\pi}{|\mathbb{q}|\left(1+\frac{\mathbb{q}\cdot\mathbf{r}_{\rm
eff}\cdot\mathbb{q}}{|\mathbb{q}|^{2}}|\mathbb{q}|\right)},$ (29)
$\displaystyle W_{\rm c}^{\rm
1D}(\mathbb{q})=\frac{4\pi}{\epsilon_{\mathrm{1D}}|\mathbb{q}|^{2}\pi
t^{2}}\Bigg{\\{}1-2I_{1}(|\mathbb{q}|t)K_{1}(|\mathbb{q}|t)\bigg{[}1-$
$\displaystyle\frac{2\epsilon_{\mathrm{1D}}\sqrt{\pi}|\mathbb{q}|tI_{1}(|\mathbb{q}|t)K_{0}(|\mathbb{q}|t)-G^{22}_{24}(|\mathbb{q}|^{2}t^{2})}{2\sqrt{\pi}|\mathbb{q}|t\big{[}\epsilon_{\mathrm{1D}}I_{1}(|\mathbb{q}|t)K_{0}(|\mathbb{q}|t)+I_{0}(|\mathbb{q}|t)K_{1}(|\mathbb{q}|t)\big{]}}\bigg{]}\Bigg{\\}},$
(30)

where $\boldsymbol{\epsilon}$ is the ion-clamped dielectric tensor of 3D
system, $\mathbb{r}_{\rm eff}\approx\boldsymbol{\epsilon}b/2$ is the in-plane
matrix of effective screening length, $b$ is the thickness of the 2D material,
$\epsilon_{\rm 1D}$ is the 1D dielectric constant, $t$ is the effective
screening radius, $I_{n}(x)$ and $K_{n}(x)$ are the $n$th-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,

$\sum_{\varkappa^{\prime}}\Big{(}\Phi^{\mathcal{S}}_{\varkappa\alpha,\varkappa^{\prime}\beta}+\Phi^{\mathcal{L}}_{\varkappa\alpha,\varkappa^{\prime}\beta}\Big{)}\tau_{\varkappa^{\prime}\gamma}=\\\
\sum_{\varkappa^{\prime}}\Big{(}\Phi^{\mathcal{S}}_{\varkappa\alpha,\varkappa^{\prime}\gamma}+\Phi^{\mathcal{L}}_{\varkappa\alpha,\varkappa^{\prime}\gamma}\Big{)}\tau_{\varkappa^{\prime}\beta},$
(31)

and

$\sum_{\varkappa\varkappa^{\prime}}\Big{(}\Phi^{\mathcal{S}}_{\varkappa\alpha,\varkappa^{\prime}\beta}+\Phi^{\mathcal{L}}_{\varkappa\alpha,\varkappa^{\prime}\beta}\Big{)}\tau_{\varkappa\varkappa^{\prime}\gamma}\tau_{\varkappa\varkappa^{\prime}\delta}=\\\
\sum_{\varkappa\varkappa^{\prime}}\Big{(}\Phi^{\mathcal{S}}_{\varkappa\gamma,\varkappa^{\prime}\delta}+\Phi^{\mathcal{L}}_{\varkappa\gamma,\varkappa^{\prime}\delta}\Big{)}\tau_{\varkappa\varkappa^{\prime}\alpha}\tau_{\varkappa\varkappa^{\prime}\beta}.$
(32)

Since the $\Phi^{\mathcal{L}}$ are analytically defined, we apply the
corrections due to the invariance conditions only on $\Phi^{\mathcal{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 short-range 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 49, 45, 46, 43 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

$[\omega_{\nu}(\mathbf{q})]^{2}e_{\nu,\kappa\alpha}(\mathbf{q})=\sum_{\kappa^{\prime}\beta}\frac{D_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{q})}{\sqrt{m_{\kappa}m_{\kappa^{\prime}}}}e_{\nu,\kappa^{\prime}\beta}(\mathbf{q}),$
(33)

where $m_{\kappa}$ is the mass of the $\kappa$-th atom within the unit cell.
The eigenvector $e_{\nu,\kappa\alpha}(\mathbf{q})$ is obtained by
diagonalizing the mass-scaled dynamical matrix
$D_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{q})/\sqrt{m_{\kappa}m_{\kappa^{\prime}}}$
with the corresponding eigenfrequency $\omega_{\nu}(\mathbf{q})$ for the
phonon mode $\nu$ at momentum $\mathbf{q}$, and $\nu$ 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 ($\mathbf{q}\rightarrow 0$), we perform a
Taylor expansion of Equation (33) in terms of $\mathbf{q}$. The resulting
long-wavelength equation for lattice vibrations is 21, 10:

$\sum_{\kappa}m_{\kappa}[\omega^{(1)}_{\nu}(\mathbf{q})]^{2}u_{\nu\alpha}=\\\
\sum_{\beta\gamma\delta}\Big{[}T^{\rm
sym}_{\alpha\beta,\gamma\delta}q_{\gamma}q_{\delta}+T^{\rm
int}_{\alpha\gamma,\beta\delta}q_{\gamma}q_{\delta}\Big{]}u_{\nu\beta},$ (34)

where
$\omega^{(1)}_{\nu}(\mathbf{q})\equiv\frac{\partial\omega_{\nu}(\mathbf{q})}{\partial\mathbf{q}}\Big{|}_{\mathbf{q}=0}\cdot\mathbf{q}$.
In this equation, $u_{\nu\alpha}$ is an arbitrary polarization vector which is
introduced via the solution for the zero-order equation of the Taylor
expansion,
$e^{(0)}_{\nu,\kappa\alpha}(\mathbf{q})=\sqrt{m_{\kappa}}u_{\nu\alpha}$. The
$T^{\rm sym}_{\alpha\beta,\gamma\delta}$ defined in Equation (Equilibrium
conditions) can be rewritten in term of a sum over primitive cells as

$T^{\rm
sym}_{\alpha\beta,\gamma\delta}=-\frac{1}{2}\sum_{\mathbf{R}}\sum_{\kappa\kappa^{\prime}}\Phi_{\kappa\alpha,\kappa^{\prime}\beta}(\mathbf{R})\tau_{\kappa\kappa^{\prime}\gamma}(\mathbf{R})\tau_{\kappa\kappa^{\prime}\delta}(\mathbf{R}),$
(35)

and $T^{\rm int}_{\alpha\gamma,\beta\delta}$ is the contribution associated
with the relaxation of internal coordinates in response to external stress
fields 10:

$\displaystyle T^{\rm int}_{\alpha\gamma,\beta\delta}=$
$\displaystyle-\sum_{\kappa\kappa^{\prime}}\sum_{\lambda\mu}\frac{\Gamma_{\kappa\lambda,\kappa^{\prime}\mu}}{\sqrt{m_{\kappa}m_{\kappa^{\prime}}}}$
$\displaystyle\times\Bigg{[}\sum_{\mathbf{R}}\sum_{\kappa^{\prime\prime}}\Phi_{\kappa\lambda,\kappa^{\prime\prime}\alpha}(\mathbf{R})\tau_{\kappa\kappa^{\prime\prime}\gamma}(\mathbf{R})\Bigg{]}$
$\displaystyle\times\Bigg{[}\sum_{\mathbf{R^{\prime}}}\sum_{\kappa^{\prime\prime\prime}}\Phi_{\kappa^{\prime}\mu,\kappa^{\prime\prime\prime}\beta}(\mathbf{R^{\prime}})\tau_{\kappa^{\prime}\kappa^{\prime\prime\prime}\delta}(\mathbf{R^{\prime}})\Bigg{]},$
(36)

with the $\nu\times\nu$ matrix $\Gamma_{\kappa\lambda,\kappa^{\prime}\mu}$
denoting the inverse of the mass-scaled IFC matrix
$\sum_{\mathbf{R}}\Phi_{\kappa\lambda,\kappa^{\prime}\mu}(\mathbf{R})/\sqrt{m_{\kappa}m_{\kappa^{\prime}}}$.
Unfortunately, the IFC matrix is singular because of the acoustic sum rules
for translational invariance and the
$\Gamma_{\kappa\lambda,\kappa^{\prime}\mu}$ is thus introduced as

$\Gamma_{\kappa\lambda,\kappa^{\prime}\mu}=\begin{cases}\hat{\Gamma}_{\kappa\lambda,\kappa^{\prime}\mu}&\kappa,\kappa^{\prime}\neq
1\\\ 0&\text{otherwise}\end{cases},$ (37)

where the pseudoinverse $\hat{\Gamma}_{\kappa\lambda,\kappa^{\prime}\mu}$ is
calculated and satisfies the relation 10

$\sum_{\mu}\sum_{\kappa^{\prime}=2}\frac{\hat{\Gamma}_{\kappa\lambda,\kappa^{\prime}\mu}}{\sqrt{m_{\kappa^{\prime}}m_{\kappa^{\prime\prime}}}}\sum_{\mathbf{R}}\Phi_{\kappa^{\prime}\mu,\kappa^{\prime\prime}\beta}(\mathbf{R})=\delta_{\lambda\beta}\delta_{\kappa\kappa^{\prime\prime}}.$
(38)

Furthermore, the tensor $T^{\rm int}_{\alpha\gamma,\beta\delta}$ satisfies the
following symmetry relations 36, 10:

$T^{\rm int}_{\alpha\gamma,\beta\delta}=T^{\rm
int}_{\gamma\alpha,\beta\delta}=T^{\rm int}_{\beta\delta,\alpha\gamma},$ (39)

which result from the rotational invariance, Equation (6), and the fact that
$\Gamma_{\kappa\lambda,\kappa^{\prime}\mu}=\Gamma_{\kappa^{\prime}\mu,\kappa\lambda}$.

For the ZA bending modes of 2D materials in the long-wavelength 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:

$[\omega^{(1)}_{\rm ZA}(\mathbf{q})]^{2}=\\\
\frac{1}{\sum_{\kappa}m_{\kappa}}\sum_{\beta}\frac{u_{\beta}}{u_{z}}\sum_{\gamma,\delta\in\\{x,y\\}}\Big{[}T^{\rm
sym}_{z\beta,\gamma\delta}+T^{\rm
int}_{z\gamma,\beta\delta}\Big{]}q_{\gamma}q_{\delta},$ (40)

where the phonon band index $\nu\equiv{\rm ZA}$ is omitted for clarity. By
using the Huang conditions and the symmetry of $T^{\rm
int}_{\alpha\gamma,\beta\delta}$, the identity $T^{\rm
sym}_{z\beta,\gamma\delta}+T^{\rm int}_{z\gamma,\beta\delta}=T^{\rm
sym}_{\gamma\delta,z\beta}+T^{\rm int}_{\gamma z,\delta\beta}$ holds. As
demonstrated in Supplementary Note 1 52, we have $T^{\rm
sym}_{\gamma\delta,z\beta}+T^{\rm int}_{\gamma z,\delta\beta}=0$ even in case
of non-purely out-of-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:

$T^{\rm sym}_{\gamma\delta,z\beta}=-T^{\rm int}_{\gamma z,\delta\beta},$ (41)

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^{\rm
sym}_{\gamma\delta,x\beta}=-T^{\rm int}_{\gamma x,\delta\beta}$ and $T^{\rm
sym}_{\gamma\delta,y\beta}=-T^{\rm int}_{\gamma y,\delta\beta}$ 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

Figure 1:  Effects of the invariance conditions on the phonon dispersion of
bulk silicon. 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 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-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.

Figure 2:  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.

As an example of infrared-active bulk materials, we calculate the phonon
dispersion of triclinic CaP3 shown in Fig. 2, with the LO-TO splitting
included. The crystal symmetry of CaP3 with the space group of $\rm{P\bar{1}}$
is low: only the identity and the spatial inversion operations. For a $2\times
2\times 2$ supercell with 64 atoms, there are 160 unique atomic pairs which
means 1440 IFC components. The total number of independent IFC parameters in
CaP3 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 $\mathbf{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 CaP3 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}\neq 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

Figure 3:  Invariance conditions applied on the phonon dispersion of graphene.
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 $\Gamma$
point. The yellow line coincides with the red one as the rotational invariance
is automatically satisfied in graphene.

As revealed by a few recent studies 20, 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 long-wavelength
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 $\Gamma$ 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 MoS2 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.

Figure 4:  Comparison between the phonon dispersions of an infrared-active 2D
MoS2 with the corrections of Born-Huang (BH) rotational invariance and Huang
(H) conditions on the short-range (SR) IFCs and on the total IFCs. The polar
Born-Huang (PBH) rotational invariance and polar Huang (PH) conditions are
also considered. The TI denotes the correction of translational invariance.
When correcting the total IFCs, the contributions from the long-range (LR)
IFCs are further removed in real space and an analytic long-range correction
(LRC) is then added back in reciprocal space to recover the LO-TO splitting.
The inset show zoom-in results around the $\Gamma$ point and the LO-TO
splitting at the $\rm E^{\prime}$ mode.

Most 2D insulators and semiconductors are infrared-active 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 MoS2 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
$\mathbf{q}$ with a finite slope 44, 54. As shown in Fig. 4, the phonon
dispersions of MoS2 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 $\Gamma$ 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 $\Gamma$
point (i.e. the short-range 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

Figure 5:  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 $\Gamma$
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. Figure 6:  Application of the invariance
conditions to the phonon dispersion of an infrared-active 1D single-wall (4,4)
boron nitride 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 $\Gamma$ point, with four acoustic branches labelled as
longitudinal acoustic (LA), twisting acoustic (TWA) and flexural acoustic (ZA,
doubly degenerate) modes, respectively.

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 39, 55, 56, 25, the frequency of the
twisting acoustic mode due to the global rotation along the nanotube axis can
become finite at the $\Gamma$ 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) CNT, i.e.
the blue line coinciding with the red one. The conditions of vanishing stress
have been already satisfied after a 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 infrared-active (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 lattice-
dynamical properties.

Figure 7:  Corrected phonon dispersions of four representative 2D materials
from the 2D material database of Ref. 57, (a) FeO2, (b) PrOI, (c) InSe and (d)
BiOCl. The inset shows a zoom-in result around the $\Gamma$ point. The black
dashed lines are the results computed from the DFPT calculations with the
translational invariance (TI) fulfilled, while the red lines denote the
results where the Born-Huang (BH) invariances and Huang (H) conditions have
also been enforced.

### 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 FeO2,
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, FeO2 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.

The updated database is available on the Materials Could archive 58, and the
phonon dispersions for all of 245 materials can be found in Supplementary
Information 52 of this study. After investigation, there are 54 materials that
are not dynamically stable, with pronounced imaginary frequencies not located
around the $\Gamma$ point in the database (see Supplementary Table 1 52).
After these truly unstable 2D structures are excluded, we still have 33
candidates that exhibit an unstable ZA branch with small imaginary frequencies
near the $\Gamma$ point after the invariance conditions are enforced (see the
detailed information of these 2D materials in Supplementary Table 2 52). 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 higher-order multipolar
interactions 43, 45, 49, 46, which we neglect here for simplicity. To reveal
the first problem, we have selected a few of cases including HfS2, NbF4, CoO2,
TaS2 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 ﬁeld in
calculating IFCs and higher-order long-range electrostatic fields in the
interpolation of IFCs, which could be important in the presence of
piezoelectricity or flexoelectricity 59, 60, 45, 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 $\Gamma$ 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 40, 15 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 plane-wave kinetic energy cutoff of 60 Ry with a
$\Gamma$-point grid computed on a $4\times 4\times 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 CaP3

The structure is relaxed with a plane-wave kinetic energy cutoff of 30 Ry and
a $10\times 10\times 10$ electronic grid. Similar to the case of silicon, we
build a $2\times 2\times 2$ supercell, and the real-space DFT calculations
based on a $2\times 2\times 2$ grid are carried out for 30 displaced
structures to get the entire IFCs. In addition, a DFPT calculation at
$\mathbb{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\times 20\times 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\times 6\times 1$ supercell and the
corresponding DFT calculations for 10 displaced supercells, with an electronic
sampling of $3\times 3\times 3$.

### 2D MoS2

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\times 14\times 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\times 6\times 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\times 1\times 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 $\Gamma$-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 Quantum ESPRESSO distribution, which will be
released in the next version of Quantum ESPRESSO.

## ACKNOWLEDGEMENTS

We would like to thank Norma Rivano and Thibault Sohier for sharing with us
the results of Ref. 47 prior to publication, and for the calculations of the
interatomic 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. 21. We also thank Junfeng Qiao, Elsa Passaro, Giovanni
Pizzi, Yanbo Li, Jian Han, Xavier Gonze, Jiongzhi Zheng, Raffaele Resta,
Massimiliano Stengel and Miquel Royo for useful discussions. This work is
supported by the Sinergia project of the Swiss National Science Foundation
(No. CRSII5_189924). S.P. acknowledges financial support from the Belgian
F.R.S.-FNRS.

## AUTHOR CONTRIBUTIONS

C.L. implemented the invariance conditions and performed the calculations.
S.P. and N.M. supervised the project. All authors analyzed the results and
contributed to writing the manuscript.

## COMPETING INTERESTS

The authors declare no competing interests.

## References

  * 1 Dove, M. T.  _Introduction to lattice dynamics_ (Cambridge university press, 1993).
  * 2 Grimvall, G.  _The electron-phonon interaction in metals_ (North-Holland Amsterdam, 1981).
  * 3 Giustino, F.  Electron-phonon interactions from first principles.  _Rev. Mod. Phys._ 89, 015003 (2017).
  * 4 Poncé, S., Li, W., Reichardt, S. & Giustino, F.  First-principles calculations of charge carrier mobility and conductivity in bulk semiconductors and two-dimensional materials.  _Rep. Prog. Phys._ 83, 036501 (2020).
  * 5 Shindé, S. L. & Srivastava, G. P.  _Length-scale dependent phonon interactions_ (Springer, 2014).
  * 6 Qian, X., Zhou, J. & Chen, G.  Phonon-engineered extreme thermal conductivity materials.  _Nat. Mater._ 20, 1188–1202 (2021).
  * 7 Einstein, A.  Die plancksche theorie der strahlung und die theorie der spezifischen wärme.  _Ann. Phys._ 327, 180–190 (1907).
  * 8 Born, M. & Von Karman, T.  Uber schwingungen im raumgittern.  _Phys. Z._ 13, 297–309 (1912).
  * 9 Debye, P.  Zur theorie der spezifischen wärmen.  _Ann. Phys._ 344, 789–839 (1912).
  * 10 Born, M. & Huang, K.  _Dynamical Theory of Crystal Lattices_ (Oxford University Press, 1954).
  * 11 Baroni, S., de Gironcoli, S., Corso, A. D. & Giannozzi, P.  Phonons and related crystal properties from density-functional perturbation theory.  _Rev. Mod. Phys._ 73, 515–562 (2001).
  * 12 Gonze, X. & Lee, C.  Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory.  _Phys. Rev. B_ 55, 10355–10368 (1997).
  * 13 Alfè, D.  PHON: A program to calculate phonons using the small displacement method.  _Comput. Phys. Commun._ 180, 2622–2633 (2009).
  * 14 Togo, A. & Tanaka, I.  First principles phonon calculations in materials science.  _Scr. Mater._ 108, 1–5 (2015).
  * 15 Giannozzi, P. _et al._ Advanced capabilities for materials modelling with quantum ESPRESSO.  _J. Phys.: Condens. Matter_ 29, 465901 (2017).
  * 16 Gonze, X. _et al._ Recent developments in the ABINIT software package.  _Comput. Phys. Commun._ 205, 106–131 (2016).
  * 17 Wang, Y., Shang, S.-L., Fang, H., Liu, Z.-K. & Chen, L.-Q.  First-principles calculations of lattice dynamics and thermal properties of polar solids.  _npj Comput. Mater._ 2, 1–10 (2016).
  * 18 Petretto, G. _et al._ High-throughput density-functional perturbation theory phonons for inorganic materials.  _Sci. Data_ 5, 1–12 (2018).
  * 19 Eriksson, F., Fransson, E. & Erhart, P.  The hiphive package for the extraction of high-order force constants by machine learning.  _Adv. Theor. Simul._ 2, 1800184 (2019).
  * 20 Carrete, J. _et al._ Physically founded phonon dispersions of few-layer materials and the case of borophene.  _Mater. Res. Lett._ 4, 204–211 (2016).
  * 21 Croy, A.  Bending rigidities and universality of flexural modes in 2D crystals.  _J. Phys.: Mater._ 3, 02LT03 (2020).
  * 22 Dresselhaus, M. S. & Eklund, P. C.  Phonons in carbon nanotubes.  _Adv. Phys._ 49, 705–814 (2000).
  * 23 Mahan, G. D. & Jeon, G. S.  Flexure modes in carbon nanotubes.  _Phys. Rev. B_ 70, 075405 (2004).
  * 24 Mingo, N., Stewart, D. A., Broido, D. A. & Srivastava, D.  Phonon transmission through defects in carbon nanotubes from first principles.  _Phys. Rev. B_ 77, 033418 (2008).
  * 25 Libbi, F., Bonini, N. & Marzari, N.  Thermomechanical properties of honeycomb lattices from internal-coordinates potentials: the case of graphene and hexagonal boron nitride.  _2D Mater._ 8, 015026 (2020).
  * 26 Kosevich, A. M.  _The crystal lattice: phonons, solitons, dislocations, superlattices_ (John Wiley & Sons, 2006).
  * 27 Leibfried, G. & Ludwig, W.  Theory of anharmonic effects in crystals.  In _Solid State Physics_ , 275–444 (Elsevier, 1961).
  * 28 Noether, E.  Invariant variation problems.  _Transp. Theory Stat. Phys._ 1, 186–207 (1971).
  * 29 Begbie, G. H. & Born, M.  Thermal scattering of x-rays by crystals - I. dynamical foundation.  _Proc. R. Soc. A_ 188, 179–188 (1947).
  * 30 Keating, P. N.  Relationship between the macroscopic and microscopic theory of crystal elasticity. I. primitive crystals.  _Phys. Rev._ 152, 774–779 (1966).
  * 31 Keating, P. N.  Effect of invariance requirements on the elastic strain energy of crystals with application to the diamond structure.  _Phys. Rev._ 145, 637–645 (1966).
  * 32 Gazis, D. C. & Wallis, R. F.  Conditions for rotational invariance of a harmonic lattice.  _Phys. Rev._ 151, 578–580 (1966).
  * 33 Sarkar, S. K. & Sengupta, S.  On born-huang invariance conditions.  _Phys. Status Solidi B_ 83, 263–271 (1977).
  * 34 Pick, R. M., Cohen, M. H. & Martin, R. M.  Microscopic theory of force constants in the adiabatic approximation.  _Phys. Rev. B_ 1, 910–920 (1970).
  * 35 Huang, K.  LXX. lattice theory of dielectric and piezoelectric constants in crystals.  _Philos. Mag._ 40, 733–747 (1949).
  * 36 Huang, K.  On the atomic theory of elasticity.  _Proc. R. Soc. A_ 203, 178–194 (1950).
  * 37 Mounet, N. _et al._ Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds.  _Nat. Nanotechnol._ 13, 246–252 (2018).
  * 38 Zhou, F., Nielson, W., Xia, Y. & Ozoliņš, V.  Compressive sensing lattice dynamics. I. general formalism.  _Phys. Rev. B_ 100, 184308 (2019).
  * 39 Mounet, N.  _Structural, vibrational and thermodynamic properties of carbon allotropes from first-principles: diamond, graphite, and nanotubes_.  Master’s thesis, Massachusetts Institute of Technology (2005).
  * 40 Giannozzi, P. _et al._ QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials.  _J. Phys.: Condens. Matter_ 21, 395502 (2009).
  * 41 Hilt, D. E. & Seegrist, D. W.  _Ridge, a computer program for calculating ridge regression estimates_ , vol. 236 (Department of Agriculture, Forest Service, Northeastern Forest Experiment, 1977).
  * 42 Poncé, S., Margine, E., Verdi, C. & Giustino, F.  EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized wannier functions.  _Comput. Phys. Commun._ 209, 116–133 (2016).
  * 43 Poncé, S. _et al._ First-principles predictions of hall and drift mobilities in semiconductors.  _Phys. Rev. Res._ 3, 043022 (2021).
  * 44 Sohier, T., Gibertini, M., Calandra, M., Mauri, F. & Marzari, N.  Breakdown of optical phonons’ splitting in two-dimensional materials.  _Nano Lett._ 17, 3758–3763 (2017).
  * 45 Royo, M., Hahn, K. R. & Stengel, M.  Using high multipolar orders to reconstruct the sound velocity in piezoelectrics from lattice dynamics.  _Phys. Rev. Lett._ 125, 217602 (2020).
  * 46 Royo, M. & Stengel, M.  Exact long-range dielectric screening and interatomic force constants in quasi-two-dimensional crystals.  _Physical Review X_ 11, 041027 (2021).
  * 47 Rivano, N., Marzari, N. & Sohier, T.  Infrared-active phonons in one-dimensional materials and their spectroscopic signatures. Preprint at {https://arxiv.org/abs/2208.09887} (2022).
  * 48 Meijer, C. S.  Über Whittakersche bzw. Besselsche Funktionen und deren Produkte.  _Nieuw Arch. Wiskd. II. Ser._ 18, 10–39 (1936).
  * 49 Brunin, G. _et al._ Electron-phonon beyond fröhlich: Dynamical quadrupoles in polar and covalent solids.  _Phys. Rev. Lett._ 125, 136601 (2020).
  * 50 Zhou, F., Sadigh, B., Åberg, D., Xia, Y. & Ozoliņš, V.  Compressive sensing lattice dynamics. II. efficient phonon calculations and long-range interactions.  _Phys. Rev. B_ 100, 184309 (2019).
  * 51 Gonze, X., Charlier, J.-C., Allan, D. & Teter, M.  Interatomic force constants from first principles: The case of $\alpha$-quartz.  _Phys. Rev. B_ 50, 13035–13038 (1994).
  * 52 See Supplementary Information for additional ﬁgures and tables, the derivation of bending conditions in 2D materials, and further discussion.
  * 53 Sluiter, M. H. F., Weinert, M. & Kawazoe, Y.  Determination of the elastic tensor in low-symmetry structures.  _Europhys. Lett._ 43, 183–188 (1998).
  * 54 Sohier, T., Calandra, M. & Mauri, F.  Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene.  _Phys. Rev. B_ 96, 075448 (2017).
  * 55 Mounet, N. & Marzari, N.  Phonon dispersions of a single-wall (8, 0) carbon nanotube: Effects of the rotational acoustic sum rule and of surface attachment.  _MRS Online Proceedings Library_ 858, 220–225 (2004).
  * 56 Ye, L.-H., Liu, B.-G., Wang, D.-S. & Han, R.  Ab initio phonon dispersions of single-wall carbon nanotubes.  _Phys. Rev. B_ 69, 235409 (2004).
  * 57 Mounet, N. _et al._ Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds.  _Materials Cloud Archive_ 2020.158, https://doi.org/10.24435/materialscloud:az-b2 (2020).
  * 58 Lin, C., Poncé, S. & Marzari, N.  General invariance and equilibrium conditions for lattice dynamics in 1D, 2D, and 3D materials.  _Materials Cloud Archive_ 2022.111, https://doi.org/10.24435/materialscloud:gf-3n (2022).
  * 59 Wu, X., Vanderbilt, D. & Hamann, D. R.  Systematic treatment of displacements, strains, and electric fields in density-functional perturbation theory.  _Phys. Rev. B_ 72, 035105 (2005).
  * 60 Ma, W. _et al._ Piezoelectricity in multilayer black phosphorus for piezotronics and nanogenerators.  _Advanced Materials_ 32, 1905795 (2020).
  * 61 Royo, M. & Stengel, M.  Lattice-mediated bulk flexoelectricity from first principles.  _Phys. Rev. B_ 105, 064101 (2022).
  * 62 Prandini, G., Marrazzo, A., Castelli, I. E., Mounet, N. & Marzari, N.  Precision and efficiency in solid-state pseudopotential calculations.  _npj Comput. Mater._ 4, 1–13 (2018).
  * 63 Perdew, J. P. _et al._ Restoring the density-gradient expansion for exchange in solids and surfaces.  _Phys. Rev. Lett._ 100, 136406 (2008).

in 1,…,0 See pages \x, of si