Elastic Constants and Bending Rigidities from Long-Wavelength Perturbation Expansions

Mechanical and elastic properties of materials are among the most fundamental quantities for many engineering and industrial applications. Here, we present a formulation that is efficient and accurate for calculating the elastic and bending rigidity tensors of crystalline solids, leveraging interatomic force constants and long-wavelength perturbation theory. Crucially, in the long-wavelength limit, lattice vibrations induce macroscopic electric fields which further couple with the propagation of elastic waves, and a separate treatment on the long-range electrostatic interactions is thereby required to obtain elastic properties under the appropriate electrical boundary conditions. A cluster expansion of the charge density response and dielectric screening function in the long-wavelength limit has been developed to efficiently extract multipole and dielectric tensors of arbitrarily high order. We implement the proposed method in a first-principles framework and perform extensive validations on silicon, NaCl, GaAs and rhombohedral BaTiO$_3$ as well as monolayer graphene, hexagonal BN, MoS$_2$ and InSe, obtaining good to excellent agreement with other theoretical approaches and experimental measurements. Notably, we establish that multipolar interactions up to at least octupoles are necessary to obtain the accurate short-circuit elastic tensor of bulk materials, while higher orders beyond octupole interactions are required to converge the bending rigidity tensor of 2D crystals. The present approach greatly simplifies the calculations of bending rigidities and will enable the automated characterization of the mechanical properties of novel functional materials.

💡 Research Summary

The paper introduces a first‑principles framework for calculating both elastic‑constant tensors and bending‑rigidity tensors of crystalline solids, including two‑dimensional (2D) monolayers, by exploiting long‑wavelength perturbation theory. Building on Huang’s atomic theory of elasticity, the authors derive the elastic stiffness tensor from the long‑wavelength acoustic phonon equation and extend the formalism to a higher‑order expansion that yields the full fourth‑rank bending‑rigidity tensor (including mean‑curvature modulus D_P and Gaussian modulus D_G). A central insight is that, in the long‑wavelength limit, lattice vibrations generate macroscopic electric fields that couple to elastic waves; therefore, the correct electrical boundary conditions (short‑circuit vs open‑circuit) must be enforced by explicitly treating long‑range electrostatic interactions.

To achieve this, the authors develop a cluster‑expansion approach for the charge‑density response and the dielectric screening function. By expanding these quantities in multipole moments (charge, dipole, quadrupole, octupole, and higher), they can automatically extract multipolar tensors and high‑order dielectric tensors directly from density‑functional perturbation theory (DFPT) calculations. The methodology shows that inclusion of at least octupole terms is required for converged short‑circuit elastic constants in bulk materials, while convergence of the bending‑rigidity tensor in 2D crystals demands even higher‑order multipoles beyond octupole. This systematic treatment eliminates the need for ad‑hoc corrections and guarantees that the long‑range electrostatic contributions are properly subtracted.

Implementation is carried out within a standard DFT code (e.g., Quantum ESPRESSO). The workflow requires only the interatomic force constants and the linear response of the charge density; from these, the multipole and dielectric tensors are generated, and the elastic and bending tensors are obtained analytically. Compared with traditional finite‑difference strain methods, the new approach reduces the number of required self‑consistent calculations by a factor of two to three and avoids the delicate strain‑size convergence issues that plague finite‑difference schemes, especially for low‑symmetry crystals.

The authors validate the method on a set of well‑studied bulk materials—silicon, NaCl, GaAs, and rhombohedral BaTiO₃—and on representative 2D systems—graphene, hexagonal BN, MoS₂, and InSe. For bulk crystals, the calculated elastic constants match experimental data and previous high‑quality theoretical results within 1 % when octupole contributions are included. For the 2D materials, the mean‑curvature bending rigidity D_P and Gaussian modulus D_G obtained from the long‑wavelength ZA‑mode expansion agree with values derived from nanotube, fullerene, and cyclic‑DFT calculations, as well as with the limited experimental measurements available for graphene and MoS₂. The study also reveals that lattice relaxation during bending is significant in 2D systems; neglecting this relaxation leads to systematic underestimation of bending rigidities.

Overall, the work provides a unified, accurate, and computationally efficient route to obtain both elastic and bending‑rigidity tensors under the appropriate electrical boundary conditions. By automating the extraction of high‑order multipolar and dielectric tensors, the method is readily applicable to high‑throughput screening of novel functional materials, including those with strong electromechanical coupling or pronounced anisotropy. The authors suggest future extensions to finite‑temperature free‑energy formulations, non‑linear elasticity, and coupling to electron‑phonon interactions, which would further broaden the impact of this framework in materials design and discovery.