Lattice equations and dielectric function for optical lattice vibrations in monolayer transition metal dichalcogenides

TL;DR

This paper develops a microscopic lattice model for optical vibrations in monolayer transition metal dichalcogenides, deriving explicit equations and dielectric functions that connect phonon properties with first-principles parameters.

Contribution

It introduces a self-consistent dipole lattice model including electronic polarization and local field effects for 2D TMDs, providing explicit equations and relations for optical phonons.

Findings

Derived explicit equations for in-plane and out-of-plane optical vibrations.

Obtained a 2D dielectric function and related phonon dispersion relations.

Studied ionic electronic polarization and local field effects thoroughly.

Abstract

Using a microscopic dipole lattice model including electronic polarization (EP) of ions and local field effects (LFEs) self-consistently, two sets of three equations are deduced for long wavelength in-plane and out-of-plane optical lattice vibrations in monolayer (ML) transition metal dichalcogenides (TMDs). Expressions for the lattice vibrational energy density are obtained for the two-dimensional (2D) crystals. The linear coefficients of the lattice equations can be determined by first-principles calculation, making these equations very useful for studying the lattice dynamical properties. The two pairs of equations describing the polar vibrations have the same forms as those for ML hexagonal BN and also resemble Huang's equations for bulk [Proc. Roy. Soc. A 208, 352 (1951)]. Each pair of the equations are solved simultaneously with the equation of electrostatics, and explicit expressions are obtained for the in-plane longitudinal and transverse optical (LO and TO) modes and out-of-plane optical (ZO) modes. The LO phonon dispersion relation is identical to the analytical expression of Sohier et al. [Nano Lett. 17, 3758 (2017)]. All the optical phonon branches except for LO are nondispersive at the long wavelengths. A 2D longitudinal lattice dielectric function $\epsilon(k,\omega)$ is deduced, allowing one to rederive the LO phonon dispersion simply from $\epsilon(k,\omega)=0$. A 2D Lyddane--Sachs--Teller relation and a frequency--susceptibility relation are obtained for in-plane and out-of-plane vibrations, respectively, through which the phonon frequencies are related to the 2D dielectric functions or susceptibilities. The explicit expressions are applied to calculate various dynamical properties when knowing three first-principles calculated parameters, and the ionic EP and LFEs are studied thoroughly.

Tables

Table 1. Table 1: Various calculated properties associated with in-plane optical vibrations in five ML TMDs and ML hBN, namely, the high-frequency susceptibility χ e subscript 𝜒 𝑒 \chi_{e} , the static susceptibility χ 0 subscript 𝜒 0 \chi_{0} , the absolute values of the Born charge e B subscript 𝑒 𝐵 e_{B} and the static effective charge e 1 subscript 𝑒 1 e_{1} of metal ions ( e 1 = 2 ​ e a subscript 𝑒 1 2 subscript 𝑒 𝑎 e_{1}=2e_{a} for ML TMDs, and e 1 = e a subscript 𝑒 1 subscript 𝑒 𝑎 e_{1}=e_{a} for ML hBN), the effective spring force constant K 1 subscript 𝐾 1 K_{1} , the LO phonon group velocity at the Γ Γ \Gamma point c l subscript 𝑐 𝑙 c_{l} , the atomic polarizability of the unit cell α 𝛼 \alpha , the total free-atom polarizability of the constituent atoms α f subscript 𝛼 𝑓 \alpha_{f} and α / s 𝛼 𝑠 \alpha/s ( s 𝑠 s is the unit cell area). These properties are calculated with the expressions derived in this article and using the parameters of first-principles calculation in Table 1 of Ref. Sohier  et al. ( 2017 ) (see text).

ML	${\chi }_e$ ($\AA$)	${\chi }_0$ ($\AA$)	$\left|e_B\right|$ ($e$)	$\left|e_1\right|$ ($e$)	$K_1$ (eV/${\AA }^2$)	$c_l$ (km/s)	$\alpha$ (${\AA }^3$)	${\alpha }_f$ (${\AA }^3$)	$\alpha /s$ ($\AA$)
MoS2	7.40	7.48	1.00	0.22	9.957	1.85	[5.25,10.54]	18.55	[0.61,1.22]
MoSe2	8.47	8.75	1.76	0.36	8.673	4.61	[6.02,12.18]	20.52	[0.64,1.29]
MoTe2	11.06	12.07	3.10	0.54	6.891	13.34	[7.41,15.20]	23.71	[0.69,1.42]
WS2	6.68	6.70	0.48	0.12	10.481	0.37	[5.18,10.31]	16.92	[0.60,1.20]
WSe2	7.75	7.86	1.15	0.26	9.053	1.60	[5.91,11.87]	18.89	[0.63,1.27]
hBN	1.22	1.67	2.71	0.46	56.495	37.10	[1.43,1.98]	4.0	[0.26,0.37]

Equations

$$\phi (z)={\int }_{-\infty }^{\infty }\hat{\phi }(q)e^{iqz}dq,$$

$$\delta (z-z_{\nu })=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{iq(z-z_{\nu })}dq.$$

$${\mathbf{E}}_{\mathbf{\rho }}(\mathbf{\rho },z)=-\frac{2\pi }{3}\frac{\mathbf{k}}{k}\mathcal{P}(\mathbf{\rho })\cdot \mathbf{k}\sum _{\nu }{e}^{-k|z-z_{\nu }|},$$

$${\mathbf{E}}_z(\mathbf{\rho },z)=-{\mathbf{e}}_z\frac{2\pi i}{3}\mathcal{P}(\mathbf{\rho })\cdot \mathbf{k}\sum _{\nu }\mathrm{sgn}(z-z_{\nu })e^{-k|z-z_{\nu }|}.$$

$$\mathbf{E}_{\bm{\rho}}(\bm{\rho},z)=-\frac{2\pi}{s}\frac{\mathbf{k}}{k}~{}\Big{[}\sum_{\nu}\mathbf{p}_{\nu}\cdot\mathbf{k}e^{-k\lvert z-z_{\nu}\rvert}\Big{]}e^{i\mathbf{k}\cdot\bm{\rho}},$$

$$\mathbf{E}_{z}(\bm{\rho},z)=-\mathbf{e}_{z}\frac{2\pi i}{s}\Big{[}\sum_{\nu}\mathbf{p}_{\nu}\cdot\mathbf{k}\operatorname{sgn}(z-z_{\nu})e^{-k\lvert z-z_{\nu}\rvert}\Big{]}e^{i\mathbf{k}\cdot\bm{\rho}}.$$

$${\mathbf{E}}_{l,1}=\mathbf{E}+Q_0{\mathbf{p}}_1+Q_1{\mathbf{p}}_2+Q_1{\mathbf{p}}_{\bar{2}},$$

$${\mathbf{E}}_{l,2}=\mathbf{E}+Q_1{\mathbf{p}}_1+Q_0{\mathbf{p}}_2+Q_2{\mathbf{p}}_{\bar{2}},$$

$${\mathbf{E}}_{l,\bar{2}}=\mathbf{E}+Q_1{\mathbf{p}}_1+Q_2{\mathbf{p}}_2+Q_0{\mathbf{p}}_{\bar{2}}.$$

$$Q_0=\sum _{(m,n)\ne (0,0)}\frac{1}{2(m^2+n^2+mn{)}^{3/2}{a}^3}\approx \frac{5.5171}{a^3},$$

$$Q_1$$

$$\displaystyle\qquad{}\Big{.}+\frac{3(m^{2}+n^{2}+mn+n+1/3)}{2[m^{2}+n^{2}+mn+n+1/3+(d/a)^{2}/4]^{5/2}}\Big{\}}\frac{1}{a^{3}},$$

$$Q_2$$

$$\displaystyle\qquad{}\Big{.}+\frac{3(m^{2}+n^{2}+mn)}{2[m^{2}+n^{2}+mn+(d/a)^{2}]^{5/2}}\Big{\}}\frac{1}{a^{3}}.$$

$${\mathbf{E}}_{u,1}=-{\mathbf{u}}_1(e_1{Q}_0+e_2{Q}_1+e_{\bar{2}}{Q}_1),$$

$${\mathbf{E}}_{u,2}=-{\mathbf{u}}_2(e_1{Q}_1+e_2{Q}_0+e_{\bar{2}}{Q}_2),$$

$${\mathbf{E}}_{u,\bar{2}}=-{\mathbf{u}}_{\bar{2}}(e_1{Q}_1+e_2{Q}_2+e_{\bar{2}}{Q}_0).$$

$${\mathbf{E}}_1=\mathbf{E}+Q_0{\mathbf{p}}_1+Q_1{\mathbf{p}}_2+Q_1{\mathbf{p}}_{\bar{2}}+2e_a(Q_1-Q_0){\mathbf{u}}_1,$$

$${\mathbf{E}}_2=\mathbf{E}+Q_1{\mathbf{p}}_1+Q_0{\mathbf{p}}_2+Q_2{\mathbf{p}}_{\bar{2}}+e_a(Q_0+Q_2-2Q_1){\mathbf{u}}_2,$$

$${\mathbf{E}}_{\bar{2}}=\mathbf{E}+Q_1{\mathbf{p}}_1+Q_2{\mathbf{p}}_2+Q_0{\mathbf{p}}_{\bar{2}}+e_a(Q_0+Q_2-2Q_1){\mathbf{u}}_{\bar{2}},$$

$${\mathbf{p}}_1=2e_a{\mathbf{u}}_1+{\alpha }_1{\mathbf{E}}_1,$$

$${\mathbf{p}}_2=-e_a{\mathbf{u}}_2+{\alpha }_2{\mathbf{E}}_2,$$

$${\mathbf{p}}_{\bar{2}}=-e_a{\mathbf{u}}_{\bar{2}}+{\alpha }_2{\mathbf{E}}_{\bar{2}}.$$

$$(1-{\alpha }_1{Q}_0){\mathbf{p}}_1-{\alpha }_1{Q}_1({\mathbf{p}}_2+{\mathbf{p}}_{\bar{2}})$$

$$=2e_a[1+{\alpha }_1(Q_1-Q_0)]{\mathbf{u}}_1+{\alpha }_1\mathbf{E},$$

$$-{\alpha }_2{Q}_1{\mathbf{p}}_1+(1-{\alpha }_2{Q}_0){\mathbf{p}}_2-{\alpha }_2{Q}_2{\mathbf{p}}_{\bar{2}}$$

$$=-e_a[1-{\alpha }_2(Q_0+Q_2-2Q_1)]{\mathbf{u}}_2+{\alpha }_2\mathbf{E},$$

$$-{\alpha }_2{Q}_1{\mathbf{p}}_1-{\alpha }_2{Q}_2{\mathbf{p}}_2+(1-{\alpha }_2{Q}_0){\mathbf{p}}_{\bar{2}}$$

$$=-e_a[1-{\alpha }_2(Q_0+Q_2-2Q_1)]{\mathbf{u}}_{\bar{2}}+{\alpha }_2\mathbf{E}.$$

$${\mathbf{p}}_1$$

$$\displaystyle\qquad{}-2e_{a}\alpha_{1}Q_{1}\big{[}1-\alpha_{2}(Q_{0}+Q_{2}-2Q_{1})\big{]}\mathbf{u}_{c}$$

$$\displaystyle\qquad{}\Big{.}+\alpha_{1}\big{[}1-\alpha_{2}(Q_{0}+Q_{2}-2Q_{1})\big{]}\mathbf{E}\Big{\}},$$

$${\mathbf{p}}_2$$

$$\displaystyle\qquad{}-e_{a}(1-\alpha_{1}Q_{0})\big{[}1-\alpha_{2}(Q_{0}+Q_{2}-2Q_{1})\big{]}\mathbf{u}_{c}$$

$$\displaystyle\qquad{}\Big{.}+\alpha_{2}\big{[}1+\alpha_{1}(Q_{1}-Q_{0})\big{]}\mathbf{E}\Big{\}}-\frac{1}{2}e_{a}\eta\mathbf{u}_{d},$$

$${\mathbf{p}}_{\bar{2}}$$

$$\displaystyle\qquad{}-e_{a}(1-\alpha_{1}Q_{0})\big{[}1-\alpha_{2}(Q_{0}+Q_{2}-2Q_{1})\big{]}\mathbf{u}_{c}$$

$$\displaystyle\qquad{}\Big{.}+\alpha_{2}\big{[}1+\alpha_{1}(Q_{1}-Q_{0})\big{]}\mathbf{E}\Big{\}}+\frac{1}{2}e_{a}\eta\mathbf{u}_{d},$$

$${\mathbf{u}}_d={\mathbf{u}}_2-{\mathbf{u}}_{\bar{2}},$$

$${\mathbf{u}}_c=({\mathbf{u}}_2+{\mathbf{u}}_{\bar{2}})/2,$$

$$D=(1-{\alpha }_1{Q}_0)[1-{\alpha }_2(Q_0+Q_2)]-2{\alpha }_1{\alpha }_2{Q}_1^2,$$

$$\eta =\frac{1-{\alpha }_2(Q_0+Q_2-2Q_1)}{1-{\alpha }_2(Q_0-Q_2)}.$$

$$\mathbf{w}=\sqrt{\frac{\bar{m}}{s}}({\mathbf{u}}_1-{\mathbf{u}}_c)=\sqrt{\frac{\bar{m}}{s}}\left[{\mathbf{u}}_1-\frac{1}{2}({\mathbf{u}}_2+{\mathbf{u}}_{\bar{2}})\right],$$

$$\mathcal{P}=a_{21}\mathbf{w}+a_{22}\mathbf{E},$$

Full text

I Introduction

Interest in two-dimensional (2D) materials such as monolayer (ML) transition metal dichalcogenides (TMDs) has surged in recent years due to their potential for novel electronic and optoelectronic devices Geim and Grigorieva (2013); Wang et al. (2012); Mak et al. (2012); Fiori et al. (2014). The ML TMDs such as MoS2 are 2D polar crystals with a honeycomb lattice structure, in which the polar lattice vibrations occur together with the nonpolar vibrations. Particularly interesting are the long wavelength polar optical vibrations; for instance, the longitudinal optical (LO) vibrations can generate a macroscopic electric field that strongly couples to electrons, and this makes the polar electron-phonon interaction dominate the electron transport processes in 2D TMDs. The 2D optical phonons have attracted great interest as they play a key role in carrier transport Kaasbjerg et al. (2012); Li et al. (2013); Danovich et al. (2017), piezoelectricity Wu et al. (2014); Zhu et al. (2014); Michel et al. (2017), valley depolarization Miller et al. (2019) and homogeneous broadening of excitonic luminescence Moody et al. (2015); Jakubczyk et al. (2016) in ML TMDs. Hence, an in-depth knowledge of the 2D lattice vibrations is very useful for the study of the electronic and optical properties of ML TMDs.

A generally used method for phonon calculations is the numerical solution of the dynamical matrix Born and Huang (1954), and with it the phonon spectra in the entire Brillouin zone can be obtained. For isotropic bulk crystals the macroscopic theory developed by Huang allows for an analytical description of the long optical lattice vibrations by a pair of equations, $\ddot{\mathbf{w}}=b_{11}\mathbf{w}+b_{12}\mathbf{E}$, $\mathbf{P}=b_{21}\mathbf{w}+b_{22}\mathbf{E}$ ($b_{12}=b_{21}$) Huang (1951a, b); Born and Huang (1954), where $\mathbf{w}$ is used for describing the optical displacement of the lattice, and $\mathbf{E}$ and $\mathbf{P}$ are the electric field and dielectric polarization as in Maxwell’s equations, both being macroscopic quantities, averaged over a lattice cell Born and Huang (1954); Kittel (2004). Compared to the Hooke’s law formula, the right-hand side of the equation of motion has an extra term due to the electric force, while in the polarization equation there is a lattice displacement contribution as well as the usual field-induced polarization. Microscopically, the field acting on any particular ion is a local field (Lorentz field), i.e., the total field of all the ions except for the ion itself, but it is eliminated by means of the Lorentz relation Born and Huang (1954); Kittel (2004), which is essential in Huang’s deduction. The electronic polarization (EP) of ions is also accounted for in Huang’s model. Using Huang’s equations not only the long wavelength optical modes but also the lattice dielectric function is obtained, and further the well-known Lyddane–Sachs–Teller (LST) relation Lyddane et al. (1941) is rederived Born and Huang (1954).

Like in bulk polar crystals, the long wavelength lattice vibrations in a 2D polar crystal are closely connected with the macroscopic and local fields in the crystal. In an isotropic bulk crystal no macroscopic field occurs with the transverse optical (TO) vibrations but the LO vibrations generate a finite long-range macroscopic field, leading to a higher LO phonon frequency and LO-TO splitting Lyddane et al. (1941). In a 2D crystal, the lattice may polarize according to the in-plane or out-of-plane lattice motion, giving rise to the macroscopic fields that differ for the in-plane and out-of-plane polarization Zhang (2019). While there is no macroscopic field accompanying the TO vibrations, like in the bulk, the macroscopic field of the LO vibrations differs substantially from that in the bulk, as it vanishes in the long wavelength limit, leading to the degeneracy of the LO and TO modes at the $\Gamma$ point Sánchez-Portal and Hernández (2002). Further, the macroscopic field due to the out-of-plane vibrations becomes extremely large close to the atomic layers as the lattice polarization is of only atomic scale in the out-of-plane direction $z$ Zhang (2019). Similarly, the local field in the 2D crystal at an ion site strongly depends on the lattice motion, making the local fields differ for the in-plane and out-of-plane lattice vibrations Mikhailov (2013). In fact these are very strong local fields Mikhailov (2013); Zhang (2019), which is different from the case of the bulk where the local fields on the ions are quite weak and largely cancelled Born and Huang (1954).

As the wavevector is parallel to the layer planes, the lattice vibrations of 2D crystals comprise the in-plane and out-of-plane vibrational modes. In ML hBN there are three branches of optical modes, all of which are polar modes. With one positive ion and two same negative ions per unit cell, a ML TMD has three branches of nonpolar optical modes corresponding to the relative motion the two negative ions, as well as three branches of polar optical modes like those of ML hBN, which are connected with the contrary motion of the positive and negative ions. Specifically these optical branches are composed of the in-plane nonpolar LO and TO modes (labeled LO1, TO1), the in-plane polar LO and TO modes (labeled LO2, TO2), the out-of-plane nonpolar modes, ZO1, and the out-of-plane polar modes ZO2. At the $\Gamma$ point of the Brillouin zone the LO1 and TO1 are degenerate and denoted as $E^{\prime\prime}$ by symmetry consideration, the doubly degenerate LO2 and TO2 modes are denoted by $E^{\prime}$, and the out-of-plane ZO1 and ZO2 modes are denoted as $A_{1}$ and $A_{2}^{\prime\prime}$, respectively Moliná-Sanchez and Wirtz (2011). $A_{2}^{\prime\prime}$ is an infrared active mode, while $A_{1}$ and $E^{\prime\prime}$ are Raman-active modes, and the $E^{\prime}$ modes are both infrared and Raman active. The $A_{1}$ and $E^{\prime}$ modes have been measured by Raman spectroscopy for a number of ML TMDs Lee et al. (2010); Tongay et al. (2012); Lanzillo et al. (2013); Rice et al. (2013); Saito et al. (2016); Lu et al. (2016).

Phonon dispersion relations of ML TMDs have been studied based on first-principles calculations Moliná-Sanchez and Wirtz (2011); Ataca et al. (2011); Sohier et al. (2017). Phonon modes calculated from first principles have been used to study mechanical stress-strain relations Li (2012), electron transport properties Kaasbjerg et al. (2012); Li et al. (2013), thermal conductivity Cai et al. (2014) and thermal expansion Cai et al. (2014); Huang et al. (2014); Çakir et al. (2014) and piezoelectricity Michel et al. (2017) in ML TMDs. These calculated phonon spectra show several common features on the optical branches: (i) apart from the degeneracy of the in-plane nonpolar modes LO1 and TO1, the polar modes LO2 and TO2 are also degenerate at the $\Gamma$ point; (ii) the LO2 and TO2 modes split at finite wavevectors, and the slow increase of LO phonon frequency leads to a small overbending (i.e., the maximum frequency occurs at a wavevector away from $\Gamma$ Sohier et al. (2017)) in the LO phonon dispersion compared to ML hBN Topsakal et al. (2009); Şahin et al. (2009); Sohier et al. (2017); (iii) the ZO1 nonpolar branch displays a nondispersive behaviour in the Brillouin zone. While most of the theoretical studies are based on first-principles calculations, there is a lack of analytical approaches and there are no equations of motion to describe the long optical lattice vibrations in ML TMDs. Using a simple model, Sohier et al. have recently derived an LO phonon dispersion relation Sohier et al. (2017) for the 2D polar crystals, including hBN and TMDs, and the parameters in the dispersion relation can be calculated from first-principles to include ionic EP and local field effects (LFEs). The dispersion relation can well describe the dispersion of the numerically calculated polar optical modes Sohier et al. (2017).

Recently we have deduced lattice equations for the long wavelength optical vibrations in ML hBN using a microscopic model, and then derived optical phonon dispersion relations in which the LO phonon dispersion relation is identical to the expression of Sohier et al. Sohier et al. (2017). Here one purpose of this study is to extend our approaches to analytically tackle the lattice vibrations of ML TMDs. To formulate a macroscopic theory of the optical vibrations, both macroscopic quantities, i.e., the macroscopic field to drive the lattice motion and the 2D dielectric polarization are needed Zhang (2019). They both need to be defined at first, and a digestion of the results on ML hBN Zhang (2019) will give us some guidance. The macroscopic field act on all the ions of a lattice cell so in principle it should have the same value for all these ions Born and Huang (1954); Kittel (2004). As ML hBN has only a single atomic layer, one of course uses the macroscopic field in the ML plane, and also defines the 2D lattice polarization straightforward Zhang (2019). Unlike ML hBN, however a ML TMD consists of three atomic layers at different $z_{i}$ ($i=1,2,\bar{2}$) positions and has a finite thickness. As each of the three layers resembles a hBN monolayer, using the macroscopic field result of Ref.Zhang (2019) each in-plane polarized layer produces a macroscopic field that varies with $z$ according to $e^{-k\lvert z-z_{i}\rvert}$ ($k$ is a small wavevector), which varies very slowly inside the ML because for any $z$ in the ML $k\lvert z-z_{i}\rvert$ is a small quantity, $k\lvert z-z_{i}\rvert\ll 1$. Therefore we use the macroscopic field at the central layer $z_{1}$ (positive ion site) to represent the macroscopic field acting on all the three layers (all the three ions of a cell). Further, this shows that the different $z_{i}$ positions of the layers has become unimportant and the 2D TMD behaves like a single layer crystal to the long lattice waves. Therefore we can define a 2D macroscopic polarization from the microscopic electric dipole density, by averaging its microscopic $\delta(z-z_{i})$ distributions at the three layers (Sec. II below). For the out-of-plane polarization, the electric fields of the three polarized atomic layers are equally divergent according to $\delta(z-z_{i})$ Zhang (2019) so the averaging treatment can be also used to handle the out-of-plane lattice motion.

Having clarified these points regarding the macroscopic field and polarization, in this paper, therefore we can extend the approaches of the previous work Zhang (2019) to study systematically the long wavelength optical vibrations in ML TMDs. We deduce lattice equations involving both macroscopic quantities using a microscopic model that accounts for ionic EP and LFEs self-consistently. Then we solve the lattice equations simultaneously with the equation of electrostatics to obtain explicit polar optical dispersion relations. In particular, we obtained an LO phonon dispersion relation which is identical to the expression of Sohier et al. Sohier et al. (2017). We also derive a 2D longitudinal lattice dielectric function, and further a 2D LST relation for in-plane motion together with a frequency–susceptibility relation for out-of-plane lattice motion. Further we use first-principles calculated quantities to study the lattice dynamical properties and also local field and polarizable ion effects in ML TMDs.

This paper is organized as follows. In Section II, a deduction of two pairs of three lattice equations for in-plane and out-of-plane motion is made from a microscopic dipole lattice model. Then the dispersion relations of all the optical branches are obtained together with the phonon group velocity and density of states. Further the lattice dielectric susceptibility and longitudinal lattice dielectric function are derived and discussed. In Section III, we present numerical results of the in-plane optical vibrations in ML TMDs. Using the first-principles calculated parameters of Ref.Sohier et al. (2017), the various lattice dynamical properties including the Born charge and the static and electronic susceptibilities are calculated and compared for several ML TMDs, and further the EP and LFEs are studied in detail. Finally, Section V summarizes the main results obtained.

II Theory

II.1 Equations of motion and lattice polarization

In this section we use ML MoS2 as representative to present our theory, but the results are applicable to all ML TMDs. Monolayer MoS2 is a 2D binary crystal (with D3h symmetry) with a hexagonal lattice structure (Fig. 2), which is composed of one sublattice of Mo and two sublattices of S, labeled with the base index $\nu=1,2,\bar{2}$, respectively. Let $a$ and $d$ be the lattice constant and the vertical separation between the two S layers, $d\approx a$ Ataca et al. (2011); Cao et al. (2012); Chang et al. (2013), and a lattice cell has an area $s=\sqrt{3}a^{2}/2$. Let $m_{\nu}$ and $e_{\nu}$ be the mass and charge of the type $\nu$ ions, and let $e_{2}=e_{\bar{2}}$ (i.e., for the two S layers), $e_{1}=-2e_{2}=2e_{a}$, where $2e_{a}$ is the static effective charge Karch and Bechstedt (1997) on a Mo atom due to the electron charge transfer. Let $\bm{\rho}_{l}$ be the 2D lattice vector, $\bm{\rho}_{l}=l_{1}\mathbf{a}_{1}+l_{2}\mathbf{a}_{2}$, $l$ being the cell index, $l=(l_{1},l_{2})$, and let $\mathbf{R}_{l\nu}$ be the position vector of the $\nu$ ion in the $l$-cell, $\mathbf{R}_{l\nu}=\bm{\rho}_{l}+\mathbf{r}_{\nu}$, where $\mathbf{r}_{\nu}=(\bm{\rho}_{\nu},z_{\nu})$, with $z_{1}=0$, $z_{2}=d/2$, $z_{\bar{2}}=-d/2$.

In the dipole lattice model Born and Huang (1954) each ion site of type $\nu$ is occupied by an electric dipole $\mathbf{p}_{\nu}$ that is due partly to the ionic displacement $\mathbf{u}_{\nu}$ and partly to the induced electric moment $\bm{\mu}_{\nu}$ on the ion. There is polarization $\mathbf{P}_{l}=\sum_{\nu}(\mathbf{p}_{\nu}/s)e^{i\mathbf{k}\cdot\mathbf{R}_{l\nu}}\delta(z-z_{\nu})$ associated with a lattice wave of 2D wavevector $\mathbf{k}$. For a long wavelength lattice wave, $\mathbf{k}\cdot\mathbf{R}_{l\nu}\approx\mathbf{k}\cdot\bm{\rho}_{l}$, then the lattice can be treated as a polarized continuum with the polarization $\mathbf{P}_{c}=[\sum_{\nu}\mathbf{p}_{\nu}\delta(z-z_{\nu})]e^{i\mathbf{k}\cdot\bm{\rho}}/s$, where the discrete lattice vector $\bm{\rho}_{l}$ has become a continuous variable $\bm{\rho}=(x,y)$. Here the $\delta$ function $\delta(z-z_{\nu})$ describes the microscopic distribution (i.e., at the atomic scale) in the vertical direction of the dipole moment density of the $\nu$ sublattice – such a microscopic distribution exists as there is no translational symmetry in the $z$ direction. When conditions are uniform over many unit cells (for small $k$) the individual layers’ specific positions $z_{\nu}$ become unimportant. To have a macroscopic quantity to effectively describe the 2D lattice polarization, therefore, we neglect the difference between the three microscopic $\delta(z-z_{\nu})$ distributions associated with the three layers of dipoles, by approximating the three $\delta$ functions with a single function independent of the base index $\nu$ which should also capture the microscopic character of the polarization’s $z$-distribution. We take the function to be an average of the three $\delta$ functions, $\bar{\delta}(z)=\sum_{\nu}\delta(z-z_{\nu})/3$, symmetric with respect to the central Mo layer due to symmetry of the 2D crystal and normalized to unity, such that the 3D polarization of the dielectric is given by $\mathbf{P}=(\sum_{\nu}\mathbf{p}_{\nu}/s)e^{i\mathbf{k}\cdot\bm{\rho}}\bar{\delta}(z)$. Evidently $\bm{\mathcal{P}}(\bm{\rho})=(\sum_{\nu}\mathbf{p}_{\nu}/s)e^{i\mathbf{k}\cdot\bm{\rho}}$ is a macroscopic quantity Born and Huang (1954) (average dipole moment per unit area of the cell, with $\bm{\mathcal{P}}_{0}=\sum_{\nu}\mathbf{p}_{\nu}/s$ hereafter) for 2D polarization whilst $\bar{\delta}(z)$ characterizes the volume polarization $\mathbf{P}$’ microscopic distribution in the $z$ direction. Physically, the averaging $\bar{\delta}(z)$ plays a role of projecting the upper and lower layers of S atoms on the central Mo layer, thus producing a similar hexagonal lattice to ML hBN [Fig. 2(a)]; this is an effective treatment for the dielectric polarization of the three-layer 2D crystals, which will become more apparent when we obtain the macroscopic fields, lattice modes and dielectric function below.

We first consider in-plane motion, with the displacements $\mathbf{u}_{\nu}$ and dipole moments $\mathbf{p}_{\nu}$ parallel to the ML. The equation of electrostatics is $\nabla\cdot(\mathbf{E}+4\pi\mathbf{P})=0$, where $\mathbf{E}$ is an irrotational macroscopic field, $\mathbf{E}=-\nabla\phi$, $\phi$ being the electric potential, $\phi(\bm{\rho},z)=\varphi(z)e^{i\mathbf{k}\cdot\bm{\rho}}$. To solve Poisson’s equation $\nabla^{2}\phi(\bm{\rho},z)=4\pi i\bm{\mathcal{P}}(\bm{\rho})\cdot\mathbf{k}\bar{\delta}(z)$, one uses the expansions,

[TABLE]

[TABLE]

One finds $\hat{\varphi}(q)=-2i\bm{\mathcal{P}}_{0}\cdot\mathbf{k}/[3(k^{2}+q^{2})][1+2\cos(qd/2)]$ and from Eq. (1) $\varphi(z)=-2\pi i\bm{\mathcal{P}}_{0}\cdot\mathbf{k}\sum_{\nu}e^{-k\lvert z-z_{\nu}\rvert}/(3k)$, and then obtains the parallel and perpendicular components of the macroscopic field:

[TABLE]

To contrast, when $\mathbf{P}_{c}$ is used for the polarization, the electric field is given by

[TABLE]

Comparing Eqs. (3a) and (4a), and also Eqs. (3b) and (4b) we see that the use of $\bar{\delta}(z)$, i.e., the averaging of the $\delta$ functions leads to the $z$-dependences of the field being substituted by the averages, that is, $e^{-k\lvert z-z_{\nu}\rvert}$ of $\mathbf{E}_{\bm{\rho}}(\mathbf{r})$ in Eq. (4a) is replaced by $\sum_{\nu}e^{-k\lvert z-z_{\nu}\rvert}/3$ of Eq. (3a), while $\operatorname{sgn}(z-z_{\nu})e^{-k\lvert z-z_{\nu}\rvert}$ of $\mathbf{E}_{z}(\mathbf{r})$ in Eq. (4b) is substituted by $\sum_{\nu}\operatorname{sgn}(z-z_{\nu})e^{-k\lvert z-z_{\nu}\rvert}/3$ of Eq. (3b), as is expected according to the principle of superposition of electric fields. As a result, the macroscopic field in the form of Eqs. (3a) and (3b) is proportional to the 2D macroscopic polarization $\bm{\mathcal{P}}(\bm{\rho})$ defined above, as is needed, and further it has captured the limiting behaviour of the electric field of Eqs. (4a) and (4b), i.e., $\mathbf{E}\rightarrow 0$, as $\mathbf{k}\rightarrow 0$. As $k\lvert z-z_{\nu}\rvert$ is very small at the long wavelengths, $k\lvert z-z_{\nu}\rvert\leq kd\ll 1$, each $e^{-k\lvert z-z_{\nu}\rvert}$ ($\nu=1,2,\bar{2}$) varies very slowly in the ML and therefore the averaging is a good treatment for the in-plane field $\mathbf{E}_{\bm{\rho}}(\mathbf{r})$. The out-of-plane field $\mathbf{E}_{z}(\mathbf{r})$ of Eq. (3b) is zero at $z=0$, and antisymmetric with respect to the central layer on which the upper and lower layers of ions are projected through the averaging $\bar{\delta}(z)$, and thus $\mathbf{E}_{z}(\mathbf{r})$ has no influence on the ionic motion. Nevertheless, only the in-plane field $\mathbf{E}_{\bm{\rho}}(\mathbf{r})$ is useful, used by the lattice equations, because it is this field that drives the lattice motion. Therefore the long wavelength macroscopic field $\mathbf{E}_{\bm{\rho}}(\mathbf{r})$ can be considered uniform within a lattice cell, having the same value on all the ions in it. Here the value of the field is taken at the center of the Mo layer, $z=0$ Cudazzo et al. (2011), and we use $\mathbf{E}=\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)=\mathbf{E}(\bm{\rho},0)$ for short.

The Lorentz local field $\mathbf{E}_{l}$ (also called the exciting field Born and Huang (1954)), which is the field acting on an ion due to all the other dipoles in a dipole lattice, is given by the sum of the macroscopic field $\mathbf{E}$ and inner field $\mathbf{E}_{in}$ Born and Huang (1954), $\mathbf{E}_{l}=\mathbf{E}+\mathbf{E}_{in}$. In the long wavelength lattice waves, the local fields at the Mo and S ion sites are given by the following expressions Mikhailov (2013),

[TABLE]

with the coefficients $Q_{0}$, $Q_{1}$ and $Q_{2}$

[TABLE]

Note that these $Q$-coefficients and also inner fields are calculated rigorously using the real three-layer crystal structure. With the approximation $d=a$, $Q_{1}$ and $Q_{2}$ are simplified to $Q_{1}\approx 1.7184/a^{3}$ and $Q_{2}\approx-0.1127/a^{3}$.

Eqs. (5a), (5b) and (5c) are 2D Lorentz relations for ML MoS2 ($\mathbf{E}$ is the macroscopic field and not an external field in a simple sense Mikhailov (2013)). These relations are valid only for the long lattice waves as in bulk but they have a more complicated form than the 3D Lorentz relation Born and Huang (1954); Kittel (2004) because of the different $Q_{i}$ ($i=0,1,2$) coefficients.

To find the Coulomb field at a $\nu$ ion we also need to consider the field change at the center of the ion due to its own displacement $\mathbf{u}_{\nu}$ Born and Huang (1954). Evidently the field is equal to the field created at the $\nu$ ion site by displacing all other ions by $-\mathbf{u}_{\nu}$, and therefore equal to the local field at that ion site in a dipole lattice with displacement dipoles $\mathbf{p}_{\nu^{\prime}}=-e_{\nu^{\prime}}\mathbf{u}_{\nu}$, where $\nu^{\prime}=1,2,\bar{2}$. Inserting this dipole expression into Eqs. (5a), (5b) and (5c) and putting $\mathbf{E}=0$, one finds the field changes at the centers of the Mo and S ions due to their displacements,

[TABLE]

The total Coulomb fields $\mathbf{E}_{1}$, $\mathbf{E}_{2}$ and $\mathbf{E}_{\bar{2}}$ at the centers of the Mo and S ions are the sums of $\mathbf{E}_{l,1}$ and $\mathbf{E}_{u,1}$ [Eqs. (5a) and (7a)], $\mathbf{E}_{l,2}$ and $\mathbf{E}_{u,2}$ [Eqs. (5b) and (7b)], and $\mathbf{E}_{l,\bar{2}}$ and $\mathbf{E}_{u,\bar{2}}$ [Eqs. (5c) and (7c)], respectively,

[TABLE]

where all the vectors are parallel to the ML.

The electronic polarization of an ion is equivalent to a point-dipole so the induced dipole moment of a $\nu$ ion is $\bm{\mu}_{\nu}=\alpha_{\nu}\mathbf{E}_{\nu}$ Born and Huang (1954), where $\alpha_{\nu}$ is the in-plane electronic polarizability of the $\nu$ ions ($\alpha_{2}=\alpha_{\bar{2}}$). Then the total dipole moments on the Mo and S ions are

[TABLE]

When the expressions for the total fields $\mathbf{E}_{1}$, $\mathbf{E}_{2}$ and $\mathbf{E}_{\bar{2}}$ are inserted, one finds

[TABLE]

Solving these equations, one can express $\mathbf{p}_{1}$, $\mathbf{p}_{2}$ and $\mathbf{p}_{\bar{2}}$ in terms of $\mathbf{u}_{1}$, $\mathbf{u}_{2}$, $\mathbf{u}_{\bar{2}}$ and $\mathbf{E}$ as follows:

[TABLE]

where $\mathbf{u}_{d}$ and $\mathbf{u}_{c}$ are the relative displacement of the two S atoms and the displacement of their centre of mass, respectively,

[TABLE]

and

[TABLE]

[TABLE]

Introducing the optical displacement $\mathbf{w}$ to describe the motion of the Mo atom relative to the two S atoms,

[TABLE]

with $\bar{m}$ being the reduced mass, $\bar{m}=2m_{1}m_{2}/(m_{1}+2m_{2})$, one finds that the areal polarization $\bm{\mathcal{P}}=(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{\bar{2}})/s$ is given by

[TABLE]

where

[TABLE]

Eq. (16) shows that the 2D macroscopic polarization $\bm{\mathcal{P}}$ of ML MoS2 reduces to a sum of two contributions, one due to the optical displacement and the other due to the macroscopic field.

It is evident from Eq. (16) that $a_{22}$ is the in-plane electronic susceptibility $\chi_{e}$ of 2D MoS2,

[TABLE]

Expression (17b) relates the dielectric susceptibility to the polarizabilities of the constituent atoms, and is a Clausius-Mossotti relation. Introducing the Born charge Gonze and Lee (1997)

[TABLE]

then the coefficient $a_{21}$ relates to $e_{B}$ through

[TABLE]

Using Eqs. (13), (17b), (18) and (19) one obtains

[TABLE]

showing that the $a_{22}$ coefficient is also related to the Born charge. $e_{B}\neq 2e_{a}$ owing to the EP of ions, i.e., $\chi_{e}>0$, and as $Q_{1}>0$, $e_{B}$ and $e_{a}$ have the same sign and $\lvert e_{B}\rvert>2\lvert e_{a}\rvert$.

Insert expressions (11a), (11b) and (11c) for $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ and $\mathbf{p}_{\bar{2}}$ into Eqs. (8a), (8b) and (8c) and express the total fields $\mathbf{E}_{1}$, $\mathbf{E}_{2}$ and $\mathbf{E}_{\bar{2}}$ in terms of the ionic displacements and macroscopic field. After a lengthy algebraic simplification the total field expressions are simplified to

[TABLE]

Further one readily finds from Eqs. (22b), (22c)

[TABLE]

The Coulomb force on a $\nu$ ion consists of two parts Born and Huang (1954), (i) the force exerted on ionic charge $e_{\nu}$ by total field $\mathbf{E}_{\nu}$, and (ii) the force exerted on induced dipole $\bm{\mu}_{\nu}$ by the field of all other ions. The latter can be obtained by using again the Lorentz relations (5a), (5b) and (5c). Subject dipole $\bm{\mu}_{\nu}$ at a particular site of ion $\nu$ a virtual displacement $\mathbf{u}$ whilst fixing all other ions in their equilibrium positions. The virtual energy is the interaction energy between the dipole and the field created equivalently at the ion site by displacing all other ions by $-\mathbf{u}$ Born and Huang (1954), which is exactly the local field in the dipole lattice with dipole moments $\mathbf{p}_{\nu^{\prime}}=-e_{\nu^{\prime}}\mathbf{u}$, ($\nu^{\prime}=1,2,\bar{2}$). Therefore the virtual energy is $-\bm{\mu}_{\nu}\cdot\mathbf{E}_{l,\nu}$, resulting in the force on the dipole $\nabla_{\mathbf{u}}(\bm{\mu}_{\nu}\cdot\mathbf{E}_{l,\nu})$. Inserting the dipole moment above into Eqs. (5a), (5b) and (5c) and putting $\mathbf{E}=0$, one obtains the local fields and then the forces on the dipoles $\bm{\mu}_{1}$, $\bm{\mu}_{2}$ and $\bm{\mu}_{\bar{2}}$ given by $2e_{a}(Q_{1}-Q_{0})\bm{\mu}_{1}$, $e_{a}(Q_{0}+Q_{2}-2Q_{1})\bm{\mu}_{2}$, and $e_{a}(Q_{0}+Q_{2}-2Q_{1})\bm{\mu}_{\bar{2}}$, respectively, where $\bm{\mu}_{\nu}=\alpha_{\nu}\mathbf{E}_{\nu}$ ($\nu=1,2,\bar{2}$).

To describe the restoring forces, let $K_{1}$ be the spring force constant between the Mo and S ions, and $K_{2}$ be the spring force constant between the two S layers. Therefore the equations of motion for the Mo and S ions are given by

[TABLE]

Substitute expression (22a) for $\mathbf{E}_{1}$ in Eq. (24a). Add Eqs. (24b) and (24c), change $\mathbf{u}_{2}+\mathbf{u}_{\bar{2}}$ to $2\mathbf{u}_{c}$, and substitute expression (23a) for $\mathbf{E}_{2}+\mathbf{E}_{\bar{2}}$. Subtract Eqs. (24b) and (24c), change $\mathbf{u}_{2}-\mathbf{u}_{\bar{2}}$ to $\mathbf{u}_{d}$, and then insert expression (23b) for $\mathbf{E}_{2}-\mathbf{E}_{\bar{2}}$. Introducing a force constant due to LFEs, $K_{e}$, for the polar optical motion,

[TABLE]

and another force constant due to LFEs, $K_{2,d}$, for the nonpolar motion,

[TABLE]

and using the Born charge $e_{B}$ [expression (19)], the equations of motion reduce to

[TABLE]

Multiplying Eqs. (27a) and (27b) by $2m_{2}$ and $m_{1}$ respectively, subtracting and then dividing by $(m_{1}+2m_{2})$, one finds

[TABLE]

Changing $\mathbf{u}_{1}-\mathbf{u}_{c}$ to $\mathbf{w}$ with Eq. (15), one obtains

[TABLE]

where

[TABLE]

$\omega_{0}$ being the intrinsic oscillator frequency (i.e., without $\mathbf{E}$) for the polar optical vibrations.

From Eqs. (20) and (30b) one finds two equal $a$-coefficients

[TABLE]

Eq. (27c) is the equation of motion for the optical vibrations of the two S ions with reduced mass $m_{2}/2$; introducing $\mathbf{w}_{d}$,

[TABLE]

the equation becomes

[TABLE]

where

[TABLE]

$\omega_{2\bar{2}}$ being the frequency of the long-wavelength in-plane nonpolar modes, i.e., the LO1 and TO1 modes, which are nondispersive, and doubly degenerate at $\Gamma$.

The equation of motion (29) and polarization equation (16) describe the in-plane polar optical vibrations of ML MoS2, the vibrational properties of which will be derived in Sec. III below. As the lattice vibrations considered here are of long wavelengths, the pair of equations (29) and (16) constitute a macroscopic description of the lattice motion. $\mathbf{E}$ appearing in these equations represents the in-plane component of the macroscopic field in the thin ML and the $\mathbf{E}$ at $z=0$ is used (as was detailed above), i.e., $\mathbf{E}=\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$. It is evident from Eqs. (27a) and (27b) that the center of mass of the three ions in the cell remains stationary (frequency $\omega=0$), yielding trivial nondynamical solutions.

Now we consider out-of-plane motion with displacements $\mathbf{u}_{\nu}$ and dipole moments $\mathbf{p}_{\nu}$ parallel to $\mathbf{e}_{z}$. Let $\alpha_{\nu}^{\prime}$ be the electronic polarizability of the $\nu$ ions; usually $\alpha_{\nu}^{\prime}\neq\alpha_{\nu}$ ($\alpha_{\nu}$ is the in-plane polarizability), both being components of the polarizability tensor Born and Huang (1954). Let $K_{1}^{\prime}$ be the spring force constant between the Mo and S ions and $K_{2}^{\prime}$ be the force constant between the two S layers. Let $2e_{a}^{\prime}$ denote the static effective charge on the Mo atoms, which may differ from the in-plane effective charge $2e_{a}$ owing to the anisotropic 3D charge density distribution Zhang (2019); thus the charges on the ions are $e_{1}^{\prime}=2e_{a}^{\prime}$, $e_{2}^{\prime}=e_{\bar{2}}^{\prime}=-e_{a}^{\prime}$. Solving the Poisson equation one finds the electrostatic potential

[TABLE]

and then the macroscopic field,

[TABLE]

Evidently the in-plane component of the field is zero in all the three atomic layers. The local fields (i.e., Lorentz relations) at the Mo and two S sites are given by

[TABLE]

respectively, where the $Q^{\prime}$ coefficients are $Q_{0}^{\prime}=-2Q_{0}$, $Q_{1}^{\prime}=-2Q_{1}$ and $Q_{2}^{\prime}=-2Q_{2}$ from a symmetry analysis Mikhailov (2013). The macroscopic field has a $\bar{\delta}(z)$ term [Eq. (36a)] but it makes no contribution to the field change experienced by an ion $\nu$ owing to its own displacement $\mathbf{u}_{\nu}$, i.e., $4\pi\mathbf{u}_{\nu}\sum_{\nu^{\prime}}e_{\nu^{\prime}}^{\prime}\bar{\delta}(z)/s=4\pi\mathbf{u}_{\nu}(e_{1}^{\prime}+e_{2}^{\prime}+e_{\bar{2}}^{\prime})\bar{\delta}(z)/s=0$. Thus the field changes at the centers of Mo and S owing to their own displacements are

[TABLE]

Similarly, the macroscopic field has no contribution to the force exerted on the induced dipole $\bm{\mu}_{\nu}$ by all other ions again because of the charge neutrality in each cell. Repeating the process as before, one finds that with replacements $Q_{0}\rightarrow Q_{0}^{\prime}$, $Q_{1}\rightarrow Q_{1}^{\prime}$, $Q_{2}\rightarrow Q_{2}^{\prime}$, $e_{a}\rightarrow e_{a}^{\prime}$, $\alpha_{\nu}\rightarrow\alpha_{\nu}^{\prime}$, $K_{1}\rightarrow K_{1}^{\prime}$ and $K_{2}\rightarrow K_{2}^{\prime}$, all the equations above for in-plane motion are applicable to the out-of-plane motion. The Born charge is

[TABLE]

where

[TABLE]

The polarization equation is

[TABLE]

where $\mathbf{w}$ is the optical displacement [Eq. (15)] and $c_{22}$ is out-of-plane electronic susceptibility $\chi_{e}^{\prime}$,

[TABLE]

and both $c_{21}$ and $c_{22}$ are related to the Born charge,

[TABLE]

The equation of motion is

[TABLE]

where

[TABLE]

$\omega_{0}^{\prime}$ being the intrinsic oscillator frequency for the polar optical vibrations.

For the out-of-plane nonpolar optical vibrations the equation of motion is given by

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

and $\omega_{2\bar{2}}^{\prime}$ being the frequency of the long-wavelength nonpolar modes, i.e., the ZO1 modes.

In the lattice equations (44) and (41) $\mathbf{E}$ represents the field in the thin ML taking $\mathbf{E}(\bm{\rho},0)$, and evidently $\mathbf{E}(\bm{\rho},0)=\mathbf{E}_{z}(\bm{\rho},0)$ for the out-of-plane vibrations.

When the equations for the in-plane motion [Eqs. (29), (16) and (33)] and out-of-plane motion [Eqs. (44), (41) and (46)] are considered simultaneously, they can be rewritten for clarity as

[TABLE]

and

[TABLE]

respectively, where $a_{12}=a_{21}$ and $c_{12}=c_{21}$, and the linear coefficients have been given by the expressions above. The two pairs of equations involving macroscopic field have the same forms as those for ML hBN Zhang (2019) and are also similar to Huang’s equations for bulk crystals Huang (1951a, b); Born and Huang (1954), whereas the other two equations show the effective forces in a similar form to the Hooke’s law formula. By means of the averaging $\bar{\delta}(z)$, not only the 2D macroscopic polarization is introduced but also the 2D lattice equations are obtained with two pairs of equal linear coefficients.

The relation $a_{12}=a_{21}$ allows us to introduce an energy density (energy per unit area) $u_{p}$ associated with the in-plane optical vibrations, as a function of $\mathbf{w}_{\bm{\rho}}(\bm{\rho})$, $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ and $\mathbf{w}_{d,\bm{\rho}}(\bm{\rho})$,

[TABLE]

Inserting this $u_{h}$ expression into $\ddot{\mathbf{w}}_{\bm{\rho}}(\bm{\rho})=-\nabla_{\mathbf{w}}u_{h}$, $\bm{\mathcal{P}}(\bm{\rho})=-\nabla_{\mathbf{E}}u_{h}$ and $\ddot{\mathbf{w}}_{d,\bm{\rho}}(\bm{\rho})=-\nabla_{\mathbf{w}_{d}}u_{h}$, where $\nabla_{\mathbf{w}}$ denotes $(\frac{\partial}{\partial w_{x}},\frac{\partial}{\partial w_{y}})$, for instance, and the vector subscripts $\mathbf{w}$ and $\mathbf{E}$ and $\mathbf{w}_{d}$ of the del $\nabla$ are shorthand notations for in-plane vectors $\mathbf{w}_{\bm{\rho}}(\bm{\rho})$, $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ and $\mathbf{w}_{d,\bm{\rho}}(\bm{\rho})$, one rederives the lattice equations (50a), (50b) and (50c). Note that, similar to the bulk case Born and Huang (1954), it is not a priori obvious that an energy density of such a simple form should exist for ML MoS2, in particular when considering $\mathbf{E}$ is the macroscopic field – not simply an externally applied field.

With the relation $c_{12}=c_{21}$, similarly, an areal energy density associated with the out-of-plane optical vibrations can be defined,

[TABLE]

from which the lattice equations (51a), (51b) and (51c) can be rederived through $\ddot{w}_{z}=-\partial u_{v}/\partial w_{z}$, $\mathcal{P}_{z}=-\partial u_{v}/\partial E_{z}$ and $\ddot{w}_{d,z}=-\partial u_{v}/\partial w_{d,z}$ [$w_{z}=w_{z}(\bm{\rho})$, $E_{z}=E_{z}(\bm{\rho},0)$ and $w_{d,z}=w_{d,z}(\bm{\rho})$].

II.2 In-plane and out-of-plane polar optical modes

The nonpolar optical modes due to the relative motion of the negative ions are dispersionless [Eqs. (34) and (48)] so we focus on the polar modes, i.e., LO2, TO2 and ZO2 (but for simplicity the index subscripts 2 are omitted hereafter, causing no confusion), corresponding to the motion of the positive ions relative to the negative ions in the 2D polar crystals. The in-plane optical modes can be obtained from Eqs. (29) and (16) plus the equation of electrostatics $\nabla\cdot(\mathbf{E}+4\pi\mathbf{P})=0$, where $\mathbf{P}$ is the polarization $\mathbf{P}=\bm{\mathcal{P}}\bar{\delta}(z)$, and $\mathbf{E}$ is an electrostatic field $\mathbf{E}=-\nabla\phi$. Let $\mathbf{w}(\bm{\rho})=\mathbf{w}_{0}e^{i\mathbf{k}\cdot\bm{\rho}}$ and electric potential $\phi(\bm{\rho},z)=\varphi(z)e^{i\mathbf{k}\cdot\bm{\rho}}$ (time dependence $e^{-i\omega t}$ is omitted). Expressing $\mathbf{E}$ in terms of $\varphi$ one has the in-plane component of $\mathbf{E}$ at $z=0$ in the ML, $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)=-i\mathbf{k}\varphi(0)e^{i\mathbf{k}\cdot\bm{\rho}}$. Applying divergence $\nabla_{\bm{\rho}}$ to Eq. (16), one finds the polarization charge density through $-\nabla\cdot\mathbf{P}=-\bar{\delta}(z)\nabla_{\bm{\rho}}\cdot\bm{\mathcal{P}}$ and then has Poisson’s equation,

[TABLE]

Inserting the expansions of $\varphi(z)$ [Eq. (1)] and $\delta(z-z_{\nu})$ [Eq. (2)] into this equation, one obtains

[TABLE]

From the integral $\varphi(0)=\int_{-\infty}^{\infty}\hat{\varphi}(q)dq$, one finds $\varphi(0)$,

[TABLE]

and then obtains the electric field in the ML,

[TABLE]

where $\gamma_{k}=(1+2e^{-kd/2})/3$.

For a normal mode with wavevector $\mathbf{k}$, Eq. (29) becomes

[TABLE]

Expression (57) admits two possibilities for $\mathbf{w}\cdot\mathbf{k}$, (i) $\mathbf{w}\cdot\mathbf{k}=0$, or (ii) $\mathbf{w}\cdot\mathbf{k}\neq 0$. In case (i) $\mathbf{w}\perp\mathbf{k}$ so the normal modes are transverse waves. According to Eqs. (55) and (56), $\hat{\varphi}(q)=0$. Then from Eq. (1) one finds $\varphi(z)=0$, and thus the macroscopic field vanishes, $\mathbf{E}(\mathbf{r})=0$. Therefore the frequency of the TO mode is

[TABLE]

independent of wavevector; that is, the long-wavelength TO modes are dispersionless.

In case (ii) the electric field $\mathbf{E}_{\bm{\rho}}(\mathbf{r})$ is nonzero, and evidently the vectors $\mathbf{w}(\bm{\rho})$, $\mathbf{E}_{\bm{\rho}}(\mathbf{r})$, $\mathbf{P}(\mathbf{r})$ associated with the mode are all longitudinal, i.e., $\mathbf{w}(\bm{\rho})\parallel\mathbf{E}_{\bm{\rho}}(\mathbf{r})\parallel\mathbf{P}(\mathbf{r})\parallel\mathbf{k}$. Inserting expression (57) into Eq. (58) yields the frequency of the longitudinal optical (LO) mode,

[TABLE]

where the $a$-coefficients are changed to $e_{B}$, $\chi_{e}$, $\omega_{0}$ by Eqs. (18), (20) and (59). For the long wavelengths, $k\ll 1/a\approx 1/d$, $\gamma_{k}\approx 1$, one obtains the general LO phonon dispersion relation,

[TABLE]

independent of the specific form of $\bar{\delta}(z)$, as the effect of the finite thickness of ML MoS2 is negligible for the long waves. $\omega_{l}$ increases monotonically with $k$, with an upper bound $\omega_{M}$ at very large $k$, $\omega_{M}=\sqrt{\omega_{0}^{2}+e_{B}^{2}/(\bar{m}s\chi_{e})}$. The dispersion relation (61) is identical to the analytical expression of Sohier et al. Sohier et al. (2017); Zhang (2019).

Expanding $\omega_{l}(k)$ [Eq. (61)] to second order in $k$, one has

[TABLE]

where $c_{l}=\pi e_{B}^{2}/(\bar{m}s\omega_{0})$ Michel and Verberck (2009) is the norm of the LO phonon group velocity corresponding to the slope of the dispersion curve at $\Gamma$ Sohier et al. (2017). For an arbitrary $k$, the norm of the group velocity,

[TABLE]

decreases as $k$ becomes larger, showing that the dispersion curve becomes flatter as $k$ increases. The LO phonon dispersion curves from first-principles calculations Sohier et al. (2017) can be described very well using Eqs. (61) and (63).

The density of states of the LO modes can be obtained from the dispersion Eq. (61),

[TABLE]

Having $\varphi(0)$ [Eq. (56)], one substitutes expression (55) back into Eq. (1) to find the electric potential and then obtains the macroscopic field associated with the LO mode,

[TABLE]

consistent with Eqs. (3a) and (3b) above.

This macroscopic field results in a higher LO frequency than the TO frequency at a finite $k$ Sánchez-Portal and Hernández (2002) with the splitting given in Eq. (60). No splitting occurs in the limit $k\rightarrow 0$ as the macroscopic field vanishes [Eq. (65a)], which is different from the case in bulk 2H-MoS2, where LO-TO splitting occurs at the $\Gamma$ point Cai et al. (2014). Thus the transparent expressions (59) and (60) describe the degeneracy at $\Gamma$ and the splitting at a finite wavevector of the LO and TO modes, well-known phenomena of the 2D semiconductors Sánchez-Portal and Hernández (2002); Wirtz et al. (2003); Michel and Verberck (2009); Topsakal et al. (2009); Ataca et al. (2011); Cai et al. (2014); Sohier et al. (2017).

The out-of-plane optical modes can be obtained from Eqs. (44) and (41) as follows. Insert the field in the ML $\mathbf{E}_{z}(\bm{\rho},0)=-\mathbf{e}_{z}\varphi^{\prime}(0)e^{i\mathbf{k}\cdot\bm{\rho}}$ into Eq. (41) and then express the polarization charge density $-\nabla\cdot\mathbf{P}$ in terms of $\mathbf{w}\cdot\mathbf{e}_{z}$ and $\varphi^{\prime}(0)$. The Poisson equation is given by

[TABLE]

When the expansions of $\varphi(z)$ and $\delta(z-z_{\nu})$ [Eqs. (1) and (2)] are inserted, one finds $\hat{\varphi}(q)$,

[TABLE]

and after the integration $\varphi^{\prime}(z)=\int_{-\infty}^{\infty}iq\hat{\varphi}(q)e^{iqz}dq$, one obtains

[TABLE]

This expression has a $\bar{\delta}(z)$ term (divergent at $z_{\nu}=0,\pm d/2$), because the atomic layers are treated as geometric planes where the ionic charge distribution and polarization density $\mathbf{P}$ have a $\delta(z-z_{\nu})$ form. Now we approximate $\delta(z-z_{\nu})$ by a Gaussian distribution with a small effective thickness $\varepsilon$ ($\varepsilon\rightarrow 0$ ), $\delta_{\varepsilon}(z-z_{\nu})=\frac{1}{\sqrt{\pi}\varepsilon}e^{-(z-z_{\nu})^{2}/\varepsilon^{2}}$, as in previous theoretical study Michel and Verberck (2009) and first-principles calculations Sohier et al. (2017). Taking $z=0$ in Eq. (68) then one finds $\varphi^{\prime}(0)$,

[TABLE]

where

[TABLE]

The electric field in the ML $\mathbf{E}_{z}(\bm{\rho},0)$ follows,

[TABLE]

Substituting this field into the equation of motion (44), one obtains the frequency of the out-of-plane mode,

[TABLE]

When the broadening $\varepsilon$ is similar to the dimensions of a lattice cell, expression (70) is simplified to a constant $\zeta=4\sqrt{\pi}[1+2e^{-(d/2\varepsilon)^{2}}]/3$, and the phonon frequency becomes independent of wavevector; for ML MoS2 when $\varepsilon\ll 4\sqrt{\pi}\chi_{e}^{\prime}/3$ the frequency reduces to

[TABLE]

The EP of ions (i.e., $\chi_{e}^{\prime}\neq 0$ ) ensures a finite frequency, otherwise $\omega_{z}$ [Eq. (72)] becomes very large and even $\omega_{z}\rightarrow\infty$ [Eq. (73)] when $\chi_{e}^{\prime}$ is neglected in the rigid ion model (i.e., without EP).

With $\varphi^{\prime}(0)$ [Eq. (69)] then one can get $\varphi^{\prime}(z)$ from Eq. (68) and also obtain $\varphi(z)$ after substituting the $\hat{\varphi}(q)$ expression (67) into Eq. (1). It follows that the macroscopic field associated with the out-of-plane polar mode is given by the following expressions [consistent with Eqs. (36a) and (36b)],

[TABLE]

Note that the in-plane component, antisymmetric with respect to the Mo layer, is negligible at very small $k$ and $\varepsilon$.

II.3 2D lattice dielectric susceptibility and dielectric function

The in-plane lattice susceptibility can be derived from Eqs. (29) and (16) by considering periodic solutions $\mathbf{E}$, $\bm{\mathcal{P}}$, $\mathbf{W}$ $\propto e^{-i\omega t}$ due to an electric disturbance with frequency $\omega$. The susceptibility, $\chi=\bm{\mathcal{P}}(\bm{\rho})/\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$, is given by

[TABLE]

where the $\omega$-terms are due to lattice vibration. From Eq. (75) the static susceptibility is

[TABLE]

Using Eq. (76), one can express $a_{12}$ [Eq. (30b)] in terms of the 2D susceptibilities as

[TABLE]

which has the same form as its 3D counterpart $b_{12}$ ($b_{21}$) that is expressed in terms of the 3D susceptibilities or dielectric constants Born and Huang (1954).

Further the 2D dynamical susceptibility $\chi(\omega)$ can be conveniently expressed in terms of the oscillator frequency $\omega_{0}$ and the high-frequency and static susceptibilities $\chi_{e}$, $\chi_{0}$ as

[TABLE]

which is similar to its counterpart of bulk crystals Born and Huang (1954).

The 2D susceptibility for vertical polarization, $\chi^{\prime}=\bm{\mathcal{P}}(\bm{\rho})/\mathbf{E}_{z}(\bm{\rho},0)$, is derived from Eqs. (44) and (41),

[TABLE]

and the static susceptibility is

[TABLE]

Using this, one can rewrite $c_{12}$ [Eq. (45b)] in terms of the 2D susceptibilities as

[TABLE]

and also transform the dynamical susceptibility $\chi^{\prime}(\omega)$ as

[TABLE]

On eliminating $e_{B}^{\prime 2}/(\bar{m}s)$ in Eq. (73) with expression (80) one obtains a simple relation

[TABLE]

The lattice DF can be derived by introducing a test charge to calculate the total potential in the electrostatic approximation as in Ref.Cudazzo et al. (2011). To find the potential let the test charge density function $\sigma$ have the form $\sigma=\sigma_{0}e^{i\mathbf{k}\cdot\bm{\rho}}\delta(z-z_{a})$, where $z_{a}$ is the $z$-coordinate of the test charge not necessarily confined in the ML, and time dependence $e^{-i\omega t}$ is omitted for clearness. Now the equation of electrostatics is $\nabla\cdot(\mathbf{E}+4\pi\mathbf{P})=4\pi\sigma$, $\mathbf{E}$ being the total field of the test charge and polarization charge. In general the field of the test charge is nonzero in the ML, and thus the lattice responds generating in-plane [Eqs. (50a) and (50b)] and out-of-plane [Eqs. (51a) and (51b)] vibrations, with all the quantities varying as $e^{i(\mathbf{k}\cdot\bm{\rho}-\omega t)}$. As the induced potential associated with the out-of-plane polarization is zero at the centre of the ML and antisymmetric with respect to the Mo layer [refer to Eq. (35) above], the out-of-plane motion makes no contribution to the DF Zhang (2019) and we only need to consider the in-plane motion, taking the dielectric polarization $\mathbf{P}=\bm{\mathcal{P}}_{\bm{\rho}}\bar{\delta}(z)$, with $\bm{\mathcal{P}}_{\bm{\rho}}=\chi(\omega)\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$.

Writing $\mathbf{E}=-\nabla\phi$ with the total potential $\phi(\mathbf{r})=\varphi(z)e^{i\mathbf{k}\cdot\bm{\rho}}$, one has Poisson’s equation,

[TABLE]

where $\phi(\bm{\rho},0)$ is the total potential in the ML Cudazzo et al. (2011), $\phi(\bm{\rho},0)=\varphi(0)e^{i\mathbf{k}\cdot\bm{\rho}}$. Using expansions (1) and (2) then one finds $\hat{\varphi}(q)$,

[TABLE]

The integration $\varphi(0)=\int_{-\infty}^{\infty}\hat{\varphi}(q)dq$ gives $\varphi(0)$,

[TABLE]

Evidently the numerator multiplied by $e^{i\mathbf{k}\cdot\bm{\rho}}$ is the electric potential in the ML of the test charge, and therefore the DF of the 2D lattice is the denominator,

[TABLE]

The DF $\epsilon(k,\omega)$ is a longitudinal DF due only to the LO vibrations (i.e., no contribution from the TO or ZO vibrations), because the polarization charge density for the in-plane motion is given by $-\bar{\delta}(z)\nabla_{\bm{\rho}}\cdot\bm{\mathcal{P}}_{\bm{\rho}}(\bm{\rho})=-\bar{\delta}(z)[a_{21}-a_{22}(\omega^{2}+a_{11})/a_{12}]i\mathbf{k}\cdot\mathbf{w}_{\bm{\rho}}\neq 0$, and also $\mathbf{w}_{\bm{\rho}}\parallel\bm{\mathcal{P}}_{\bm{\rho}}\parallel\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)\parallel\mathbf{k}$. Without test charge ($\sigma=0$), a finite electric field $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ and potential $\phi(\bm{\rho},0)$ occur due to the LO vibrations, and from Eq. (86) therefore the LO modes are the solutions to

[TABLE]

which is very similar to the bulk crystal result Born and Huang (1954). Combining expressions (87) and (88) yields LO phonon dispersion $\omega_{l}(k)$ indeed identical with expression (60).

The static and high-frequency DFs are $\epsilon_{0}(k)=1+2\pi\chi_{0}k\gamma_{k}$, and $\epsilon_{\infty}(k)=1+2\pi\chi_{e}k\gamma_{k}$, respectively. The lattice DF of the 2D crystal can be transformed as

[TABLE]

i.e., in a very similar form to the lattice DF of bulk crystals Born and Huang (1954), the difference being that here the LO phonon frequency and DFs are both wavevector-dependent. Expression (89) shows that the TO phonon frequency $\omega_{t}$ is a pole of $\epsilon(k,\omega)$ whereas at LO phonon frequencies $\omega_{l}(k)$ the DF is zero. Further in bulk there is the Lyddane–Sachs–Teller (LST) relation Lyddane et al. (1941), $\omega_{l}^{2}/\omega_{t}^{2}=\epsilon_{0}/\epsilon_{\infty}$; for the 2D crystal a similar relation can be obtained from expression (89),

[TABLE]

The extended LST relation (90) shows a clear physical significance: it connects the phonon frequencies with the static and high-frequency DFs; given the former, then ratio of the latter is known, and vice versa. Dispersion occurs in the 2D crystal with the polar LO mode and DFs strongly dependent on the wavevector. Knowing the LO phonon dispersion $\omega_{l}(k)$ [$\omega_{l}(0)=\omega_{t}$], for instance, from an experimental measurement or theoretical calculation, then the ratio between the two DFs $\epsilon_{0}(k)$ and $\epsilon_{\infty}(k)$ is obtained. Meanwhile, as the LO and TO mode frequencies differ due solely to the macroscopic field, the relation is very useful for evaluating the phonon frequency increase brought by the field. Similarly, for the out-of-plane motion the frequency–susceptibility relation (83) measures the effect of the macroscopic field on the phonon frequency.

For the long wavelengths, taking $\gamma_{k}=1$, Eq. (87) becomes a general DF expression independent of the specific form of $\bar{\delta}(z)$; it is also the same as the DF $\epsilon(k,\omega)$ of ML hBN Zhang (2019), and accordingly the $\epsilon_{0}(k)$ and $\epsilon_{\infty}(k)$ expressions are simplified to the DF formula for a normal 2D dielectric deduced by Cudazzo et al. Cudazzo et al. (2011).

III Numerical results and further discussions

III.1 Various lattice dynamical properties of ML TMDs

In what follows we confine ourselves to in-plane motion, and specifically the polar LO and TO modes, corresponding to the vibrations of the positive ions relative to the negative ions. These modes at the zone center are denoted as $E^{\prime}$; for ML MoS2, for instance, the first-principles calculations by density-functional perturbation theory (DFPT) find $\omega(E^{\prime})$=380.2 cm*-1* in Ref.Ataca et al. (2011) and $\omega(E^{\prime})$=402.7 cm*-1* in Ref.Cai et al. (2014). In a recent study Zhang (2019) we found that the atomic polarizability of the unit cell of ML hBN is significantly reduced compared to the total free-atom polarizablity of B and N. The polarizabilities $\alpha_{1}$ and $\alpha_{2}$ of the constituent atoms in ML TMDs may differ significantly from the free-atom values and are unknown quantities so the values of the $a$-coefficients can not be determined with expressions such as Eqs. (17a) and (17b). A set of three values such as $\omega_{0}$, $e_{B}$, $\chi_{e}$ can be obtained from first-principles calculations and therefore the three mutually independent $a$-coefficients $a_{11}$, $a_{12}$ (or $a_{21}$), $a_{22}$ of the lattice equations can be determined through the expressions (18), (30a) and (30b). Furthermore, as the three expressions (21), (59) and (76) relate the three quantities of a ML TMD, namely the two macroscopic susceptibilities $\chi_{e}$, $\chi_{0}$ and the collective vibration frequency $\omega_{0}$ to the three microscopic quantities, i.e., the static effective charge $e_{a}$, the Born charge $e_{B}$ and the spring force constant $K_{1}$, one can calculate $e_{a}$, $K_{1}$, $\chi_{0}$ using these expressions. Of the four microscopic quantities $\alpha_{1}$, $\alpha_{2}$, $e_{a}$, $K_{1}$ of our dipole lattice model on which the $a$-coefficients originally depend [refer to Eqs. (17a), (17b), (19), (25) and (30a)], only the two atomic polarizabilities $\alpha_{1}$, $\alpha_{2}$ are unknown. The adoption of the set of known quantities $\omega_{0}$, $e_{B}$, $\chi_{e}$ calculated from first principles facilitates the use of the deduced equations by circumventing the unknowns $\alpha_{1}$ and $\alpha_{2}$, as we shall see below.

The atomic polarizabilities $\alpha_{1}$ and $\alpha_{2}$ can not be determined from Eqs. (17b) and (19) even when $e_{B}$ and $\chi_{e}$ are known, because $e_{B}$ and $\chi_{e}$ are related via Eq. (21) and not independent of each other. Interestingly, the sum $\alpha=\alpha_{1}+2\alpha_{2}$, i.e., the total atomic polarizability of the unit cell is limited in an interval (given $\chi_{e}$), which can be determined using the Clausius-Mossotti relation (17b). Changing variable $\alpha_{2}$ to $\alpha_{1}$ for a $\alpha$ transforms Eq. (17b) into a quadratic equation in variable $\alpha_{1}$. Then follow the lines of the derivation given in Ref.Zhang (2019) (Appendix A). By requiring the discriminant $\Delta_{\alpha}\geq 0$ as well as $\alpha_{1}\geq 0$ and $\alpha_{2}\geq 0$, one obtains the interval for $\alpha$,

[TABLE]

In a recent study Sohier et al. (2017) the LO phonon dispersion of a ML TMD or hBN is given by $\omega_{l}^{2}=\omega_{0}^{2}+\mathcal{S}k/(1+r_{eff}k)$, where the $\mathcal{S}$ parameter relates to the Born charge via $\mathcal{S}=2\pi e_{B}^{2}/(\bar{m}s)$ Sohier et al. (2017); Zhang (2019), and $r_{eff}$ is an effective screening length, given by by $r_{eff}=\epsilon_{p}t/2$ with an effective medium model Sohier et al. (2016, 2017), $\epsilon_{p}$ and $t$ being effective dielectric constant and effective thickness of the ML material. Both parameters $\mathcal{S}$ and $r_{eff}$ are computed by first-principles calculation Sohier et al. (2017). A comparison with Eq. (61) gives the screening length $r_{eff}=2\pi\chi_{e}$. Using the $\mathcal{S}$ and $r_{eff}$ values (listed in Table 1 of Ref.Sohier et al. (2017) together with $\omega_{0}$) for a ML TMD or hBN, one finds $\lvert e_{B}\rvert$ and $\chi_{e}$, and adding $\omega_{0}$, one has the values for the set of three quantities ($\lvert e_{B}\rvert$, $\chi_{e}$, $\omega_{0}$), and further calculates the quantities $\chi_{0}$, $\lvert e_{a}\rvert$, $K_{1}$ together with the interval of $\alpha$, as detailed above. The result is presented in Table 1 for five ML TMDs and ML hBN for comparison.

In Table 1 we give the absolute values of the Born charges as their signs cannot be determined from the $\mathcal{S}$ values of Ref.Sohier et al. (2017). Note that the effective charges $e_{B}$ in ML TMDs have been calculated in first-principles studies, but both positive ($e_{B}>0$) Ataca et al. (2011) and negative ($e_{B}<0$) Michel et al. (2017); Danovich et al. (2017) values have been found. Considering that the key quantities such as $\omega_{l}(k)$ [Eq. (61)] and $\chi_{0}$ [Eq. (76)] depend on $e_{B}^{2}$, and $e_{B}$ pairs with $e_{a}$ in expressions (21) and (59) (recall the pair have the same sign, i.e., $e_{B}e_{a}=\lvert e_{B}e_{a}\rvert$), with the absolute values we still find some interesting results. First, Comparing the $\chi_{e}$ and $\chi_{0}$ values of the ML TMDs, we see that the susceptibility is only slightly increased when the lattice vibrational contribution is included: there is a 0.3% increase for ML WS2, 9% increase for ML MoTe2 and 1-3% increase for the other ML TMDs, and these increases are smaller than the 37% susceptibility increase of ML hBN. The 2D susceptibility is a key parameter of the effective 2D interaction in ML TMDs, i.e., the Keldysh potential Keldysh (1978); Cudazzo et al. (2011), for which either $\chi_{e}$ or $\chi_{0}$ can be used, and as the above result indicates, they will not cause much difference to the evaluation of such 2D interaction. Second, the Born charge values are quite small in the ML TMDs, and they are all smaller than 2$e$ and the $e_{B}$ value of ML hBN except for ML MoTe2 which has a Born charge of 3.1$e$, thus resulting in a small vibrational contribution to the susceptibility. All these $\lvert e_{B}\rvert$ values are very close to recent density functional theory (DFT) Danovich et al. (2017) and DFPT Michel et al. (2017) calculations (compare to Table I of Ref.Danovich et al. (2017) and Table II of Ref.Michel et al. (2017)) and are within a 7% of deviation except for ML WS2 with a 9% deviation in Born charge. Compared to 2D hBN, in ML TMDs the ions also carry a small static effective charge. For ML MoS2, for instance, the charge on the Mo ions is 0.22$e$, less than half the charge on the B ions of ML hBN. No $e_{a}$ values have been reported for ML TMDs so far but in bulk 2H-MoS2 first-principles calculation in conjunction with Mulliken analysis yields a charge transfer of 0.43 electrons from one Mo atom to two S atoms Ataca et al. (2011). Third, in a ML TMD the atomic polarizability of the unit cell $\alpha$ has a limited range of values calculated with inequality (91), and similar to ML hBN these values are much smaller than the values of the total free-atom polarizability $\alpha_{f}$ (ninth column), $\alpha_{f}=\alpha_{1,f}+2\alpha_{2,f}$ for ML TMDs and $\alpha_{f}=\alpha_{1,f}+\alpha_{2,f}$ for ML hBN, where the free-atom polarizabilities of the constituent atoms $\alpha_{1,f}$ and $\alpha_{2,f}$ are taken from Ref. Schwerdtfeger (2006). For ML MoS2, for instance, the upper and lower bounds of $\alpha$ are 43% and 72% smaller than the free-atom polarizability $\alpha_{f}$, respectively. The numerical result also shows that for each ML TMD the lower bound value of $\alpha$ is about half of its upper bound value. In addition, $\alpha$’s two bounds, $\alpha_{f}$ and $\chi_{e}$ vary with the ML material in quite a similar manner–the atoms that have a larger atomic polarizability $\alpha_{f}$ form a ML TMD with a larger dielectric susceptibility $\chi_{e}$. Fourth, a smaller spring force constant $K_{1}$ combined with a larger reduced mass $\bar{m}$ makes the ML TMDs have a lower mode frequency $\omega_{0}$ than 2D hBN. Five, as the group velocity $c_{l}$ corresponds to the slope of the LO phonon dispersion curve at $\Gamma$, a small $c_{l}$ for ML WS2 means its flat phonon dispersion while a large $c_{l}$ for ML MoTe2 corresponds to a steep slope of the dispersion near $\Gamma$.

III.2 Local field and polarizable ion effects in ML TMDs

Knowing the values of the microscopic quantities now we use them to evaluate the EP of ions and LFEs on the polar vibrations. Neglecting EP of ions ($\alpha_{1}=\alpha_{2}=0$), there is no electronic susceptibility, i.e., $\chi_{e}=0$, and the static susceptibility is also reduced substantially with only a small contribution due to lattice vibration. When LFEs are not accounted for (equivalent to $Q_{0}=Q_{1}=Q_{2}=0$), the high-frequency dielectric susceptibility of ML TMDs becomes $(\alpha_{1}+2\alpha_{2})/s=\alpha/s$ [refer to Eqs. (13), (17b) and (18)], the intervals of which are shown in Table 1 (last column). Clearly these are small values compared to the $\chi_{e}$ including the LFEs (second column), with the upper bounds being only 13-18% of $\chi_{e}$ and the lower bounds only 6-9% of $\chi_{e}$. Without LFEs, the Born charge becomes smaller, equal to the static effective charge, i.e., $e_{B}=2e_{a}$ [refer to Eqs. (13) and (19)], which is the same as in the case when ionic EP is neglected because $Q_{i}$ ($i=0,1,2$) and $\alpha_{j}$ ($j=1,2$) form products in the $e_{B}$’s expression (19). The LO phonon group velocity $c_{l}$ is also decreased substantially, which is only 3-6% of the $c_{l}$ value calculated with the polarizable ion model (PIM) including the LFEs. The intrinsic oscillator frequency $\omega_{0}$ is determined by $K_{1}-K_{e}$ [Eq. (30a)] rather than by $K_{1}$ alone, i.e., $\omega_{0}^{2}=2(K_{1}-e_{a}e_{B}Q_{1})/\bar{m}$, and thus excluding LFEs raises the $\omega_{0}$ value with a percentage increase given by $1/(K_{1}/K_{e}-1)$. We find that the percentage increase in $\omega_{0}^{2}$ of ML TMDs is between 0.2% for ML WS2 and 7.5% for ML MoTe2–much smaller than the 31% increase for ML hBN–indicating that the LFEs have only a small influence on the intrinsic frequency $\omega_{0}$ of ML TMDs.

Further the polarizable ion and local field effects on LO phonon dispersion can be evaluated using the values of Table 1. In the PIM including the LFEs, the LO phonon frequency $\omega_{l}(k)$ can be calculated from Eq. (61) using $e_{B}$, $\chi_{e}$ in Table 1 and $\omega_{TO}$ in Table 1 of Ref.Sohier et al. (2017), or equally by the dispersion formula of Sohier et al. Sohier et al. (2017) using $\mathcal{S}$, $r_{eff}$, $\omega_{TO}$ in their Table 1. When LFEs are neglected in the PIM, $\omega_{l}(k)$ becomes dependent on the unit-cell atomic polarizability $\alpha$, given by $\omega_{l}^{2}(k)=2K_{1}/\bar{m}+8\pi e_{a}^{2}k/[\bar{m}(s+2\pi\alpha k)]$, which was calculated using two values of $\alpha$, i.e., its lower bound $\alpha_{l}$ and upper bound $\alpha_{u}$ in Table 1. In the rigid ion model (RIM) accounting for the LFEs, the LO phonon dispersion is given by $\omega_{l}^{2}(k)=2(K_{1}-2e_{a}^{2}Q_{1})/\bar{m}+8\pi e_{a}^{2}k/(\bar{m}s)$, and when neglecting LFEs the LO phonon dispersion is simplified to $\omega_{l}^{2}(k)=2K_{1}/\bar{m}+8\pi e_{a}^{2}k/(\bar{m}s)$. A comparison of these phonon dispersion relations is demonstrated with ML MoS2 and ML MoTe2 in Figs. 2(a) and 2(b), respectively. We see that the phonon dispersion curves of the two ML TMDs are similar. The RIM without LFEs yields the highest phonon frequency, which increases linearly with the wavevector in the long wavelength region (top line). The LFEs cause a reduction to the frequencies but the dispersion remains linear (lower dashed line). When the EP of ions is included, the dispersion becomes nonlinear with a smaller slope at a larger wavevector in both cases of excluding (dotted and dot-dashed lines) and including the LFEs (solid line). Although the values of $\alpha_{1}$ and $\alpha_{2}$ are unknown, using the intervals of the total $\alpha$ the PIM without LFEs yields phonon frequencies that are limited in a very narrow range (i.e., dotted and dot-dashed lines), 0.5 cm*-1* for ML MoS2 and 2 cm*-1* for ML MoTe2. Clearly the phonon frequencies are reduced significantly due to the ionic EP and LFEs (solid line), and in particular, there is a steep slope on the small wavevector side corresponding to a substantially increased phonon group velocity. The phonon dispersion curves of both monolayers become flatter at larger wavevectors [refer to Eq. (63)]. A wider frequency range of ML MoTe2 as shown in Fig. 2, however means that it has a larger LO phonon group velocity $c_{l}$ than ML MoS2, consistent with the result in Table 1.

For out-of-plane optical vibrations, the polar mode ($A_{2}^{\prime\prime}$ mode) and nonpolar mode ($A_{1}$ mode) frequencies have been calculated for ML MoS2 using DFPT, $\omega(A_{2}^{\prime\prime})$=490.5 cm*-1* and $\omega(A_{1})$=423.9 cm*-1* in Ref.Cai et al. (2014); $\omega(A_{2}^{\prime\prime})$=465.0 cm*-1* and $\omega(A_{1})$=406.1 cm*-1* in Ref.Ataca et al. (2011). So far there has been no report on Born charge or 2D susceptibility, and when these basic quantities are available, the dynamical properties can be calculated using our expressions in a similar way to the in-plane vibration calculation.

IV Summary and Conclusions

We have deduced two sets of three equations [Eqs. (50a), (50b), (50c) and Eqs. (51a), (51b), (51c)] to describe the long wavelength in-plane and out-of-plane lattice vibrations, respectively, in monolayer TMDs using a microscopic dipole lattice model accounting for the LFEs and EP self-consistently. The two pairs of equations [Eqs. (50a), (50b) and Eqs. (51a), (51b)], which have the same forms as those for ML hBN Zhang (2019) and are also similar to Huang’s equations for bulk crystals, describe the polar optical lattice vibrations, whereas the other two equations [Eqs. (50c) and (51c)], in which the effective forces are proportional to the relative ionic displacements and have a similar form to the Hooke’s law formula, describe the nonpolar optical vibrations in the 2D polar crystals. The averaging of the microscopic distributions in the $z$ direction of the dipole moment density makes it possible to define a 2D macroscopic polarization $\bm{\mathcal{P}}(\bm{\rho})$ and express the volume polarization of the ML TMDs as a product of $\bm{\mathcal{P}}(\bm{\rho})$ and the averaged microscopic distribution $\bar{\delta}(z)$. These together with the 2D Lorentz relations are fundamental to deducing the lattice equations from the atomic theory. As the linear coefficients of the equations are related to the quantities that can be obtained from first principles calculation, the lattice equations are very useful for studying the lattice dynamical properties analytically. We have also obtained the expressions for the areal energy densities associated with the in-plane and out-of-plane optical vibrations, respectively.

The lattice equations [Eqs. (50a), (50b) or Eqs. (51a), (51b)] are solved simultaneously with the equation of electrostatics to deduce the polar optical modes. Explicit expressions have been obtained for the frequencies, the macroscopic fields associated with the polar LO and ZO modes, and also the LO phonon group velocity and density of states. The LO phonon dispersion relation, applicable to both ML TMDS and ML hBN, is identical to the analytical expression of Sohier et al., and it evidently shows that the LO and TO modes are degenerate at $\Gamma$ and split up at finite wavevectors due to the long-range macroscopic field, characteristic of these 2D polar crystals. The transparent LO phonon dispersion relation also describes the first-principles calculation results very well: the LO mode frequency increases with wavevector but the dispersion becomes flatter at larger wavevectors. The phonon frequency expressions show that apart from the LO branch all the other five optical branches are dispersionless at the long wavelengths. It is also found that the ZO phonon frequency is finite due to ionic EP, but otherwise with the rigid ion model the ZO frequency becomes extremely large.

The in-plane and out-of-plane lattice dielectric susceptibilities have been deduced from the lattice equations. The 2D longitudinal lattice dielectric function $\epsilon(k,\omega)$ has also been derived, allowing the LO phonon dispersion to be rederived from $\epsilon(k,\omega)=0$, similar to the case of bulk crystals. Further, the 2D Lyddane–Sachs–Teller relation (90) and frequency–susceptibility relation (83) have been derived for in-plane and out-of-plane vibrations, respectively, which are very useful for evaluating the effects of the macroscopic field on the phonon frequencies.

We have demonstrated an application of the analytical expressions by using them to study the lattice dynamical properties for in-plane vibrations, when knowing a set of three quantities $\omega_{0}$, $\chi_{e}$, $e_{B}$ from first-principles calculation. In ML TMDs except ML MoTe2 and MoSe2 the effective charges are quite small, causing a small vibrational contribution (less than 1.4%) to the dielectric susceptibility. Although the individual atomic polarizabilities in ML TMDs are unknown, the total atomic polarizability of the unit cell is limited in an interval which is obtained from the Clausius-Mossotti relation, and the intervals for ML TMDs and hBN have been calculated and also compared with the total free-atom polarizabilities. The LFEs and EP should be included simultaneously; otherwise, neglecting either or both underestimates substantially the calculated properties: the Born charge decreases simply to the value of the static effective charge (which becomes 78% smaller for ML MoS2, for instance), the dielectric susceptibilities $\chi_{e}$ and $\chi_{0}$ are underestimated by over 80%, and the LO phonon group velocity is underestimated by over 90%. With no EP or LFE, the LO modes have very small dispersion, showing linear dispersion for the two RIMs and a very narrow frequency range (0.5-2.0 eV) for the PIM without LFE, which is distinct from the LO phonon dispersion calculated after including both ionic EP and LFEs.

Acknowledgements.

We acknowledge support from the New Energy and Materials Collaboration project of the School of Physics, and the Natural Science Research Funds (Nos. 419080500175 & 419080500260) of Jilin University.