3D continuum phonon model for group-IV 2D materials

A general three-dimensional continuum model of phonons in two-dimensional materials is developed. Our first-principles derivation includes full consideration of the lattice anisotropy and flexural modes perpendicular to the layers and can thus be applied to any two-dimensional material. In this paper, we use the model to not only compare the phonon spectra among the group-IV materials but also to study whether these phonons differ from those of a compound material such as molybdenum disulfide. The origin of quadratic modes is clarified. Mode coupling for both graphene and silicene is obtained, contrary to previous works. Our model allows us to predict the existence of confined optical phonon modes for the group-IV materials but not for molybdenum disulfide. A comparison of the long-wavelength modes to density-functional results is included.


Introduction
Phonon spectra in two-dimensional (2D) nanomaterials have almost exclusively been computed using density-functional theory (DFT) based codes. One of the earliest applications to group-IV elemental 2D materials was for the important prediction of the stability of silicene and germanene [1]. These are complex calculations and prone to qualitative errors due to the various approximations such as convergence criteria and use of approximate functionals [2]. Even the stability or not of a given structure could be incorrectly inferred. For example, borophene and indiene have been predicted to be unstable in one paper [3], though other calculations (and in the case of borophene, even experiments [4]) have obtained opposite results [5,6]. In addition to obtaining a spectrum, it is often also useful to be able to predict and/or interpret properties of the spectrum based upon either microscopic or symmetric arguments. An excellent example is the prediction of a Dirac cone for silicene on the basis of symmetry [7] when DFT calculations either failed to recognize it [8] or were unable to explain it [1].
An alternative model of lattice vibrations is a classical continuum model, which is expected to reproduce most accurately phonons with wavelengths longer than lattice separations, i.e., near k = 0. One of the earliest such models applied to an ionic crystal slab is that of Fuchs and Kliewer [9], from which they deduced that transverse optical (TO) and longitudinal optical (LO) modes have different frequencies at k ≈ 0 as well as the existence of surface optical (SO) modes with an exponential dependence away from the slab. Slightly different models have been introduced by a number of authors for graphene [10][11][12][13]. One commonality is to treat the sheet as strictly two-dimensional. Additionally, instead of deriving the phonon dispersion relations from first principles, they all assumed the known results that there are in-plane acoustic modes with linear dispersions and out-of-plane transverse acoustic modes with quadratic dispersions (the latter being consistent with the elastic theory of thin plates) to construct either a Lagrangian [10] or equations of motion [11]. Goupalov also considered optical phonons but simply parameterized the dispersion relations to match the experimental data. In all of the above, the out-ofplane vibrations were assumed decoupled from the in-plane ones.
In this paper, a continuum theory of acoustic and optical phonons in 2D nanomaterials is derived from first principles, contrary to earlier approaches, by starting with the elastic and electric equations, and taking into account the full crystalline symmetry and piezoelectric couplings when allowed by symmetry. We apply the theory to obtain the phonons in group-IV elemental 2D materials. Given that there are two fundamental structures for the free-standing sheets (the flat hexagonal structure of graphene and the buckled hexagonal structure of the other elements -we will refer to silicene as the prototypical example of the latter), we will consider both of them. It should be recognized that, to date, silicene has been grown on substrates in different reconstruction state. The reconstruction is an atomic-scale distinction that is not describable by the current continuum model. Substrate effects on the phonons can be considered in an extended model that would need to be solved numerically. This can be studied in a future publication as it is beyond the analytical solutions sought for in the current manuscript. Furthermore, we have included a study of a well-known compound 2D material (molybdenum disulfide MoS 2 ) in order to further understand the properties derived for the elemental materials.

Continuum model
The 2D materials will be treated as 3D thin-plate materials in the following, well aware that the out-of-plane dimension contains one or a few atomic layers. In general, this would allow one to study multilayers though we will only focus on monolayers in this paper. Nonetheless, it is important to keep the third coordinate in the analysis to reveal the true phonon dispersions as observed experimentally and in DFT calculations.
The general 3D elastic equations are given by the equation of motion [14] (1) where T ik is the stress tensor, ρ is the mass density, and u i is the displacement. Equation 1 contains all the physics of the problem. In order to simplify and then solve it, it is necessary to specify the crystal symmetry of the vibrating system.
The three problems considered in this article are graphene, silicene and MoS 2 . The Bravais lattice symmetries of the singlelayer graphene (D 6h ≡ 6/mmm) and MoS 2 ( ) structures belong to the hexagonal system, while silicene belongs to the trigonal system (point group D 3d ). Graphene and silicene are non-piezoelectric materials because of the inversion symmetry of the unit cell, while MoS 2 is piezoelectric because its unit cell exhibits inversion asymmetry.

Application: graphene
For both graphene and MoS 2 , the general form of the stiffness tensor for hexagonal structures is (2) and the stress-strain relations T I = c IJ S J for graphene become (3) where T I and S J denote stress and strain, respectively. Here we have used Voigt notation for tensors. The latter equations are different for MoS 2 due to the presence of piezoelecticity.

Elastic equations
Combining Equation 1 and Equation 3, we obtain for graphene (4) where ω is the vibrational (angular) frequency. The displacements u x and u y are in-plane displacements, and u z is the out-ofplane displacement.
Expressing the latter set of equations in the displacements, we get (7) (8) (9) A solution of the combined system (Equation 7-Equation 9) with appropriate boundary conditions allows for the determination of the displacements u x , u y , u z . Note that the displacements in the three directions are coupled. This is different from the 2D model assumed before, which led to a decoupling of the out-ofplane vibrations. The coupling is a consequence of the finite thickness of the sheet with no mirror symmetry imposed. Thus, our model is sufficiently general to apply to multilayers. Earlier DFT calculations [15] had argued that there is no coupling between the in-plane and out-of-plane modes due to the mirror symmetry of graphene, leading to fewer scattering channels and, therefore, a higher thermal conductivity compared to, e.g., silicene. Our model reveals that such mode coupling, even when present, would occur for large k z values (due to the small thicknesses) and, hence, would be unlikely to have a significant impact.

Acoustic phonons
The phonons are normal mode solutions to the equations of motion. To proceed from the graphene elastic equations, we make the following plane wave ansatz (10) (11) (12) where k x and k y are wave numbers. Inserting Equation 10-Equation 12 into Equation 7-Equation 9 yields the following matrix expression in the unknown functions f x , f y , and f z : (13) where , and A semi-analytic solution can be easily obtained for the case when k y = 0. In this case, C = D 2 = 0 and the u y displacement decouples from the displacements u x and u z . It follows that f y obeys the wave equation (15) The solution to this differential equation can be found immediately by imposing the vacuum boundary condition (16) i.e., where 2h is the graphene layer thickness and −h,h define the graphene layer boundaries.
A 4 × 4 matrix equation in β i is finally obtained by invoking the mechanical boundary conditions (graphene in vacuum) Observe that (26) Solving the determinantal equation for the latter 4 × 4 matrix equation as a function of k x specifies a discrete set of (band) eigenfrequencies ω i (k x ) and the corresponding eigenmodes f x , f z where i denotes the band index. A numerical solution reveals the out-of-plane mode to be quadratic in nature.

Surface optical phonons
Surface optical phonons can be derived by solving for the electrostatic potential via the Maxwell-Poisson equation for the displacement field D i .
For graphene, Using the fact that the dielectric function is given by (28) we find We note in passing that the latter equation agrees with the result of Licari and Evrard [16] for interface optical phonon modes in cubic crystals where symmetry forces ε xx = ε zz ≡ ε. In the cubic case, however, confined optical phonon modes also exist at the LO phonon frequency since ε(ω LO ) = 0.

Application: silicene
We refer to silicene as the canonical example of the other group-IV materials even though the following derivation and qualitative properties apply to all of them since they all have the same symmetry of the buckled hexagonal structure, leading to the trigonal D 3d point group.
The general form of the stiffness tensor for D 3d trigonal structures is Again, the differences compared to graphene are the c 14 terms. The electric equation for trigonal structures is exactly the same as for hexagonal structures since the permittivity matrix is the same. Hence, the general method for finding surface optical phonon modes repeats the description of graphene. Qualitatively, there is therefore no difference among the phonons for all the group-IV elemental 2D materials.

Application: MoS 2
While the primary focus of this paper is on the properties of phonons of the group-IV elemental materials, it is important to know if there are properties of the phonons that are due to these being elemental. Hence, we will now consider the phonons for MoS 2 as a prototypical compound 2D material because of its extensively studied properties. As already mentioned, the lack of inversion symmetry leads to the new phenomenon of piezoelectricity.

Elastic and electric equations for MoS 2
The stiffness tensor for MoS 2 is the same as for graphene, Equation 2, because of the same hexagonal symmetry. However, the stress-strain relations are different because of the presence of piezoelectricity. Specifically, there are additional contributions to the stress-strain constitutive relations:

Phonon modes
In the case of MoS 2 , the general solutions are more complicated. One can still make the plane wave ansatz, where ε 0 is the vacuum permittivity. Solving the determinantal equation for the 8 × 8 matrix equation as a function of k x specifies a discrete set of (band) eigenfrequencies ω i (k x ) and the corresponding eigenmodes f x , f z and , where i denotes the band index. For each term m, this latter condition, in general, would lead to a different ω m . However, a normal mode is characterized by a unique frequency ω. Hence, for confined optical phonon modes, only one term in m is allowed in the general Fourier series expansions above. Further, imposing continuity in the normal electric displacement component at z = ±h by use of the third equation in Equation 55 gives (87) and

Confined phonon modes
unless accidentally ε zz (ω m ) = 0 (treated separately in the next paragraph). For ε zz (ω m ) ≠ 0, it follows from Equation 80 that A x,m = 0 and Equation 81 yields A z,m = 0. The conclusion is that confined optical phonons in the general case cannot exist in MoS 2 ! Note that for graphene, since it is a nonpiezoelectric material, confined optical phonon modes do exist. The difference between the two materials lies in Equation 77, specifically in the term containing the piezoelectric coefficient.
If ε zz = 0, continuity of the normal electric displacement component requires, as before, and we return to the situation treated in the previous subsection about confined optical phonons. We need to consider two possible cases: (a) ε xx (ω) = 0 when ε zz (ω) = 0 and (b) ε xx (ω) ≠ 0 when ε zz (ω) = 0. In case (a), we immediately obtain from continuity in the normal electric displacement. can be used and, unless accidental degeneracy applies, the obtained ω m values do not fulfill ε zz (ω) = 0. In conclusion, confined optical phonon modes do not exist in the case ε zz = 0.

Computed dispersion relations
In order to illustrate the validity of the phonon dispersion relations obtained from our model, we compare them to DFT calculations. The continuum theory will require as input elasticity constants, piezoelectric coefficients, and dielectric functions.

DFT
We first give the standard phonon dispersion relation as obtained from DFT calculations (Figure 1). They are obtained from first principles calculations using the Vienna ab initio simulation package (VASP) [17] with a kinetic energy cut-off of 500 eV in the expansion of the electronic wave functions. Four C and six Mo and S valence electrons are considered. The generalized gradient approximation of the exchange-correlation potential in the Perdew-Burke-Ernzerhof flavor is employed [18]. Brillouin-zone integrations are performed using the tetrahedron method with Blöchl corrections [19]. We construct cells consisting of the monolayer and an approximately 15 Å thick vacuum region along the c direction. Structure optimizations are performed on Γ-centered 24 × 24 × 1 k-meshes. A direct method based on 4 × 4 × 1 supercells is used for obtain- ing the phonon dispersions within the harmonic approximation [20]. Forces are evaluated on 3 × 3 × 1 k-meshes, including long-range dipole contributions to the dynamical matrix following the method of [21]. Born effective charges and dielectric tensors are obtained within perturbation theory [22] and elastic constants by the homogeneous deformation method [23]. The elastic stiffness tensor and frequency-dependent dielectric tensor (independent particle approximation [24]) are calculated on Γ-centered 36 × 36 × 1 k-meshes. For the permittivity data we used the DC values ε xx = 4.4ε 0 and ε zz = 1.3ε 0 .
In Figure 2, we show the frequency vs wavenumber (ω-k x ) dispersion in the vicinity of the Γ point (k y = 0). Evidently, one mode shows a parabolic dispersion and one mode is linear. There are two higher-order modes originating from the bound-ary conditions along the z coordinate. The dispersion curves associated with the u y vibrations decoupled from the u x − u z vibrations are not shown in the plots. We emphasize that for graphene optical and acoustic phonon modes decouple and are computed separately as described above. In Figure 3, the dispersion curves for optical phonon modes near the Γ point are shown.  [25] and the piezoelectric constant e x1 is from Ref. [26]. In the case of MoS 2 , all optical and acoustic phonon modes are computed by solving the combined set of elastic and electric equations and associated slab boundary conditions as described above. We note in passing that both confined and surface acousto-optical phonon modes are found using the present formalism. In Figure 4, we show the frequency vs wavenumber (ω-k x ) dispersion in the vicinity of the Γ point (k y = 0).

Conclusion
A three-dimensional, first-principles, continuum elasticity theory was developed for two-dimensional materials. Piezoelectric materials required the simultaneous consideration of the electrostatics equations. Application to graphene, silicene and MoS 2 revealed a number of interesting results. The out-of-plane vibrations were coupled to the in-plane motion for graphene (and for silicene as well), contrary to previous results using infinitely thin sheets. Acoustic modes with linear and quadratic dispersions were obtained, in agreement with experimental results and other models. We predict the existence of confined optical modes in all of the elemental group-IV materials and the nonexistence of confined optical modes in MoS 2 . Our model Note again the parabolic dispersion of one predominantly acoustic mode away from the Γ point besides a linear dispersion predominant acoustic mode. We also obtain a second nonlinear mode starting at (k,ω) = (0,0), which stems from the rather complicated frequency-dependent permittivity of MoS 2 .
will be applied to other 2D materials as well as to multilayers in future work. Additionally, the results of this model can be combined with the solution to the electron problem to compute electron-phonon scattering rates. The latter would be useful for understanding a number of physical properties as well as for applications.