Bending rigidities and universality of flexural modes in 2D crystals

The existence of flexural modes with a quadratic phonon-dispersion is a distinguishing property of two-dimensional materials and has important consequences for their properties. Here, we deduce theoretically within the harmonic approximation the conditions for which orthotropic two-dimensional materials display a flexural mode. Further, we derive formulae for the calculation of the corresponding bending rigidities using the equilibrium structure and the second-order force constants as input. This completes the description of the elasticity of 2D crystals. Our findings are exemplarily validated by ab initio calculations of the phonon dispersions of four representative materials.

Since the advent of graphene [1], two-dimensional materials (2DMs) have attracted a substantial amount of interest and triggered a wide range of discoveries [2][3][4][5][6]. Among the many outstanding and distinguishing properties of 2DMs, the existence of flexural modes (also called bending or ZA modes), which display a quadratic phonon dispersion for small wave-vectors, is of fundamental importance [7][8][9][10]. After all, the mechanical properties of 2DMs depend on the interplay of in-plane and out-of-plane phonons (in the form of ripples) [11][12][13] -eventually evading the consequences of the Mermin-Wagner theorem [14].
Up to now hundreds of ab initio calculations have been performed to characterize existing and potential 2DMs [15][16][17][18]. However, the nature of the flexural modes is hard to pinpoint computationally, which led to some controversy about the universality of observing flexural modes for 2DMs (see reference [19] for a list of publications reporting linear or quadratic ZA dispersions for the same 2DM) despite their indisputable existence in the continuum limit. Moreover, in contrast to 3D crystals, the knowledge of the elastic constants is not sufficient to obtain a description of the long-wavelength behavior of 2DMs. In addition to the elastic constants, the bending rigidities are required, which also determine the dispersion of the flexural modes [8,11,13,20]. While general expressions for the elastic constants of crystals are long known [21], there are no such general expressions for the bending rigidities. In this Letter, we derive formulae for the calculation of the curvatures of ZA dispersions, i.e. the bending rigidities, of 2DMs using their equilibrium structure and the second-order force constants (FCs) as an input. Further, we theoretically show that all (freestanding) orthotropic 2DMs exhibit a flexural mode provided that the material is free of homogeneous stresses. The derivations emphasize the role of sum rules for the FCs which reflect the underlying fundamental symmetries of the crystals [19,22]. Our findings are corroborated by considering four representative 2DMs for which we calculate the phonon dispersions from first principles using VASP [23] and PHONOPY [24].
Taking the normal vector of the 2DM to be parallel to the z-axis, the primitive vectors ⃗ a 1 and ⃗ a 2 lie in the x − y plane with vanishing z-component (cf figure 1). The equilibrium positions of the k = 1, …, N atoms in the primitive cell ℓ are denoted by ⃗ x  with the ('C-type') dynamical matrix [25] C αβ Here, ⃗ x is the distance vector between the atoms (ℓ ′ , k) and (0, k ′ ), m k is the mass of atom k and Φ are the FCs. Greek indices will always indicate cartesian coordinates (x, y, z) in the following. The j = 1, …, 3 N eigenvalues of equation (1) are denoted by ω 2 (j⃗ q) and the respective The FCs obey a set of sum rules, which result from the invariance of the total energy against translations and rotations [21]. The exact form of the sum rules for the latter category has been subject of debate in the past [26,27], but for our purposes, we only need the following three relations The first two sum rules are related to translational invariance and the last one is due to rotational invariance.
To extract the behavior for long wave-lengths (small ⃗ q) perturbation theory can be used [21]. Accordingly, we set ⃗ q → ϵ⃗ q and the frequencies, the eigenfunctions and the dynamical matrix in the secular equation (1) are expanded in terms of ε. Specifically, for the frequencies one has Keeping terms up to ε 2 in the expansion one finds an equation for the long-wavelength acoustic waves [21] ( ∑ The square and round brackets denote coefficients which are given in terms of the FCs and the equilibrium atomic positions, and δ αγ represents the usual Kronecker delta. The matrix Γ thus acts as the inverse of the FCs. The brackets fulfill the symmetry relations which follow from translational and rotational invariance given by equations (3). Details of the derivation can be found in reference [21] and in the SI. Here we only note, that for Bravais lattices with one atom in the primitive cell the round brackets vanish since Γ = 0. Otherwise, the renormalization in equation (5) due to the round brackets stems from the rearrangement of the atoms in the unit cell in response to external strain. This is the same mechanism leading to the violation of the Cauchy-Born rule for 2DMs [28,29], which is also important for a quantitative description of the strain-engineering of the band-gap of 2D semiconductors [30]. It is crucial for the emergence of flexural modes as we will show below.
More generally, the softening of the elastic constants represented by the round brackets can be related to 'nonaffine displacements' [31][32][33], which have been shown to be important for the quantitative description of non-centrosymmetric crystals and amorphous solids [33,34]. For a detailed discussion of the relation of nonaffine displacements to the long-waves theory of Born and Huang, the reader is referred to reference [35].
From the theory of elasticity, one expects the brackets to be related to the tensor of elastic constants if the material initially is free of stresses [27]. Under this premise the (Huang) compatibility conditions have to be fulfilled [21]. Based on those relations it can be shown that the elements of the tensor of elastic constants C are given by where A cell is the area of the unit cell [21]. Up to this point, all the equations also hold for 3D materials.
Focusing on 2DMs, we only consider wave-vectors ⃗ q = (q x , q y , 0) in the following.
In order to observe a flexural mode, the frequency contribution ω (1) for that particular mode has to vanish. In this case, the next order ω (2) will determine the dispersion. For orthotropic 2DMs, the polarization of the flexural mode is purely in z-direction 1 . In this case, one sees from equation (5) where only the symmetry relations (cf equation (8b)) of the round brackets have been used (see also SI). As we show explicitly in the SI, the sum of the first two terms vanishes for all orthotropic 2DMs. This set of bending-mode conditions reads In particular, they are fulfilled when both sides are independently zero. This is the case if the material is atomically flat, like graphene or hBN. One finds [γλ, zz] = 0 and (γz, λz) = 0, which are identically zero since = 0 for all atoms 2 . For all other cases, the bending conditions imply that the relaxation of the internal degrees of freedom partially compensates the external strain. If in addition to the bending-mode conditions the compatibility conditions given by equations (9) hold, the last two terms in the sum in equations (11) cancel and ω (1) ZA vanishes. In this case, a flexural mode emerges and the next orders in the ε-expansion have to be considered, which is done below. One can show that the difference between the square brackets, [zz, γλ] − [γλ, zz], is proportional to the anisotropic stress components [21,27]. Consequently, if the material initially is subject to homogeneous stress, the flexural modes are suppressed and the ZA-dispersion will become linear in q to leading order. This shows that a strict geometry optimization in ab initio calculations is crucial to obtain quadratic dispersions. It further corroborates the observation that the translational and rotational invariance of the FCs are necessary but not sufficient conditions for finding flexural modes [19,22].
To get an expression for the bending frequency ω (2) ZA , the perturbation expansion leading to equation (5) has to be continued up to ε 4 . The derivation follows the one leading to equation (5) and is presented in the SI. The resulting ZA dispersion is given in terms of generalized brackets, The new round brackets contain the renormalization effect due to the internal degrees of freedom (see SI), while the square brackets are in general non-vanishing for Bravais lattices.
Analogous to the relation of the brackets with the tensor of elastic constants in equation (10), the sixteen coefficients multiplying the components of the wave vectors in equation (13) are related to the bending rigidities D γλξζ . The symmetric rank-four tensor D containing the bending rigidities has in general nine distinct elements [36]. For orthotropic materials, this number reduces further to four elements. In terms of those elements, the dispersion of the flexural mode becomes  (5) and (13). Lower row: tensor of elastic constants obtained from equation (10). Values smaller than 10 −3 N m −1 are suppressed. Table 1. Elastic properties obtained from the tensors of elastic constants and bending-rigidities: Lamé parameters λ and µ, Young modulus Y and Poisson ratio ν, and principal bending rigidity D.
where ρ is the 2D mass density. For (D xxyy + D yyxx + 4D xyxy ) 2 ≤ 4D xxxx D yyyy , the right hand-side of the equation is positive and the crystal is stable. If the x − and y − axes coincide with the principal bending axes, the two sides of the positivity condition are equal and D xxxx and D yyyy correspond to the principal bending rigidities. In the case of an isotropic material, there is only one principal bending-rigidity D = D xxxx = D yyyy = (D xxyy + D yyxx + 4D xyxy )/2. The so-called Gaussian rigidity is given by D G = 4D xyxy . The bending rigidities and the elastic constants together allow the formulation of a continuum model for long-wavelength deformations of isotropic [8,11,13] and anisotropic 2DMs [20].
To corroborate our findings we have performed DFT calculations as as implemented in VASP [23] for graphene, silicene, 2 H −MoS 2 and InSe. The generalized gradient approximation (GGA) according to Perdew, Burke, and Ernzerhof (PBE) [37] was used with projector-augmented wave-potentials [38]. The phonon dispersions are obtained using the finite-displacement method of PHONOPY [24] for supercells of sizes up to 10 × 10. Further details of the calculations are given in the SI. The resulting low-frequency branches of the phonon dispersions are shown in the upper part of figure 2 (full lines) together with the dispersions obtained from equations (5) and (13) (dashed lines) for the LA/TA and ZA modes, respectively. The quantities entering those relations, the FCs and the equilibrium geometry, are obtained from the DFT calculations using PHONOPY. Overall, the agreement of the analytic results and the numerical values is quite good. One observes a systematic overestimation of the curvature which can be traced back to non-vanishing elements of the tensor of elastic constants involving the out-of-plane direction, see the lower part of figure 2. Those are on the order of 1 N m −1 for the largest unit-cells and thus much smaller than the other elastic constants. The calculated elastic properties are summarized in table 1 and are compatible with values reported in the literature [6].
In conclusion, we have shown within the harmonic approximation that orthotropic 2DMs universally display a flexural mode with a quadratic phonon dispersion if they are initially free of stress. We have further obtained expressions for the bending rigidities in terms of the FCs and the equilibrium structure. Those formulae allow for a complete characterization of the elastic properties of 2DMs and are instrumental for high-throughput calculations.