Long optical lattice vibrations and dielectric function of monolayer hexagonal boron nitride

TL;DR

This paper develops a microscopic lattice model to analyze optical vibrations and dielectric properties of monolayer hexagonal boron nitride, incorporating electronic polarization and local field effects for accurate predictions.

Contribution

It introduces a self-consistent lattice equation framework that explicitly accounts for electronic polarization and local field effects in 2D hexagonal BN.

Findings

Derived explicit phonon dispersion relations matching previous numerical results.

Established a 2D dielectric function based solely on LO phonon vibrations.

Highlighted the importance of including both ionic EP and LFEs for accurate dynamical properties.

Abstract

Using two pairs of lattice equations deduced from a microscopic dipole lattice model taking into account electronic polarization (EP) of ions and local field effects (LFEs) self-consistently, in-plane and out-of-plane optical vibrations in two-dimensional (2D) hexagonal BN are studied theoretically. The three mutually independent coefficients of either pair of lattice equations are determined by a set of three generally known quantities such as the 2D electronic and static susceptibilities and phonon frequency, making the lattice equations very useful for calculating the lattice dynamical properties. Explicit expressions are obtained for lattice vibrational energy density, and phonon dispersion, group velocity and density of states. The transparent phonon dispersion relations describe the previous numerical calculations very well, and the longitudinal optical (LO) phonon dispersion relation is identical to the analytical expression of Sohier {\it et al}. [Nano Lett. {\bf 17}, 3758 (2017)]. The out-of-plane phonon frequency is finite owing to ionic EP. A 2D lattice dielectric function $\epsilon(k,\omega)$ is derived--due solely to the LO vibrations--which also allows the LO phonon dispersion to be rederived simply from $\epsilon(k,\omega)=0$, similar to the bulk case. A 2D Lyddane--Sachs--Teller relation and a frequency--susceptibility relation are obtained for the in-plane and out-of-plane vibrations, respectively, connecting the phonon frequencies to the 2D dielectric functions or susceptibilities. Using three first-principles calculated parameters, the lattice dynamical properties are studied comprehensively, particular attention being paid to the EP and LFEs. The ionic EP and LFEs should be included simultaneously, but otherwise neglecting either or both causes large discrepancies to the calculated dynamical properties.

Tables

Table 1. Table 1: Two groups of calculated lattice-dynamical quantities associated with the in-plane optical vibrations in ML hBN, namely, the intrinsic oscillator frequency ω 0 subscript 𝜔 0 \omega_{0} , the 2D electronic susceptibility χ e subscript 𝜒 𝑒 \chi_{e} and static susceptibility χ 0 subscript 𝜒 0 \chi_{0} , the Born charge e B subscript 𝑒 𝐵 e_{B} , the static effective charge e a subscript 𝑒 𝑎 e_{a} , the effective spring force constant K 𝐾 K , the force constant due to LFEs K e subscript 𝐾 𝑒 K_{e} , and the LO phonon group velocity c l subscript 𝑐 𝑙 c_{l} . In the upper (lower) row, the set of three quantities ω 0 subscript 𝜔 0 \omega_{0} , χ e subscript 𝜒 𝑒 \chi_{e} , χ 0 subscript 𝜒 0 \chi_{0} ( ω 0 subscript 𝜔 0 \omega_{0} , χ e subscript 𝜒 𝑒 \chi_{e} , e B subscript 𝑒 𝐵 e_{B} ) are obtained from first-principles calculations of Refs. Ferrabone  et al. ( 2011 ); Erba  et al. ( 2013 ) (Ref. Sohier  et al. ( 2017 ) ), and then used to calculate the other quantities in this table with our expressions accounting for both EP and LFEs (see text).

${\omega }_0$ (cm-1)	${\chi }_e$ ($\AA$)	${\chi }_0$ ($\AA$)	$e_B$ ($e$)	$e_a$ ($e$)	$K$ (eV/${\AA }^2$)	$K_e$ (eV/${\AA }^2$)	$c_l$ (km/s)
1371a	0.85	1.31	2.70	0.61	59.796	17.613	37.24
1387b	1.22	1.67	2.71	0.46	56.495	13.317	37.10

Table 2. Table 2: Expressions for 2D electronic susceptibility χ e subscript 𝜒 𝑒 \chi_{e} , 2D static susceptibility χ 0 subscript 𝜒 0 \chi_{0} , Born charge e B subscript 𝑒 𝐵 e_{B} , TO phonon frequency ω t subscript 𝜔 𝑡 \omega_{t} , LO phonon dispersion ω l ​ ( k ) subscript 𝜔 𝑙 𝑘 \omega_{l}(k) , and LO phonon group velocity c l subscript 𝑐 𝑙 c_{l} associated with the in-plane optical vibrations of ML hBN in the RIM and PIM, obtained with or without taking LFEs into account (see text). Below the expressions are the values of these physical quantities and also the specific quantitative dispersion ω l ​ ( k ) subscript 𝜔 𝑙 𝑘 \omega_{l}(k) , obtained with a set of first-principles values χ e = 0.85 ​ Å subscript 𝜒 𝑒 0.85 italic-Å \chi_{e}=0.85~{}\AA , χ 0 = 1.31 ​ Å subscript 𝜒 0 1.31 italic-Å \chi_{0}=1.31~{}\AA , ω 0 subscript 𝜔 0 \omega_{0} =1371 cm -1 (put in the last row); values of χ e subscript 𝜒 𝑒 \chi_{e} and χ 0 subscript 𝜒 0 \chi_{0} are in Å italic-Å \AA , e B subscript 𝑒 𝐵 e_{B} in e 𝑒 e , ω t subscript 𝜔 𝑡 \omega_{t} and ω l subscript 𝜔 𝑙 \omega_{l} in cm -1 and c l subscript 𝑐 𝑙 c_{l} in km/s.

model		${\chi }_e$	${\chi }_0$	$e_B$	${\omega }_t$	${\omega }_l\left(k\right)$	$c_l$
RIM	no LFEs	0	$\frac{e_a^2}{sK}$	$e_a$	$\sqrt{\frac{K}{\overline{m}}}$	$\sqrt{\frac{K}{\overline{m}}+\frac{2\pi e_a^{2}k}{\overline{m}s}}$	$\frac{\pi e_a^2}{s\sqrt{\overline{m}K}}$
		0	0.017	0.61	1632	$\sqrt{{1632.4}^2+{837.4}^2\tilde{k}}$	1.61
	LFEs	0	$\frac{e_a^2}{s\left(K-e_a^2{Q}_1\right)}$	$e_a$	$\sqrt{\frac{K-e_a^2{Q}_1}{\overline{m}}}$	$\sqrt{\frac{K-e_a^2{Q}_1}{\overline{m}}+\frac{2\pi e_a^{2}k}{\overline{m}s}}$	$\frac{\pi e_a^2}{s\sqrt{\overline{m}\left(K-e_a^2{Q}_1\right)}}$
		0	0.018	0.61	1577	$\sqrt{{1577}^2+{837.4}^2\tilde{k}}$	1.67
PIM	no LFEs	$\frac{{\alpha }_1+{\alpha }_2}{s}$	$\frac{{\alpha }_1+{\alpha }_2}{s}+\frac{e_a^2}{sK}$	$e_a$	$\sqrt{\frac{K}{\overline{m}}}$	$\sqrt{\frac{K}{\overline{m}}+\frac{2\pi e_a^{2}k}{\overline{m}\left[s+2\pi \left({\alpha }_1+{\alpha }_2\right)k\right]}}$	$\frac{\pi e_a^2}{s\sqrt{\overline{m}K}}$
		[0.24,0.32]†	[0.26,0.34]†	0.61	1632	${\sqrt{{1632.4}^2+\frac{{837.4}^2\tilde{k}}{1+\gamma \tilde{k}}}}^{\oplus }$	1.61
	LFEs	Eq. (19)	Eq. (74)		Eq. (47)	Eq. (49)	Eq. (67)
		0.85	1.31	2.70	1371	$\sqrt{{1371}^2+\frac{{3691.5}^2\tilde{k}}{1+13.42\tilde{k}}}$	37.24

Table 3. Table 3: Comparison of intrinsic oscillator frequency ω 0 ′ superscript subscript 𝜔 0 ′ \omega_{0}^{\prime} , static effective charge e a ′ superscript subscript 𝑒 𝑎 ′ e_{a}^{\prime} , Born charge e B ′ superscript subscript 𝑒 𝐵 ′ e_{B}^{\prime} , force constant due to LFEs K e ′ superscript subscript 𝐾 𝑒 ′ K_{e}^{\prime} and effective force constant K ′ superscript 𝐾 ′ K^{\prime} associated with the out-of-plane optical modes in ML hBN calculated with the PIM including LFEs, for two values of frequency ω z subscript 𝜔 𝑧 \omega_{z} taken from previous first-principles calculations (Refs. Erba  et al. ( 2013 ); Topsakal  et al. ( 2009 ); Miyamoto  et al. ( 1995 ); Wirtz  et al. ( 2003 ) ), one experimental ω z subscript 𝜔 𝑧 \omega_{z} value (Ref. Rokuta  et al. ( 1997 ) ) and one low frequency ω z subscript 𝜔 𝑧 \omega_{z} of 405 cm -1 , using 2D electronic susceptibility χ e ′ = 0.151 superscript subscript 𝜒 𝑒 ′ 0.151 \chi_{e}^{\prime}=0.151 Å italic-Å \AA and static susceptibility χ 0 ′ = 0.164 superscript subscript 𝜒 0 ′ 0.164 \chi_{0}^{\prime}=0.164 Å italic-Å \AA calculated from first principles (Ref. Ferrabone  et al. ( 2011 ) ).

${\omega }_z$ (cm-1)	${\omega }_0^{\prime }$ (cm-1)	$e_a^{\prime }$ ($e$)	$e_B^{\prime }$ ($e$)	$K_e^{\prime }$ (eV/${\AA }^2$)	$K^{\prime }$ (eV/${\AA }^2$)
836a	802	1.26	-0.27	7.136	21.575
800b	768	1.20	-0.25	6.535	19.758
734c	704	1.11	-0.23	5.501	16.630
405d	389	0.61	-0.13	1.675	5.063

Table 4. Table 4: Expressions for 2D electronic susceptibility χ e ′ superscript subscript 𝜒 𝑒 ′ \chi_{e}^{\prime} , 2D static susceptibility χ 0 ′ superscript subscript 𝜒 0 ′ \chi_{0}^{\prime} , Born charge e B ′ superscript subscript 𝑒 𝐵 ′ e_{B}^{\prime} , intrinsic oscillator frequency ω 0 ′ superscript subscript 𝜔 0 ′ \omega_{0}^{\prime} and phonon frequency ω z subscript 𝜔 𝑧 \omega_{z} associated with the out-of-plane optical vibrations of ML hBN in the RIM and PIM, obtained with or without taking LFEs into account (see text). Below the expressions are the values of these physical quantities obtained with a set of first-principles values, χ e ′ = 0.151 superscript subscript 𝜒 𝑒 ′ 0.151 \chi_{e}^{\prime}=0.151 Å italic-Å \AA , χ 0 ′ = 0.164 superscript subscript 𝜒 0 ′ 0.164 \chi_{0}^{\prime}=0.164 Å italic-Å \AA , ω z = 836 subscript 𝜔 𝑧 836 \omega_{z}=836 cm -1 (given in the last row); values of χ e ′ superscript subscript 𝜒 𝑒 ′ \chi_{e}^{\prime} and χ 0 ′ superscript subscript 𝜒 0 ′ \chi_{0}^{\prime} are in Å italic-Å \AA , e B ′ superscript subscript 𝑒 𝐵 ′ e_{B}^{\prime} in e 𝑒 e , and ω 0 ′ superscript subscript 𝜔 0 ′ \omega_{0}^{\prime} and ω z subscript 𝜔 𝑧 \omega_{z} in cm -1 .

model		${\chi }_e^{\prime }$	${\chi }_0^{\prime }$	$e_B^{\prime }$	${\omega }_0^{\prime }$	${\omega }_z$
RIM	no LFEs	0	$\frac{e_a^{\prime 2}}{s{K}^{\prime }}$	$e_a^{\prime }$	$\sqrt{\frac{K^{\prime }}{\overline{m}}}$	$\infty$
		0	0.196	1.26	981	$\infty$
	LFEs	0	$\frac{e_a^{\prime 2}}{s\left(K^{\prime }+2e_a^{\prime 2}{Q}_1\right)}$	$e_a^{\prime }$	$\sqrt{\frac{K^{\prime }+2e_a^{\prime 2}{Q}_1}{\overline{m}}}$	$\infty$
		0	0.076	1.26	1571	$\infty$
PIM	no LFEs	$\frac{{\alpha }_1^{\prime }+{\alpha }_2^{\prime }}{s}$	$\frac{{\alpha }_1^{\prime }+{\alpha }_2^{\prime }}{s}+\frac{e_a^{\prime 2}}{s{K}^{\prime }}$	$e_a^{\prime }$	$\sqrt{\frac{K^{\prime }}{\overline{m}}}$	$\sqrt{\frac{K^{\prime }}{\overline{m}}+\frac{e_a^{\prime 2}}{\overline{m}\left({\alpha }_1^{\prime }+{\alpha }_2^{\prime }\right)}}$
		[0.36,0.48]†	[0.55,0.67]†	1.26	981	[1165,1220]†
	LFEs	Eq. (36b)	Eq. (78)		Eq. (38a)	$\sqrt{\frac{K^{\prime }+2e_a^{\prime }{e}_B^{\prime }{Q}_1}{\overline{m}}+\frac{e_B^{\prime 2}}{\overline{m}s{\chi }_e^{\prime }}}$
		0.151	0.164	-0.27	802	836

Equations

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

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

$${\mathbf{E}}_{\mathbf{\rho }}(\mathbf{\rho },z)=-\frac{2\pi }{s}\frac{\mathbf{k}}{k}\mathbf{p}\cdot \mathbf{k}{e}^{-k|z|}{e}^{i\mathbf{k}\cdot \mathbf{\rho }},$$

$${\mathbf{E}}_z(\mathbf{\rho },z)=-{\mathbf{e}}_z\frac{2\pi i}{s}\mathbf{p}\cdot \mathbf{k}\mathrm{sgn}(z)e^{-k|z|}{e}^{i\mathbf{k}\cdot \mathbf{\rho }}.$$

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

$${\mathbf{E}}_{l,2}=\mathbf{E}+Q_1{\mathbf{p}}_1+Q_0{\mathbf{p}}_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=\sum _{m,n}\frac{1}{2(m^2+n^2+mn+n+1/3{)}^{3/2}{a}^3}\approx \frac{11.5753}{a^3}.$$

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

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

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

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

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

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

$$(1-{\alpha }_1{Q}_0){\mathbf{p}}_1-{\alpha }_1{Q}_1{\mathbf{p}}_2=e_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=-e_a[1+{\alpha }_2(Q_1-Q_0)]{\mathbf{u}}_2+{\alpha }_2\mathbf{E}.$$

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

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

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

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

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

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

$$D=1-({\alpha }_1+{\alpha }_2)Q_0-{\alpha }_1{\alpha }_2(Q_1^2-Q_0^2).$$

$$\mathcal{P}=({\mathbf{p}}_1+{\mathbf{p}}_2)/s,$$

$$\mathbf{w}=\sqrt{\frac{\bar{m}}{s}}({\mathbf{u}}_1-{\mathbf{u}}_2),$$

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

$$a_{21}=\frac{e_{a}}{D\sqrt{\bar{m}s}}\big{[}1+\alpha_{1}(Q_{1}-Q_{0})\big{]}\big{[}1+\alpha_{2}(Q_{1}-Q_{0})\big{]},$$

$$a_{22}=\frac{1}{sD}\big{[}(\alpha_{1}+\alpha_{2})+2\alpha_{1}\alpha_{2}(Q_{1}-Q_{0})\big{]}.$$

$$a_{22}={\chi }_e.$$

$$e_{B}=\frac{e_{a}}{D}\big{[}1+\alpha_{1}(Q_{1}-Q_{0})\big{]}\big{[}1+\alpha_{2}(Q_{1}-Q_{0})\big{]},$$

$$a_{21}=\frac{e_B}{\sqrt{\bar{m}s}}.$$

$$a_{22}={\chi }_e=\frac{1}{s{Q}_1}\left(\frac{e_B}{e_a}-1\right),$$

$$\mathbf{E}_{1}=\frac{1}{D}\big{[}1+\alpha_{2}(Q_{1}-Q_{0})\big{]}\big{[}e_{a}Q_{1}(\mathbf{u}_{1}-\mathbf{u}_{2})+\mathbf{E}\big{]},$$

$$\mathbf{E}_{2}=\frac{1}{D}\big{[}1+\alpha_{1}(Q_{1}-Q_{0})\big{]}\big{[}e_{a}Q_{1}(\mathbf{u}_{1}-\mathbf{u}_{2})+\mathbf{E}\big{]}.$$

$$m_1{\ddot{\mathbf{u}}}_1=-K({\mathbf{u}}_1-{\mathbf{u}}_2)+e_a{\mathbf{E}}_1+e_a{\alpha }_1(Q_1-Q_0){\mathbf{E}}_1,$$

$$m_2{\ddot{\mathbf{u}}}_2=-K({\mathbf{u}}_2-{\mathbf{u}}_1)-e_a{\mathbf{E}}_2-e_a{\alpha }_2(Q_1-Q_0){\mathbf{E}}_2.$$

$$K_e=e_a{e}_B{Q}_1,$$

$$m_1{\ddot{\mathbf{u}}}_1=(K_e-K)({\mathbf{u}}_1-{\mathbf{u}}_2)+e_B\mathbf{E},$$

Full text

I Introduction

Since the discovery of graphene two-dimensional (2D) materials have become a subject of intense research due to their novel mechanical, electronic and optical properties Geim and Grigorieva (2013). Monolayer (ML) hexagonal Boron nitride (hBN) and transition metal dichalcogenides (TMDs) such as MoS2 are two types of prominent 2D binary crystals with a similar honeycomb lattice structure to graphene. While graphene is a purely covalent material, ML hBN is a 2D polar crystal with mixed covalent and ionic bonding. Also distinct from graphene, a semimetal, 2D hBN is an insulator with a large bandgap $\sim$7 eV Rasmussen et al. (2016); Şahin et al. (2009). Further, piezoelectricity occurs in ML hBN owing to inversion symmetry breaking. In recent years ML hBN has attracted intense interest due to its strong piezoelectric response Naumov et al. (2009); Droth et al. (2016); Michel et al. (2017), large in-plane stiffness and strong resistance to stretching Şahin et al. (2009); Song et al. (2010); Andrew et al. (2012); Thomas et al. (2016), and also its capability to combine with graphene due to lattice match as a key material in layered graphene/hBN heterostructure electronic devices Ponomarenko et al. (2011); Haigh et al. (2012).

Optical phonons are key scattering partners of electrons Ridley (2013) and play crucial roles in carrier transport in electronic devices. In a polar crystal long optical lattice vibrations are closely connected to the electric fields Born and Huang (1954). While the macroscopic electric field $\mathbf{E}$ is an average field, with the averaging being made to the total field of all the ions over a lattice cell Born and Huang (1954); Kittel (2004), the local field on any particular ion $\mathbf{E}_{l}$, as was introduced by Lorentz, is the total electric field after deduction of the contribution due to the ion itself. In a bulk polar crystal the Lorentz relation connects the local and macroscopic fields with the macroscopic dielectric polarization $\mathbf{P}$, $\mathbf{E}_{l}=\mathbf{E}+4\pi\mathbf{P}/3$ Born and Huang (1954); Kittel (2004), where the last term is the inner field Born and Huang (1954). The electric field acting on each ion of the crystal polarizes it by inducing an electric dipole moment, thus causing electronic polarization (EP) to the ion, or concisely ionic polarization. A treatment neglecting (including) such EP is termed rigid (polarizable) ion model, i.e., RIM and PIM. Using the Lorentz relation whilst considering the ions to be polarizable with a microscopic model, Huang deduced a pair of equations describing the long optical vibrations of isotropic bulk crystals, $\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}$ describes the optical displacement of the unit cell, $\mathbf{w}=\sqrt{\bar{m}/v_{a}}(\mathbf{u}_{1}-\mathbf{u}_{2})$, $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ being the displacements of the positive and negative ions, $\bar{m}$ being their reduced mass and $v_{a}$ the cell volume. The lattice dielectric function can be deduced directly from Huang’s equations, allowing one to conveniently express the $b$-coefficients in terms of experimentally measurable quantities such as the static and high-frequency dielectric constants $\epsilon_{0}$ $\epsilon_{\infty}$ and infrared dispersion frequency $\omega_{0}$ Born and Huang (1954). Solving Huang’s equations and the equation of electrostatics yields the longitudinal and the transverse optical (LO and TO) modes, with TO phonon frequency $\omega_{t}$ equal to $\omega_{0}$. Further, there is a macroscopic electrostatic field associated with the LO modes making their frequency $\omega_{l}$ higher than $\omega_{t}$, with the frequency ratio given by the Lyddane–Sachs–Teller (LST) relation Lyddane et al. (1941), $\omega_{l}/\omega_{t}=\sqrt{\epsilon_{0}/\epsilon_{\infty}}$, which can also be rederived from Huang’s equations Born and Huang (1954).

In experimental study Rokuta et al. (1997), phonon spectra for ML hBN on Ni and Pt have been measured by electron energy loss spectroscopy. Phonon spectra of ML hBN are usually calculated by diagonalizing the dynamical matrix which is obtained by first-principles calculation, or a tight-binding or an empirical model. Phonon spectra have been studied for ML hBN within the local-density approximation in density functional theory (DFT) Miyamoto et al. (1995); Illera et al. (2017). As a ML of hBN can be rolled up to form a BN nanotube, first-principles Ferrabone et al. (2011) and tight-binding Sánchez-Portal and Hernández (2002) calculations have been performed to make a comparative study of dielectric polarizabilities Ferrabone et al. (2011) as well as phonon spectra Sánchez-Portal and Hernández (2002) of hBN MLs and nanotubes. Phonon modes of 2D hBN calculated from first principles Serrano et al. (2007); Topsakal et al. (2009); Şahin et al. (2009); Michel et al. (2017); Sohier et al. (2017) or an empirical force constant model Michel and Verberck (2009, 2011a, 2011b) have been compared with those of three-dimensional (3D) bulk hBN Serrano et al. (2007); Topsakal et al. (2009); Michel and Verberck (2011a, b) or other 2D honeycomb materials such as TMDs Michel et al. (2017); Sohier et al. (2017) and group III nitrides Şahin et al. (2009). In these studies, there are several common features present in the phonon spectra for ML hBN Sánchez-Portal and Hernández (2002); Serrano et al. (2007); Topsakal et al. (2009); Şahin et al. (2009); Michel and Verberck (2009); Sohier et al. (2017): (i) the LO and TO modes are degenerate at the $\Gamma$ point but split up at a finite wavevector with the LO mode having a higher frequency; (ii) overbending occurs in the LO phonon dispersion curve so the LO modes have maximum frequency not at $\Gamma$, but at an intermediate point away from the Brillouin zone edges; (iii) both TO and ZO (i.e., out-of-plane optical) modes show a nondispersive character at long wavelengths with a nearly constant frequency. While most of these studies are performed by numerical methods, only few studies have used analytical approaches Sánchez-Portal and Hernández (2002); Michel and Verberck (2009); Sohier et al. (2017). The degeneracy has been proved to be due to the macroscopic field’s in-plane component vanishing at zero wavevector Sánchez-Portal and Hernández (2002), and analytical expressions have been obtained later for the long-wavelength dispersion of LO and TO modes Michel and Verberck (2009); Sohier et al. (2017). In Ref.Sohier et al. (2017) an LO phonon dispersion relation is derived using a simple model, in which the relationship between the squared LO and TO phonon frequencies for bulk materials Cochran and Cowley (1962); Giannozzi et al. (1991) is generalized and used for the 2D materials, and the parameters in the dispersion relation are obtained from their first-principles calculation so both ionic EP and LFEs are taken into account. The analytical theory of Michel et al. Michel and Verberck (2009) tackles the dynamical matrix, which is obtained based on a microscopic RIM, to find the vibrational modes; while the eigensolutions of the dynamical matrix are numerically calculated for an arbitrary wavevector, the analytical phonon dispersions are obtained for small wavevectors. In the study Michel and Verberck (2009) ionic polarization has not been accounted for, as the inclusion of EP presents a challenge and makes it more difficult to obtain analytical solutions to the polar optical vibrations. Equally important are the local fields on the ions, which have been found to be very strong in 2D hBN Mikhailov (2013). Further, the local fields and EP are interdependent Born and Huang (1954) and should be included in a self-consistent manner, and therefore the addition of local field effects (LFEs) makes the lattice-dynamical solutions for a general wavevector more complicated. For the long wavelengths, however one expects that the 2D lattice motion can be described on a macroscopic basis, i.e., using macroscopic quantities such as the macroscopic field and dielectric polarization, by lattice equations like Huang’s equations for bulk crystals. Such equations need to be deduced from a microscopic model so as to include the intricate EP and LFEs. So far as we know, there are no equations of motion for the macroscopic description of the 2D lattice vibrations, and ionic EP or LFEs on the lattice dynamical properties of the 2D crystals have not been studied systematically. From the viewpoint of basic research, ML hBN, the structurally simplest 2D polar crystal, provides a model system for analytically studying the 2D optical vibrations and further local field and polarizable ion effects.

In this paper, we study long wavelength optical lattice vibrations and local field and polarizable ion effects for ML hBN with an analytical approach. For this purpose, we deduce lattice equations for optical vibrations using a microscopic model that includes the ionic EP and LFEs. We make the deduction with Huang’s approach, by introducing the macroscopic field into the equation of motion whilst constructing an equation for the macroscopic dielectric polarization by adding up the two contributions due to the lattice displacement and the induced electric polarization. We solve the simultaneous lattice equations and equation of electrostatics, rather than solve the dynamical matrix equations as in previous studies, to obtain explicit expressions for the optical modes, which can describe the key features of the 2D phonon modes. We derive a 2D longitudinal lattice dielectric function $\epsilon(k,\omega)$ which also allows one to rederive the LO phonon dispersion simply from $\epsilon(k,\omega)=0$. We also deduce a LST relation, a 2D counterpart of the LST relation in bulk, for in-plane motion and a frequency–susceptibility relation for out-of-plane motion. Using first-principles calculated quantities we study the lattice dynamical properties for 2D hBN and discuss in great detail the local field and polarizable ion effects.

This paper is organized as follows. In Section II, a deduction of 2D lattice equations for in-plane and out-of-plane motion is made, through a 2D Lorentz relation, from a microscopic dipole lattice model including LFEs and EP. From the lattice equations, the dispersion relations of the optical modes are deduced, followed by a derivation of the phonon group velocity and density of states (DOS). Then the dynamical lattice dielectric susceptibilities are derived yielding the relations relating the static and electronic susceptibilities to the $a_{12}$ and $c_{12}$ coefficients of the lattice equations. Further, a 2D longitudinal lattice dielectric function is deduced after considering the general test charge distributions, and from it a 2D LST relation follows. In Section III, we present results of the in-plane and out-of-plane optical vibrations in ML BN. A comparison of microscopic quantities such as effective charges and spring force constants is made, which are calculated when knowing three quantities including the electronic and static susceptibilities from independent first-principles calculations. Then various lattice-dynamical quantities, such as the Born charge, the phonon dispersion and the static and electronic susceptibilities, are compared with those obtained from a RIM, with or without LFEs taken into account, to study the EP and LFEs on the lattice dynamical properties. Finally, Section V summarizes the main results obtained. In Appendix A, we show that using the 2D Clausius-Mossotti relation obtained, the unit-cell atomic polarizability falls in an interval for in-plane or out-of-plane polarization, which is used to evaluate the LFEs on the phonon dispersion and 2D dielectric susceptibilities. In Appendix B we show that in terms of macroscopic theory the relations $a_{12}=a_{21}$ and $c_{12}=c_{21}$, which connect the coefficients of the lattice equations, follow from the principle of energy conservation. Further we obtain a lattice-vibrational energy density as a function of the optical displacement and electric field for in-plane or out-of-plane vibrations.

II Theory

II.1 Equations of motion and lattice polarization

ML hBN, a 2D binary crystal with point group D3h, is composed of two sublattices of B and N (as shown in Fig. 3), labeled with $\kappa=1,2$, respectively. Let $m_{\kappa}$ and $e_{\kappa}$ be the mass and charge of the type $\kappa$ ions , and let $e_{1}=-e_{2}=e_{a}$, where $e_{a}$ is the static effective charge Karch and Bechstedt (1997) due to electron charge transfer $-e_{a}$ from B to N in 2D h-BN, $e_{a}>0$ Topsakal et al. (2009); Michel et al. (2017). The masses of the B and N atoms are $m_{1}$=10.81 Da and $m_{2}$=14.01 Da. In the dipole lattice model Born and Huang (1954) each ion site of type $\kappa$ is occupied by an electric dipole $\mathbf{p}_{\kappa}$ which arises due partly to the ionic displacement $\mathbf{u}_{\kappa}$ and partly to the induced electric moment $\bm{\mu}_{\kappa}$ on the ion. Associated with the long wavelength optical modes, there is macroscopic dielectric polarization $\mathbf{P}=(\mathbf{p}/s)e^{i\mathbf{k}\cdot\bm{\rho}}\delta(z)$, where $\mathbf{p}$ is the total dipole moment $\mathbf{p}=\mathbf{p}_{1}+\mathbf{p}_{2}$ of a unit cell with area $s=\sqrt{3}a^{2}/2$, $a$ being the lattice constant $a$=2.5 $\AA$ (see Fig. 3). The $\delta$ function describes the dependence of the polarization on $z$ for a ML in the plane $z=0$. $\mathbf{k}$ is the 2D wave vector, and $\bm{\rho}=(x,y)$ is the position vector in the plane parallel to the ML.

We first consider in-plane optical vibrations; that is, the displacements $\mathbf{u}_{\kappa}$ and dipole moments $\mathbf{p}_{\kappa}$ lie in the ML plane. The macroscopic electric field $\mathbf{E}$ due to the charge density $-\nabla\cdot\mathbf{P}$ is given by the equation of electrostatics $\nabla\cdot(\mathbf{E}+4\pi\mathbf{P})=0$, the electric field $\mathbf{E}$ being an irrotational field, $\mathbf{E}=-\nabla\phi$ [$\phi$ is the electrostatic potential, $\phi(\bm{\rho},z)=\varphi(z)e^{i\mathbf{k}\cdot\bm{\rho}}$]. To solve the corresponding Poisson’s equation $\nabla^{2}\phi(\bm{\rho},z)=4\pi i\mathbf{p}\cdot\mathbf{k}e^{i\mathbf{k}\cdot\bm{\rho}}\delta(z)/s$, we expand $\varphi(z)$ and $\delta(z)$,

[TABLE]

[TABLE]

We find $\hat{\varphi}(q)=-2i\mathbf{p}\cdot\mathbf{k}/[s(k^{2}+q^{2})]$ and from Eq. (1) $\varphi(z)=-2\pi i\mathbf{p}\cdot\mathbf{k}e^{-k\lvert z\rvert}/(sk)$, and then obtain for the macroscopic field its components in the directions parallel and perpendicular to the ML:

[TABLE]

Here $\operatorname{sgn}(z)$ is the sign function and in particular the $z$-component of the field is zero in the ML.

Recall that in bulk polar crystals the macroscopic field is strongly nonanalytical close to zero wavevector, i.e., its limiting value depending on the direction along which zero wavevector is approached Born and Huang (1954), and the LO vibrations have a finite field at small wavevectors. Distinct from the field in 3D crystals, the macroscopic field [Eqs. (3a) and (3b)] of 2D hBN vanishes at a very small wavevector for in-plane motion Sánchez-Portal and Hernández (2002), independent of the direction of the wavevector.

In a 2D dipole lattice a dipole $\mathbf{p}$ on the ion at a lattice point $\bm{\rho}_{i}=(x_{i},y_{i})$ gives rise to a field at the origin equal to $-\mathbf{p}/\rho_{i}^{3}+3\mathbf{p}\cdot\bm{\rho}_{i}\bm{\rho}_{i}/\rho_{i}^{5}$. The Lorentz local field $\mathbf{E}_{l}$, also called the exciting field Born and Huang (1954), is the electric field acting on an ion due to all the other dipoles oscillating in the lattice, and can be written as the macroscopic field $\mathbf{E}$ plus another field $\mathbf{E}_{in}$ (namely, the inner field) Born and Huang (1954), $\mathbf{E}_{l}=\mathbf{E}+\mathbf{E}_{in}$. In a long wavelength lattice wave, the local fields at the B and N lattice points are given by the following expressions Mikhailov (2013), respectively,

[TABLE]

with the coefficients $Q_{0}$ and $Q_{1}$ Mikhailov (2013); Dell’Anna and Merano (2016)

[TABLE]

Let dipole $\mathbf{p}_{1}$ or equally $\mathbf{p}_{2}$ have 0.01 $e\AA$, for instance, corresponding to a displacement of 0.01 $\AA$ of the ions with charge $e$, then from these expressions one finds very strong local fields with a magnitude 1000 kV/cm Mikhailov (2013).

The LFEs are included through coefficients $Q_{0}$ and $Q_{1}$. Note that $\mathbf{E}$ is the macroscopic field [Eqs. (3a) and (3b)] rather than an external field in a simple sense Mikhailov (2013). Such a simple relation between the local field and the macroscopic field is valid only for the long lattice waves, while the general expression for an arbitrary wavelength is quite complicated with both coefficients $Q_{0}$ and $Q_{1}$ dependent on the wavevector. In 3D iostropic polar crystals Born and Huang (1954) the difference between the local and macroscopic fields is proportional to the macroscopic polarization directly (the Lorentz relation), $\mathbf{E}_{l}=\mathbf{E}+4\pi\mathbf{P}/3$, as the $Q$ coefficients are equal to $4\pi/(3v_{a})$ at $\mathbf{k}=0$ Born and Huang (1954). When approximating $Q_{0}=Q_{1}$, Eqs. (4a) and (4b) can be transformed into a simple Lorentz relation involving macroscopic areal polarization $(\mathbf{p}_{1}+\mathbf{p}_{2})/s$, and also a Clausius-Mossotti relation can be deduced for 2D BN (Appendix A). Expressions (4a) and (4b), the Lorentz relations for 2D BN, show that the finite local fields occur in ML hBN at $\mathbf{k}=0$, different from those in the 3D polar crystals where local fields vanish in the long wavelength limit Born and Huang (1954).

Apart from the macroscopic field and local fields, we also need to find the field change at the center of an ion of type $\kappa$ owing to its own displacement $\mathbf{u}_{\kappa}$ Born and Huang (1954). The field is evidently equal to the field created at the ion $\kappa$ site by displacing all other ions by $-\mathbf{u}_{\kappa}$. Hence it is equal to the local field at that ion site in a dipole lattice with displacement dipoles $\mathbf{p}_{\kappa^{\prime}}=-e_{\kappa^{\prime}}\mathbf{u}_{\kappa}$, where type $\kappa^{\prime}=1,2$. Substituting this dipole expression into Eqs. (4a) and (4b) and putting $\mathbf{E}=0$ as wavevector $\mathbf{k}=0$, we find the field changes at the centers of the B and N ions due to their own displacements respectively,

[TABLE]

The total Coulomb fields $\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ at the centers of the B and N ions are simply the sums of $\mathbf{E}_{l,1}$ and $\mathbf{E}_{u,1}$ [Eqs. (4a) and (6a)], and $\mathbf{E}_{l,2}$ and $\mathbf{E}_{u,2}$ [Eqs. (4b) and (6b)], respectively,

[TABLE]

where all the vectors are in the layer plane.

Assuming that the electronic polarization of an ion is equivalent to a point-dipole Born and Huang (1954), the electronic (i.e., induced) dipole moment of the ion $\kappa$ is then given by $\bm{\mu}_{\kappa}=\alpha_{\kappa}\mathbf{E}_{\kappa}$, where $\alpha_{\kappa}$ is the in-plane electronic polarizability of the ion $\kappa$. Then the total dipole moments on the B and N ions are

[TABLE]

Inserting the expressions for the total fields $\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ into Eqs. (8a) and (8b) and then rearranging the terms, we find

[TABLE]

Solving Eqs. (9a) and (9b) then we can express $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ in terms of $\mathbf{u}_{1}$, $\mathbf{u}_{2}$ and $\mathbf{E}$ as follows:

[TABLE]

where

[TABLE]

Define the areal polarization $\bm{\mathcal{P}}$

[TABLE]

and introduce the optical displacement $\mathbf{w}$

[TABLE]

where $\bar{m}$ is the reduced mass, $\bar{m}=m_{1}m_{2}/(m_{1}+m_{2})$.

When expressions (10a) and (10b) are substituted for $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$, we obtain

[TABLE]

where

[TABLE]

Eq. (14) shows that the macroscopic quantity $\bm{\mathcal{P}}$ of the 2D crystal, which has a clear physical meaning as given by expressions (12), is simplified 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. (14) that $a_{22}$ is the in-plane component $\chi_{e}$ of the electronic susceptibility of the 2D crystal,

[TABLE]

When the Born charge Gonze and Lee (1997)

[TABLE]

is introduced, then the coefficient $a_{21}$ relates simply to the Born charge $e_{B}$ by

[TABLE]

From Eqs. (15b), (16) and (17) we find

[TABLE]

showing that apart from $a_{21}$, the coefficient $a_{22}$ is also related to the Born charge. $e_{B}\neq e_{a}$ owing to the electronic polarization of the ions, and further considering $Q_{1}>0$, $e_{B}$ is greater than $e_{a}$.

Inserting expressions (10a) and (10b) for $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ into Eqs. (7a) and (7b), we express the total fields $\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ in terms of the ionic displacements and macroscopic field. After collecting like terms we can simplify the total field expressions to

[TABLE]

The Coulomb force acting on the ion $\kappa$ consists of two parts Born and Huang (1954), (i) the force exerted on the ionic charge $e_{\kappa}$ by the total field $\mathbf{E}_{\kappa}$, and (ii) the force exerted on the induced dipole $\bm{\mu}_{\kappa}$ by the field of all other ions. The latter electric force can be sought by using again the Lorentz relations (4a) and (4b) as follows. Imagine that we subject the dipole $\bm{\mu}_{\kappa}$ at a particular site of ion $\kappa$ a virtual displacement $\mathbf{u}$ while keeping all other ions in their undisplaced positions. The virtual energy is then simply the interaction energy between the dipole and the field at the ion site which is created equivalently by displacing all other ions by $-\mathbf{u}$ Born and Huang (1954), and which is thus the local field in a dipole lattice with dipole moments $\mathbf{p}_{\kappa^{\prime}}=-e_{\kappa^{\prime}}\mathbf{u}$, (type $\kappa^{\prime}=1,2$). Therefore the virtual energy is $-\bm{\mu}_{\kappa}\cdot\mathbf{E}_{l,\kappa}$, and the force on the dipole is the negative gradient of this virtual energy with respect to virtual displacement $\mathbf{u}$, given by $\nabla_{\mathbf{u}}(\bm{\mu}_{\kappa}\cdot\mathbf{E}_{l,\kappa})$. Inserting the above dipole moment expression into Eqs. (4a) and (4b) and recalling the macroscopic field $\mathbf{E}=0$, we find the local fields and then obtain for the forces on the dipoles $\bm{\mu}_{1}$ and $\bm{\mu}_{2}$ the expressions $e_{a}(Q_{1}-Q_{0})\bm{\mu}_{1}$ and $-e_{a}(Q_{1}-Q_{0})\bm{\mu}_{2}$, respectively, where $\bm{\mu}_{\kappa}=\alpha_{\kappa}\mathbf{E}_{\kappa}$ ($\kappa=1,2$).

Apart from the electric forces, there are also restoring forces caused by the overlap potential between the positive and negative ions in the 2D polar crystal. The restoring forces on the B and N ions are $-K(\mathbf{u}_{1}-\mathbf{u}_{2})$, and $-K(\mathbf{u}_{2}-\mathbf{u}_{1})$, respectively, where the spring force constant $K$ for in-plane motion is a simple scalar rather than a tensor due to the hexagonal symmetry.

Therefore the equations of motion for the B and N ions are given by

[TABLE]

Substitute expressions (20a) and (20b) for $\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ respectively in the two equations above. Using the Born charge $e_{B}$ [expression (17)] and introducing a force constant due to LFEs $K_{e}$,

[TABLE]

the equations of motion reduce to

[TABLE]

Multiplying Eqs. (23a) and (23b) by $m_{2}$ and $m_{1}$ respectively, subtracting and then dividing by $(m_{1}+m_{2})$, we find

[TABLE]

On expressing $(\mathbf{u}_{1}-\mathbf{u}_{2})$ in terms of $\mathbf{w}$ [Eq. (13)], we obtain

[TABLE]

where

[TABLE]

$\omega_{0}$ being the intrinsic oscillator frequency, i.e., in the absence of macroscopic field $\mathbf{E}$.

Comparing Eqs. (18) and (26b) we find the relation

[TABLE]

The equation of motion (25) and the polarization equation (14) describe the in-plane polar optical vibrations of 2D BN, the frequencies of which will be derived in Sec. III below. The lattice vibrations we are dealing with are of long wavelengths, and the pair of equations (25) and (14) constitute a macroscopic description of the lattice motion. It is shown in Appendix B that, from the viewpoint of the macroscopic theory, the relation (27) is due to the principle of energy conservation, and further this relation makes it possible to define an areal energy density [expression (121)] from which the lattice equations (25) and (14) can be rederived. It is evident from Eqs. (23a) and (23b) that the center of mass of the two atoms remains stationary (frequency $\omega=0$), yielding trivial nondynamical solutions. We note that for clearness $\mathbf{E}$ appearing in all the equations above is used to represent the in-plane component of the macroscopic field in the ML, i.e., $\mathbf{E}=\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$.

Now we consider ionic motion perpendicular to the layer plane, in which case the displacements $\mathbf{u}_{\kappa}$ and dipole moments $\mathbf{p}_{\kappa}$ are parallel to $\mathbf{e}_{z}$. Let $\alpha_{\kappa}^{\prime}$ denote the electronic polarizability of the type $\kappa$ ions; note that in general $\alpha_{\kappa}^{\prime}\neq\alpha_{\kappa}$ (the latter is the polarizability for in-plane motion), as they are simply components of the polarizability tensor Born and Huang (1954). Let $K^{\prime}$ be the force constant associated with the perpendicular motion. Considering the anisotropic 3D charge density distribution, another static effective charge $e_{a}^{\prime}$ exists likewise, which may differ from $e_{a}$ for in-plane motion (anisotropy), and is needed in the point-ion model for the perpendicular motion; thus the charges on the ions are $e_{1}^{\prime}=-e_{2}^{\prime}=e_{a}^{\prime}$.

Solving the Poisson equation yields the electrostatic potential

[TABLE]

and then the $z$- and in-plane components of the macroscopic field follow,

[TABLE]

We note that (i) there is a $\delta(z)$ term in the $z$-component of the macroscopic field, which reflects the microscopic character of the polarization in terms of its atomic scale when the dipoles point to the $z$ direction. This term follows the dielectric polarization $\mathbf{P}$ according to $-4\pi\mathbf{P}$, which is also true when the $\delta(z)$ function of $\mathbf{P}$ is generalized to an arbitrary function $f(z)$; (ii) the in-plane component of the field is zero in the ML.

The local fields at the B and N sites are given by

[TABLE]

respectively, where the coefficients $Q_{0}^{\prime}$ and $Q_{1}^{\prime}$ are negative, $Q_{0}^{\prime}=-2Q_{0}$, and $Q_{1}^{\prime}=-2Q_{1}$ Mikhailov (2013). The long wavelength macroscopic field has a $\delta(z)$ form [Eq. (29a)] but in calculating the field change experienced by a type $\kappa$ ion owing to its own displacement $\mathbf{u}_{\kappa}$, the contribution of this macroscopic field is zero, i.e., $4\pi\mathbf{u}_{\kappa}\sum_{\kappa^{\prime}}e_{\kappa^{\prime}}^{\prime}\delta(z)/s=4\pi\mathbf{u}_{\kappa}(e_{1}^{\prime}+e_{2}^{\prime})\delta(z)/s=0$. Thus the field changes at the centers of B and N owing to their own displacements are

[TABLE]

respectively. Similarly, in calculating the force exerted on the induced dipole $\bm{\mu}_{\kappa}$ at a type $\kappa$ ion by the field of all other ions, the macroscopic field makes no contribution to this force again because the net charge per cell vanishes. Then repeating the process as before, we find that with replacements $Q_{0}\rightarrow Q_{0}^{\prime}$, $Q_{1}\rightarrow Q_{1}^{\prime}$, $e_{a}\rightarrow e_{a}^{\prime}$, $\alpha_{\kappa}\rightarrow\alpha_{\kappa}^{\prime}$ and $K\rightarrow K^{\prime}$, the equations and expressions above for the 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 as given by Eq. (13) and $c_{22}$ is the out-of-plane component $\chi_{e}^{\prime}$ of the electronic susceptibility of the 2D crystal,

[TABLE]

and $c_{21}$ and $c_{22}$ relate to the Born charge via

[TABLE]

The equation of motion is given by

[TABLE]

where

[TABLE]

$\omega_{0}^{\prime}$ being the intrinsic oscillator frequency.

In the lattice equations (37) and (34) $\mathbf{E}$ is the field in the ML, $\mathbf{E}(\bm{\rho},0)$, and evidently $\mathbf{E}(\bm{\rho},0)=\mathbf{E}_{z}(\bm{\rho},0)$ for the out-of-plane vibrations. The areal energy density associated with the out-of-plane optical vibrations is given by expression (136).

When the equations for in-plane motion [Eqs. (25) and (14)] and out-of-plane motion [Eqs. (37) and (34)] are considered simultaneously, they can be rewritten for clarity as

[TABLE]

and

[TABLE]

respectively, where $a_{12}=a_{21}$ and $c_{12}=c_{21}$. These equations have similar forms to Huang’s equations for bulk crystals Huang (1951a); Born and Huang (1954).

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

The in-plane optical vibration modes can be obtained from Eqs. (25) and (14) in conjunction with the equation of electrostatics $\nabla\cdot(\mathbf{E}+4\pi\mathbf{P})=0$, where $\mathbf{P}$ is the dielectric polarization (namely, a dipole moment per unit volume), $\mathbf{P}=\bm{\mathcal{P}}\delta(z)$, and $\mathbf{E}$ is an irrotational field, $\mathbf{E}=-\nabla\phi$. Let $\mathbf{w}(\bm{\rho})=\mathbf{w}_{0}e^{i\mathbf{k}\cdot\bm{\rho}}$ and electrostatic potential $\phi(\bm{\rho},z)=\varphi(z)e^{i\mathbf{k}\cdot\bm{\rho}}$ (time dependence $e^{-i\omega t}$ is omitted for clearness). Expressing the field $\mathbf{E}=-\nabla\phi$ in terms of $\varphi$ gives the in-plane component $\mathbf{E}_{\bm{\rho}}(\bm{\rho},z)=-i\mathbf{k}\varphi(z)e^{i\mathbf{k}\cdot\bm{\rho}}$ and further this field in the ML $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)=-i\mathbf{k}\varphi(0)e^{i\mathbf{k}\cdot\bm{\rho}}$. Then apply divergence $\nabla_{\bm{\rho}}$ to Eq. (14), and we have the polarization charge density $-\nabla\cdot\mathbf{P}=-\delta(z)\nabla_{\bm{\rho}}\cdot\bm{\mathcal{P}}=-\delta(z)[a_{21}i\mathbf{k}\cdot\mathbf{w}+a_{22}k^{2}\varphi(0)e^{i\mathbf{k}\cdot\bm{\rho}}]$ and obtain Poisson’s equation,

[TABLE]

To solve Eq. (41), we insert the expansions of $\varphi(z)$ [Eq. (1)] and $\delta(z)$ [Eq. (2)], yielding

[TABLE]

When this expression is substituted for $\hat{\varphi}(q)$ of the following equation,

[TABLE]

we find

[TABLE]

and then obtain the electric field in the ML,

[TABLE]

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

[TABLE]

The expression (45) admits two possibilities for $\mathbf{w}\cdot\mathbf{k}$, namely, case (i) $\mathbf{w}\cdot\mathbf{k}=0$, or case (ii) $\mathbf{w}\cdot\mathbf{k}\neq 0$. In case (i), as $\mathbf{w}\cdot\mathbf{k}=0$ the normal modes are transverse waves. Eq. (44) gives $\varphi(0)=0$, and on account of Eq. (42) we have $\hat{\varphi}(q)=0$. Then it follows from Eq. (1) that $\varphi(z)=0$, and therefore in the electrostatic approximation the macroscopic field vanishes identically, $\mathbf{E}(\mathbf{r})=0$. The frequency of the TO mode of wavevector $\mathbf{k}$ is given by the solution of Eq. (46) with the field equal to zero,

[TABLE]

which is independent of wavevector; that is, the long-wavelength TO modes are dispersionless, consistent with previous tight-binding Sánchez-Portal and Hernández (2002) and first-principles Wirtz et al. (2003); Topsakal et al. (2009); Sohier et al. (2017) calculations.

In case (ii) the electric field $\mathbf{E}_{\bm{\rho}}(\mathbf{r})$ is nonzero, and from the equations 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., parallel to wave vector $\mathbf{k}$, $\mathbf{w}(\bm{\rho})\parallel\mathbf{E}_{\bm{\rho}}(\mathbf{r})\parallel\mathbf{P}(\mathbf{r})\parallel\mathbf{k}$. Inserting expression (45) for $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ into Eq. (46) we find the frequency of the longitudinal optical (LO) mode,

[TABLE]

Expressing the $a$-coefficients in terms of $e_{B}$, $\chi_{e}$, $\omega_{0}$ by Eqs. (16), (18) and (47), then $\omega_{l}(k)$ can be rewritten as

[TABLE]

Physically, Eq. (49) is valid for small wavevectors; mathematically $\omega_{l}$ is a monotonically increasing function of $k$ and has an upper limit $\omega_{M}$ for very large $k$,

[TABLE]

In the theoretical study by Sohier et al. Sohier et al. (2017) the LO phonon dispersion is expressed as $\omega_{l}^{2}=\omega_{0}^{2}+\mathcal{S}k/(1+r_{eff}k)$. From their Eqs. (2) and (3)111There is a misprint in Eq. (2) of Ref.Sohier et al. (2017): the square should be for each summand rather than for the total sum. one finds that the parameter $\mathcal{S}$ relates to the Born charge via $\mathcal{S}=2\pi e_{B}^{2}/(\bar{m}s)$. $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), where $\epsilon_{p}$ and $t$ are 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). According to the defintion in Refs.Cudazzo et al. (2011); Berkelbach et al. (2013), the effective screening length is given by $2\pi\chi_{e}$ or $2\pi\chi_{0}$ when including the vibrational contribution [$\chi_{0}$ is the static susceptibility Eq. (74)]. As it is the high-frequency dielectric constant that determines the difference between the squared LO and TO phonon frequencies Cochran and Cowley (1962); Giannozzi et al. (1991); Gonze and Lee (1997) so $r_{eff}$ takes the former, i.e., $r_{eff}=2\pi\chi_{e}$. Therefore the LO phonon dispersion Eq. (49) is identical to the analytical expression of Sohier et al..

Having $\varphi(0)$ [Eq. (44)] for the LO mode $\omega_{l}(k)$, we substitute expression (42) back into Eq. (1) to give the electric potential and then obtain the macroscopic field associated with the LO mode, expressed in terms of mode displacement $\mathbf{w}$ [in a form consistent with Eqs. (3a) and (3b)],

[TABLE]

The long-range electrostatic interactions cause a higher LO frequency than the TO frequency for a finite $k$ Sánchez-Portal and Hernández (2002) with the splitting determined by the latter term in the square brackets of Eq. (49). This term also shows that there is no splitting in the limit $k\rightarrow 0$ as the macroscopic field vanishes [Eq. (51a)], which is different from the situation in bulk BN where LO-TO splitting occurs also at the $\Gamma$ point Ohba et al. (2001); Topsakal et al. (2009). Therefore, the transparent expressions (47) and (49) show the degeneracy of the LO and TO modes at $\Gamma$ and their splitting at a finite wavevector, 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); Sohier et al. (2017).

The out-of-plane optical vibrations, namely, the ZO modes, can be obtained from Eqs. (37) and (34) with the electrostatic approach as follows. Recall that the $z$-component of the field [Eq. (29a)] has a $\delta(z)$ term, divergent at $z=0$, because the ML is treated as a geometric plane where the ionic charge distribution and polarization density $\mathbf{P}$ have a $\delta(z)$ form. We approximate $\delta(z)$ by a Gaussian distribution with a small thickness $\varepsilon$ ($\varepsilon\rightarrow 0$ ), $\delta_{\varepsilon}(z)=\frac{1}{\sqrt{\pi}\varepsilon}e^{-z^{2}/\varepsilon^{2}}$, as in the theoretical study Michel and Verberck (2009). A similar treatment is used also in first-principles calculations Sohier et al. (2017). The $z$-component of the field in the ML is $\mathbf{E}_{z}(\bm{\rho},0)=-\mathbf{e}_{z}\varphi^{\prime}(0)e^{i\mathbf{k}\cdot\bm{\rho}}$, and after inserting this field into Eq. (34) the polarization charge density can be expressed as $-\nabla\cdot\mathbf{P}=-\delta_{\varepsilon}^{\prime}(z)\mathbf{e}_{z}\cdot\bm{\mathcal{P}}=-[c_{21}\mathbf{w}\cdot\mathbf{e}_{z}-c_{22}\varphi^{\prime}(0)e^{i\mathbf{k}\cdot\bm{\rho}}]\delta_{\varepsilon}^{\prime}(z)$. Thus Poisson’s equation is given by

[TABLE]

Expanding $\varphi(z)$ [Eq. (1)] and $\delta_{\varepsilon}(z)$

[TABLE]

then we have

[TABLE]

[TABLE]

Inserting the expansions (1) and (55) into Poisson’s equation (52) we find $\hat{\varphi}(q)$,

[TABLE]

and then carrying out the integration in Eq. (54) we obtain

[TABLE]

where $\operatorname{erf}(x)$ is the error function, $\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt$.

Taking $z=0$ in Eq. (57) we find $\varphi^{\prime}(0)$,

[TABLE]

where

[TABLE]

$\operatorname{erfc}(x)$ being the complementary error function, $\operatorname{erfc}(x)=1-\operatorname{erf}(x)$.

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

[TABLE]

When this field is substituted into the equation of motion (37), we obtain the frequency of the ZO mode,

[TABLE]

Considering that both $\varepsilon$ and $k$ are small quantities, expression (59) reduces to a constant $\zeta=4\sqrt{\pi}$, and the ZO phonon frequency becomes independent of wavevector,

[TABLE]

for a planar material when taking $\varepsilon\ll 4\sqrt{\pi}\chi_{e}^{\prime}$ the frequency becomes

[TABLE]

The ionic polarization (i.e., $\chi_{e}^{\prime}\neq 0$) ensures a finite frequency; otherwise $\omega_{z}$ [Eq. (62)] becomes extremely large and even $\omega_{z}\rightarrow\infty$ [Eq. (63)] when $\chi_{e}^{\prime}$ is neglected in the RIM (also see Table 4). The result here that the long wavelength ZO modes are dispersionless agrees with previous calculations which show a very flat dispersion curve at small wavevectors Sánchez-Portal and Hernández (2002); Wirtz et al. (2003); Michel and Verberck (2009); Topsakal et al. (2009).

Having $\varphi^{\prime}(0)$ [Eq. (58)] for a ZO mode we then get $\varphi^{\prime}(z)$ from Eq. (57) and also obtain $\varphi(z)$ after substituting the $\hat{\varphi}(q)$ expression (56) into Eq. (1). It follows that the macroscopic field associated with the ZO mode is given by

[TABLE]

which for small $k$ and $\varepsilon$ reduce to the following expressions [compare to Eqs. (29a) and (29b)],

[TABLE]

In particular we note that the in-plane component is an odd function of $z$, and negligible at very small $k$ and $\varepsilon$.

II.3 Phonon group velocity and density of states

Expanding $\omega_{l}(k)$ [Eq. (49)] to second order in wavevector $k$, the LO phonon frequency near the $\Gamma$ point is given by

[TABLE]

where

[TABLE]

which is identical to that given in Ref.Michel and Verberck (2009).

From the dispersion Eq. (49) we find straightforward the DOS of the LO modes,

[TABLE]

where $\omega_{M}$ is the upper bound of LO phonon frequency [expression (50)].

Similarly, expanding $\omega_{z}(k)$ [Eq. (61)] we find the ZO phonon frequency near $\Gamma$,

[TABLE]

where $\omega_{z}$ is the ZO mode frequency given by Eq. (62), and

[TABLE]

Dropping the terms due to the EP of ions, expansions (66) and (69) then reduce to a form as given by the RIM of Ref.Michel and Verberck (2009) [Eqs. (57) and (59) therein]. The nearly flat dispersion of the ZO modes allows their DOS to be approximated by Michel and Verberck (2009)

[TABLE]

Evidently $c_{z}$ is close to zero as $\varepsilon\ll 4\sqrt{\pi}\chi_{e}^{\prime}$ for the 2D crystal, resulting in the ZO modes being nondispersive.

The group velocity is $\nabla_{\mathbf{k}}\omega_{i}(k)=\omega_{i}^{\prime}(k)\mathbf{k}/k$, where $i$ indexes a phonon branch of LO, TO or ZO modes so the $c_{l}$ and $c_{z}$ above are simply the norms of the group velocities of the LO and ZO modes at $\Gamma$, respectively, corresponding to the slopes Sohier et al. (2017) of the phonon dispersion curves. For the long lattice waves, as the wavevector $k$ increases, clearly $\omega_{l}$ increases and $\omega_{z}$ decreases, whereas $\omega_{t}$ stays flat near $\Gamma$ [Eq. (47)].

II.4 2D lattice dielectric susceptibility and dielectric function

The lattice susceptibility of ML hBN for in-plane polarization can be deduced from Eqs. (25) and (14) by considering periodic solutions, $\mathbf{W}$, $\mathbf{E}$, $\bm{\mathcal{P}}$ $\propto e^{-i\omega t}$, due to an external charge or electromagnetic field with oscillation frequency $\omega$. Then Eq. (25) reduces simply to Eq. (46). Inserting Eq. (46) into Eq. (14) to eliminate $\mathbf{W}$, one finds

[TABLE]

and subsequently the 2D dielectric susceptibility, defined by $\chi=\bm{\mathcal{P}}(\bm{\rho})/\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$,

[TABLE]

where $\chi_{e}$ is the high-frequency (i.e., clamped-ion, or electronic) susceptibility.

The static dielectric susceptibility $\chi_{0}$ follows from expression (73),

[TABLE]

where the second term on the right-hand side is the vibrational (also called ionic) contribution. With Eq. (74), $a_{12}$ [Eq. (26b)] can be expressed in terms of the 2D susceptibilities,

[TABLE]

The 2D dynamical susceptibility $\chi(\omega)$ can be expressed in terms of three quantities that are usually used, namely, the intrinsic oscillator frequency $\omega_{0}$ and the high-frequency and static susceptibilities $\chi_{e}$, $\chi_{0}$,

[TABLE]

which has a similar form to its counterpart of bulk polar crystals Born and Huang (1954).

Similarly, the 2D dielectric susceptibility for the vertical $z$-polarization $\chi^{\prime}$ (the prime does not indicate a derivative), defined by $\chi^{\prime}=\bm{\mathcal{P}}(\bm{\rho})/\mathbf{E}_{z}(\bm{\rho},0)$ for any particular $\omega$, is deduced from Eqs. (37) and (34),

[TABLE]

with the static dielectric susceptibility $\chi_{0}^{\prime}$,

[TABLE]

Combining Eqs. (78) and (38b), one can express $c_{12}$ in terms of the 2D susceptibilities,

[TABLE]

The 2D dynamical susceptibility $\chi^{\prime}(\omega)$ can also be expressed in terms of three quantities $\chi_{e}^{\prime}$, $\chi_{0}^{\prime}$, $\omega_{0}^{\prime}$,

[TABLE]

On eliminating $e_{B}^{\prime 2}/(\bar{m}s)$ in Eq. (63) with the help of Eq. (78) we find a simple relation

[TABLE]

The 2D high-frequency and static dielectric susceptibilities can be calculated from first principles, and knowing their values the susceptibility expressions (74) and (78) obtained here will be used below in Sec. III to determine the model parameters such as the effective charges and spring force constants, and the $a$- or $c$-coefficients of the lattice equations [Eqs. (39a), (39b), (40a) and (40b)].

The lattice dielectric function (DF) can be obtained by considering the response of the 2D lattice to a test charge in the electrostatic approximation. Now the equation of electrostatics is given by $\nabla\cdot(\mathbf{E}+4\pi\mathbf{P})=4\pi\sigma$, where $\sigma$ is the test charge density function dependent on $\mathbf{r}$, $t$, and $\mathbf{E}$ is the total field of the test charge plus the polarization charge. To find the dielectric response let $\sigma$ have the form $\sigma=\sigma_{0}e^{i\mathbf{k}\cdot\bm{\rho}}f(z)$, where $f(z)$ is an arbitrary function to describe the charge distribution in the $z$ direction, for the general case that the external charge is not necessarily confined in the ML, and time dependence $e^{-i\omega t}$ has been absorbed into $\sigma_{0}$ for clearness.

First, let us consider a charge distribution that is asymmetric with respect to the ML, i.e., $f(-z)\neq f(z)$, for instance, when the test charge is simply put above or below the ML. Then the field due to the charge is nonzero in the ML, $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)\neq 0$, $\mathbf{E}_{z}(\bm{\rho},0)\neq 0$, and therefore the lattice responds creating both in-plane [Eqs. (39a) and (39b)] and out-of-plane [Eqs. (40a) and (40b)] vibrations, with all quantities such as $\mathbf{w}_{\bm{\rho}}$, $\mathbf{w}_{z}$, $\bm{\mathcal{P}}_{\bm{\rho}}$, $\bm{\mathcal{P}}_{z}$, $\mathbf{E}$ varying according to $e^{i(\mathbf{k}\cdot\bm{\rho}-\omega t)}$. Now the dielectric polarization is given by the sum of the in-plane and out-of-plane contributions, $\mathbf{P}=(\bm{\mathcal{P}}_{\bm{\rho}}+\bm{\mathcal{P}}_{z})\delta(z)$, where $\bm{\mathcal{P}}_{\bm{\rho}}=\chi(\omega)\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ and $\bm{\mathcal{P}}_{z}=\chi^{\prime}(\omega)\mathbf{E}_{z}(\bm{\rho},0)$, with $\chi$ and $\chi^{\prime}$ being the 2D susceptibilities [Eqs. (73) and (77)]. Writing $\mathbf{E}=-\nabla\phi$ with the total electrostatic potential $\phi(\bm{\rho},z)=\varphi(z)e^{i\mathbf{k}\cdot\bm{\rho}}$, one expresses the polarization charge density $-\nabla\cdot\mathbf{P}$ in terms of the derivatives of potential $\phi$, and then Poisson’s equation follows from the equation of electrostatics,

[TABLE]

where $\phi(\bm{\rho},0)$ is the total potential in the ML, $\phi(\bm{\rho},0)=\varphi(0)e^{i\mathbf{k}\cdot\bm{\rho}}$. Inserting the expansion of $f(z)$,

[TABLE]

and expansions (1) and (2) into Poisson’s equation (82) one finds $\hat{\varphi}(q)$,

[TABLE]

where the three terms on the right-hand side represent the respective contributions due to the external charge, in-plane and out-of-plane lattice polarization. Here $\hat{\varphi}(q)$ is expressed in terms of $\varphi(0)$ and $\varphi^{\prime}(0)$ as they are proportional to the field components in the ML, i.e., $\varphi(0)\propto\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ and $\varphi^{\prime}(0)\propto\mathbf{E}_{z}(\bm{\rho},0)$. Integrating $\hat{\varphi}(q)$ over $q$ according to Eq. (43), one then obtains $\varphi(0)$,

[TABLE]

where

[TABLE]

The out-of-plane (ZO) motion makes no contribution to $\varphi(0)$ as the integration value of the third term of $\hat{\varphi}(q)$ is zero. Clearly $\varphi_{\sigma}(0)e^{i\mathbf{k}\cdot\bm{\rho}}$ is the electrostatic potential in the ML due to the test charge alone, and therefore the dielectric function of the 2D lattice, which is defined as the ratio of this test charge potential to the total potential in the ML $\varphi(0)e^{i\mathbf{k}\cdot\bm{\rho}}$ (extension of the 2D wavevector-dependent dielectric function $\epsilon(k)$ of Ref.Cudazzo et al. (2011)), is given by

[TABLE]

showing a linear dependence on wavevector as the $\epsilon(k)$ of Ref.Cudazzo et al. (2011)).

In expression (84) for $\hat{\varphi}(q)$, $\varphi(0)$ is known [Eq. (85)], and $\varphi^{\prime}(0)$ can be determined through expression (54) for $\varphi^{\prime}(z)$ as follows. Substituting expression (84) for $\hat{\varphi}(q)$ in Eq. (54) one finds

[TABLE]

The square bracketed part is a familiar form that has appeared in Eq. (29a). To find $\varphi^{\prime}(0)$ we introduce an effective thickness $\varepsilon$ to $\delta(z)$ to approach it with $\delta_{\varepsilon}(z)$ as was done above in Sec. III. Let $\Lambda=\delta_{\varepsilon}(0)$; $\Lambda=1/(\sqrt{\pi}\varepsilon)$, for instance, when $\delta(z)$ is approximated by a Gaussian Michel and Verberck (2009), $\delta_{\varepsilon}(z)=e^{-z^{2}/\varepsilon^{2}}/(\sqrt{\pi}\varepsilon)$. Now taking $z=0$ in Eq. (88), the term containing $\chi(\omega)\varphi(0)$ vanishes, i.e., the in-plane motion makes no contribution to $\varphi^{\prime}(0)$, which is given by

[TABLE]

Having $\hat{\varphi}(q)$ [expression (84)] then the total potential $\phi$ is known, and the total field $\mathbf{E}$ can be obtained straightforward,

[TABLE]

where $\varphi(0)$ and $\varphi^{\prime}(0)$ are given by expressions (85) and (89), respectively.

Recalling $\mathbf{E}_{\bm{\rho}}(\bm{\rho},z)=-i\mathbf{k}\phi(\mathbf{r})$ with Eq. (90a) for $\mathbf{E}_{\bm{\rho}}(\bm{\rho},z)$, evidently the induced potential associated with the ZO polarization is proportional to $\chi^{\prime}(\omega)\varphi^{\prime}(0)\operatorname{sgn}(z)e^{-k\lvert z\rvert}$, which is zero in the ML plane and antisymmetric with respect to it, consistent with Eq. (28) above. The in-plane component of the field associated with the ZO motion [the term containing $\chi^{\prime}(\omega)\varphi^{\prime}(0)$ in Eq. (90a)] is zero in the ML and thus the ZO motion does not influence the LO motion; meanwhile the LO motion has zero $z$-component of its field [the $\chi(\omega)\varphi(0)$ term in Eq. (90b)] in the ML and thus no influence on the ZO motion, making the LO and ZO vibrations driven by the exernal charge essentially independent of each other. For small thickness $\varepsilon$, $\varphi^{\prime}(0)$ is small from Eq. (89) and there is only a small ZO component, which becomes negligible in the limit $\varepsilon\rightarrow 0$, in the driven lattice motion [refer to Eq. (82)]. For a symmetric test charge distribution $f(-z)=f(z)$, $\mathbf{E}_{z}(\bm{\rho},0)=0$, and this causes only in-plane vibration, with displacement $\mathbf{w}$ and polarization $\mathbf{P}$ both parallel to the ML plane, and evidently yields the same DF $\epsilon(k,\omega)$ as given by expression (87) above.

As the polarization charge density associated with the in-plane motion is $-\delta(z)\nabla_{\bm{\rho}}\cdot\bm{\mathcal{P}}_{\bm{\rho}}(\bm{\rho})=-\delta(z)[a_{21}-a_{22}(\omega^{2}+a_{11})/a_{12}]i\mathbf{k}\cdot\mathbf{w}_{\bm{\rho}}\neq 0$, and $\mathbf{w}_{\bm{\rho}}\parallel\bm{\mathcal{P}}_{\bm{\rho}}\parallel\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)\parallel\mathbf{k}$, the DF $\epsilon(k,\omega)$ is due only to the LO vibrations and thus is a longitudinal DF. In the absence of a test charge ($\sigma=0$), there is still a finite electric field $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ and potential $\phi(\bm{\rho},0)$ due to the LO vibrations, and therefore it follows from Eq. (85) that the LO modes are the solutions to

[TABLE]

From Eqs. (87) and (91) one indeed obtains the LO phonon dispersion $\omega_{l}(k)$ of expression (49). In bulk crystals the lattice DF is zero at the bulk LO mode frequency Born and Huang (1954). From expression (87) the static DF is $\epsilon_{0}(k)=1+2\pi\chi_{0}k$, while at high frequencies $\omega\gg\omega_{0}$ the vibrational contribution is negligible yielding the DF $\epsilon_{\infty}(k)=1+2\pi\chi_{e}k$, with $2\pi\chi_{0}$ and $2\pi\chi_{e}$ being the effective screening lengths Cudazzo et al. (2011); Berkelbach et al. (2013). Both static and high-frequency DFs have the same form as that deduced by Cudazzo et al. Cudazzo et al. (2011).

In terms of the three key quantities $\omega_{t}$, $\omega_{l}(k)$, $\epsilon_{\infty}(k)$ the lattice DF of the 2D crystal can be expressed in the form

[TABLE]

which is similar to the lattice DF of bulk polar crystals Born and Huang (1954), the difference being that here both the LO phonon frequency and DF are dependent on wavevector. In bulk crystals there is the 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 deduced from expression (92),

[TABLE]

with the three quantities dependent on wavevector.

The extended LST relation (93) connects the phonon frequencies to the two DFs; given the former, then the ratio of the latter is known, and vice versa. Furthermore, since the difference between the LO and TO frequencies is due entirely to the macroscopic field, the relation can be used to measure the phonon frequency change caused by the field. Similarly, for out-of-plane motion the susceptibilities $\chi_{0}^{\prime}$ and $\chi_{e}^{\prime}$ are related to phonon frequencies $\omega_{z}$ and $\omega_{0}^{\prime}$ via the frequency–susceptibility relation (81), which conveniently quantifies the effect of the macroscopic field on the ZO phonon frequency.

III Lattice dynamical properties: local field and polarizable ion effects

III.1 In-plane optical modes

The 2D clamped-ion susceptibility $\chi_{e}$, a key quantity entering the derived expressions for in-plane motion, can be obtained from first-principles calculation in four ways as follows. The 2D DF $\epsilon(q)$ was calculated within the random-phase approximation by a first-principles supercell approach for a number of 2D crystals such as h-BN Rasmussen et al. (2016), phosphorene and TMDs Berkelbach et al. (2013); Hüser et al. (2013); Rasmussen et al. (2016). In the first method, from the calculated DF the $\chi_{e}$ values of 2D TMDs were extracted in Ref.Berkelbach et al. (2013) by employing the relation $\epsilon(L)=1+4\pi\chi_{e}/L$ Cudazzo et al. (2011), where $L$ is the interlayer separation for a supercell containing two MLs of TMD, and $\epsilon(L)$ is the in-plane dielectric constant due to electronic polarization. For 2D MoS2, the obtained $\chi_{e}$=6.6 $\AA$ Berkelbach et al. (2013) is nearly equal to the susceptibility 6.5 $\AA$ extracted from bulk MoS2 simply by taking $L=c/2$ ($c$ is the lattice constant of bulk MoS2) whilst making $\epsilon(L)$ equal to the in-plane dielectric constant of bulk MoS2. When the experimental lattice constant $c$=6.66 $\AA$ Solozhenko et al. (1995) and in-plane high-frequency dielectric constant $\epsilon=4.95$ (Ref.Geick et al. (1966)) of bulk h-BN are put into the above relation, the susceptibility of 2D BN is estimated to be $\chi_{e}$=1.0 $\AA$. The 2D susceptibility $\chi_{e}$ can also be obtained by employing $\epsilon(q)=1+2\pi\chi_{e}q$ Cudazzo et al. (2011) as the slope of the $\epsilon(q)$ curve Hüser et al. (2013); Rasmussen et al. (2016) in the small-wavevector $q$ region equals the screening length $r_{eff}=2\pi\chi_{e}$. From Fig. 1 of Ref.Rasmussen et al. (2016) one finds a susceptibility $\chi_{e}$=5.9 $\AA$ for 2D MoS2 and $\chi_{e}$=0.8 $\AA$ for 2D hBN. The 2D susceptibility can also be calculated when knowing the screening length $r_{eff}$. For 2D hBN the calculated screening length is 7.64 $\AA$ Sohier et al. (2017), thus corresponding to a susceptibility $\chi_{e}$ of 1.216 $\AA$. In the fourth method, the (clamped ion) electronic energy in the unit cell is calculated for crystals in an external electric field Ferrero et al. (2008); Ferrabone et al. (2011), and then the first derivative of the electronic energy with respect to the field yields the induced dipole moment per cell, namely, the induced macroscopic polarization corresponding, for 2D hBN, to the last term in Eqs. (14) and (34), and the second derivative of the energy with respect to the field yields the electronic polarizability per cell. For ML hBN an electronic polarizability 4.591 $\AA^{3}$ per cell (Table IV of Ref.Ferrabone et al. (2011)) was obtained from a coupled perturbed Kohn-Sham (CPKS) calculation with the ab initio CRYSTAL09 code. Dividing this value by the area of the unit cell $s$ we find a 2D susceptibility $\chi_{e}$=0.848 $\AA$ which is very close to the value 0.8 $\AA$ of the second method. Similarly, using the static polarizability 7.111 $\AA^{3}$ per cell given in the same reference Ferrabone et al. (2011) we find a static dielectric susceptibility (i.e., including the lattice contribution) $\chi_{0}$=1.314 $\AA$. Unless otherwise stated, these 2D susceptibilities $\chi_{e}$ and $\chi_{0}$ from Ref.Ferrabone et al. (2011) are used in this study.

Through the three expressions (19), (47) and (74) we obtained, the three quantities of the 2D crystal, namely, the two macroscopic susceptibilities $\chi_{e}$, $\chi_{0}$ and the collective vibration frequency $\omega_{0}$ are related to the three microscopic quantities, $e_{a}$, $e_{B}$, $K$. Therefore, the values of the three microscopic quantities can be determined by using a set of three values $\chi_{e}$, $\chi_{0}$, $\omega_{0}$ calculated from first principles. Further, all the three mutually independent $a$-coefficients $a_{11}$, $a_{21}$ (or $a_{12}$), $a_{22}$ in the equation of motion (25) and polarization equation (14) can also be determined through Eqs. (16), (26a) and (75). Of the four microscopic quantities $\alpha_{1}$, $\alpha_{2}$, $e_{a}$, $K$ on which the $a$-coefficients originally depend [refer to Eqs. (15a), (15b), (26a), (22) and (17)], the polarizabilities of the constituent atoms of ML hBN $\alpha_{1}$, $\alpha_{2}$ are unknown. Therefore, the adoption of the three known quantities $\chi_{e}$, $\chi_{0}$, $\omega_{0}$ facilitates the use of the deduced equations for in-plane motion by circumventing the two unknowns $\alpha_{1}$ and $\alpha_{2}$. We note that, although the values of $\alpha_{1}$, $\alpha_{2}$ are unknown, their sum $\alpha_{1}+\alpha_{2}$, the atomic poarizability of the unit cell, is found to fall in an interval, expressed in terms of $\chi_{e}$ in inequality (103), which is obtained with a generalized Clausius-Mossotti relation connecting the microscopic polarizabilities $\alpha_{1}$, $\alpha_{2}$ and the macroscopic dielectric susceptibility $\chi_{e}$ (refer to Appendix A). Taking $\chi_{e}=0.85~{}\AA$ Ferrabone et al. (2011) yields $1.3084~{}\AA^{3}\leq\alpha_{1}+\alpha_{2}\leq 1.7530~{}\AA^{3}$, which are smaller than the total free-atom polarizability 4.0 $\AA^{3}$ (Appendix A). This interval will be used below to evaluate the LFEs on the LO phonon frequency and 2D dielectric susceptibilities (Table 2 below).

Using the susceptibilities $\chi_{e}=0.848~{}\AA$ and $\chi_{0}=1.314~{}\AA$ from CPKS calculation (Ref.Ferrabone et al. (2011)) and the phonon frequency $\omega_{0}$=1371 cm*-1* calculated from first principles by the same group (Ref.Erba et al. (2013)), we calculated the three microscopic quantities $e_{a}$, $e_{B}$, $K$ from Eqs. (19), (47) and (74). We also calculated the force constant due to LFEs $K_{e}$ [Eq. (22)] and the LO phonon group velocity $c_{l}$ [Eq. (67)] and present the result in Table 1 (upper row). In the theoretical study by Sohier et al. Sohier et al. (2017) the $\mathcal{S}$ parameter, related to $e_{B}$ by $\mathcal{S}=2\pi e_{B}^{2}/(\bar{m}s)$, is calculated by density-functional perturbation theory (DFPT). For 2D hBN the value $\mathcal{S}=8.4\times 10^{-2}$ eV{}^{2}\cdot$$\AA Sohier et al. (2017) so we find $e_{B}=2.71e$. In that study, the calculated frequency $\omega_{0}$ is 1387.2 cm*-1* (Table 1 of Ref.Sohier et al. (2017)), and recall that a susceptibility $\chi_{e}$=1.216 $\AA$ has been already calculated above through $r_{eff}$=7.64 $\AA$ Sohier et al. (2017). Now having these three values of $\chi_{e}$, $e_{B}$, $\omega_{0}$, again by using Eqs. (19), (47) and (74) we calculated the quantities $K_{e}$, $c_{l}$ as well as $\chi_{0}$, $e_{a}$, $K$, which are also listed in Table 1 (lower row) for comparison. First, we see two nearly equal values of Born charge, 2.70$e$ and 2.71$e$, the latter being also equal to another recent DFPT calculation Michel et al. (2017) as well as the $e_{B}$ value of the bulk hBN Ohba et al. (2001). According to Eq. (19), $e_{B}$ is greater than $e_{a}$ as $\chi_{e}>0$ due to EP of ions, and the numerical result shows that $e_{B}$ is four times larger than $e_{a}$. Also the two effective charge $e_{a}$ values are close to previous first-principles results 0.56$e$ Grad et al. (2003) and 0.43$e$ Topsakal et al. (2009) obtained respectively with Bader’s method and Löwdin’s analysis on charge transfer. Second, both susceptibilities $\chi_{e}$ and $\chi_{0}$ are larger in the lower row but the vibrational contributions to the static susceptibility namely the $\chi_{0}-\chi_{e}$ values are nearly equal in the two cases, which are 0.46 and 0.45 $\AA$ for the upper and lower rows, respectively. Third, the group velocity $c_{l}$ corresponds to the slope of the LO phonon dispersion curve at $\Gamma$; the two values of $c_{l}$ are very close with a difference of only 0.38% (Table 1) but they are one order of magnitude larger than the value 1.94 km/s calculated without taking into account EP of ions Michel and Verberck (2009). Fourth, in Ref.Michel and Verberck (2009) the force constant parameters are generated by making a percentage reduction of those of graphene, and from their table I we work out a spring force constant $K$=41.86 eV/$\AA^{2}$ ($K^{\prime}$=14.48 eV/$\AA^{2}$ for out-of-plane motion), which is 26% and 30% smaller than the force constant $K$ values 56.5 and 59.8 eV/$\AA^{2}$ of the present Table 1. Fifth, the intrinsic oscillator frequency $\omega_{0}$ is determined by $K-K_{e}$ [Eq. (26a)] rather than by $K$ alone, i.e., $\omega_{0}^{2}=(K-K_{e})/\bar{m}$; then the ratio $K_{e}/K$ expresses the percentage reduction due to the LFEs, which from Table 1 is $K_{e}/K$=24% and 29%. Sixth, the two quantities $c_{l}/\omega_{0}$ and $4\pi\chi_{e}$ form a factor of the $k^{2}$ term in the expansion of $\omega_{l}(k)$ [Eq. (66)], and from Table 1 the ration is $\frac{c_{l}}{\omega_{0}}/4\pi\chi_{e}$=13.5% and 9.3%, indicating the EP of ions makes the dominant contribution. Further, this result, together with the third point above, shows that neglecting ionic polarization such as RIM Michel and Verberck (2009) cannot properly describe the LO phonon dispersion. Below we shall use the first row of parameters to evaluate the influence of EP and LFEs on the lattice vibrations (as is presented in Table 2).

In Sec. II the electronic polarization of all the ions is included through $\alpha_{\kappa}$, and the model is referred to as PIM hereafter, in particular when comparing to the results calculated with the RIM ($\alpha_{\kappa}=0$). In the above calculations for Table 1, therefore both the EP and LFEs have been taken into account. We now look at what happens when EP or LFEs are not included. To do this, we compare various lattice-dynamical quantities such as $\chi_{e}$, $\chi_{0}$, $e_{B}$, $\omega_{t}$, $\omega_{l}(k)$, $c_{l}$ obtained with the RIM and PIM when the LFEs are neglected or taken into account, the expressions of which are listed in Table 2, for the four approaches in total. Given $\chi_{e}$, $\chi_{0}$, $\omega_{0}$, their expressions (19), (47) and (74), as noted in Table 2, are used to calculate the microscopic quantities such as $e_{a}$, $e_{B}$, $K$, $\alpha_{1}+\alpha_{2}$ (see Appendix A) needed in these approaches. Just below the expressions in the table the specific value-substituted LO phonon dispersion $\omega_{l}(k)$ and the values of the other five quantities are also listed for comparison: for the PIM including the LFEs (last row), the susceptibilities $\chi_{e}$, $\chi_{0}$ and frequency $\omega_{t}$ (equal to $\omega_{0}$) are input values taken from first-principles calculations Ferrabone et al. (2011); Erba et al. (2013), while the others were calculated using this set of $\chi_{e}$, $\chi_{0}$, $\omega_{0}$ values, as we did in the above calculations for Table 1; for the other three approaches, the quantities are transformed to depend on $e_{a}$ and/or $K$ as their expressions show, and are therefore calculated using the $e_{a}$ and $K$ parameters that are obtained with the same set of $\chi_{e}$, $\chi_{0}$, $\omega_{0}$ values (as was given in the first row of Table 1). Several points can be made by comparing the results obtained with the four approaches. First, there is no electronic polarization in the RIM and therefore $\chi_{e}=0$. In the PIM with no LFEs the electronic susceptibility of ML hBN is simply $\chi_{e}=(\alpha_{1}+\alpha_{2})/s$, and taking the interval we obtained above for $\alpha_{1}+\alpha_{2}$ gives $0.2417~{}\AA\leq\chi_{e}\leq 0.3239~{}\AA$, values that are even smaller than half the 0.85 $\AA$ obtained when the LFEs are taken into account. Second, in the RIM—for both cases of excluding and including the LFEs—and also in the PIM when neglecting LFEs, the Born charge reduces simply to the static charge $e_{a}=0.61e$, a value only 23% of the $e_{B}$ that we calculated with the PIM including the LFEs. Third, the ions in motion contribute $e_{B}^{2}/[s(K-e_{a}e_{B}Q_{1})]$ to the static susceptibility $\chi_{0}$, and with the RIM a smaller Born charge of 0.61$e$ causes this contribution to be only $\sim$0.018 $\AA$, thus underestimating the vibrational contribution by 96% compared to the 0.46 $\AA$ which was calculated with the PIM including the LFEs. Fourth, the LO phonon group velocity $c_{l}$ evaluated with the three simpler approaches, 1.61-1.67 km/s, is similar to the 1.94 km/s of a RIM result Michel and Verberck (2009) but is one order of magnitude smaller than the 37.24 km/s given by the PIM including the LFEs.

Furthermore, the LO phonon dispersion relations given by the four quantitative expressions of Table 2 are compared in Fig. 3, where wavevector $k$ is given in units of $\lvert\Gamma-K\rvert$ Sohier et al. (2017), the distance between the $\Gamma$ and $K$ points in the Brillouin zone. For the PIM without LFEs, the parameter $\gamma$ is simply taken to be the lower bound 3.8172 as the numerical calculation indicates that the use of the upper bound 5.1144, for instance, reduces the LO phonon frequencies by less than 0.08% (not shown), a change indiscernible to the dispersion curve. When LFEs are neglected, the two models PIM and RIM yield very close phonon frequencies that fall in a very narrow range from 1631 to 1653 cm*-1* (upper two lines). The two curves touch each other at the $\Gamma$ point with the same tangent and slope $c_{l}$=1.61 km/s; the curve obtained by the RIM displays a linear dispersion relation through the long-wavelength region $k\leq 0.15\lvert\Gamma-K\rvert$, while the PIM curve becomes flatter at larger wavevectors in the same region of $k$. When the LFEs are taken into account, the phonon frequencies are reduced significantly. For the dispersion calculated with the RIM (middle line), the dependence of frequency on wavevector remains linear, with a small slope similar to the case with no LFEs. For the dispersion curve obtained with the PIM (solid line), in contrast, owing to the EP of ions, a steep slope appears on the small wavevector side, corresponding to an substantially increased phonon group velocity 37.24 km/s, and as the wavevector increases the increase of the phonon frequency becomes slower (i.e. with a smaller $c_{l}$), deviating significantly from its linear dependence near the $\Gamma$ point. The LO phonon dispersion of Ref.Sohier et al. (2017), $\omega_{l}=[\omega_{0}^{2}+\mathcal{S}k/(1+r_{eff}k)]^{1/2}$, with first-principles calculated values $\omega_{0}$=1387.2 cm*-1*, $\mathcal{S}=8.4\times 10^{-2}$ eV{}^{2}\cdot$$\AA, and $r_{eff}$=7.64 $\AA$, is also shown (dotted line) and is very close to the dispersion curve obtained with the our PIM. According to Eq. (49), the upper limit of $\omega_{l}$ at a sufficiently large $k$ is $\omega_{M}$=1701.6 cm*-1*; at wavevector $k=0.15\lvert\Gamma-K\rvert$, the LO phonon frequency $\omega_{l}$ is 1569 cm*-1*, which is quite close to its limiting value (only 7.8% smaller). For wavevectors $k>0.02\lvert\Gamma-K\rvert$, as $\omega$ increases, the decrease of the group velocity with wavevector $k$ causes a rapid increase to the LO phonon DOS, as is clearly seen from Fig. 3 [also refer to the DOS expression (68)]. Overbending is a prominent feature of the numerically calculated LO phonon dispersion of ML hBN Miyamoto et al. (1995); Sánchez-Portal and Hernández (2002); Wirtz et al. (2003); Serrano et al. (2007); Michel and Verberck (2009); Topsakal et al. (2009). For instance, a maximum LO phonon frequency of 1533 cm*-1* (Ref.Wirtz et al. (2003)) or 1570 cm*-1* (Ref.Topsakal et al. (2009)) appears at a wavevector $\sim$0.30$\lvert\Gamma-K\rvert$. Evidently a positive slope $c_{l}$ is essential for the overbending of the dispersion curve $\omega_{l}(k)$. Therefore our result above indicates that due to the EP of ions, the increase of group velocity $c_{l}$ has enhanced the overbending substantially.

III.2 Out-of-plane optical modes

Values of phonon frequency $\omega_{z}$ Erba et al. (2013); Topsakal et al. (2009); Miyamoto et al. (1995); Wirtz et al. (2003), 2D electronic susceptibility $\chi_{e}^{\prime}$ and static susceptibility $\chi_{0}^{\prime}$ Ferrabone et al. (2011) associated with out-of-plane motion have been calculated from first principles. Having these then the intrinsic oscillator frequency $\omega_{0}^{\prime}$ is known through the relation (81). Given three quantities $\omega_{0}^{\prime}$, $\chi_{e}^{\prime}$, $\chi_{0}^{\prime}$, the three microscopic quantities namely the static effective charge $e_{a}^{\prime}$, the Born charge $e_{B}^{\prime}$ and the effective force constant $K^{\prime}$ follow from the three equations (36b), (38a) and (78). Further, all the three mutually independent $c$-coefficients $c_{11}$, $c_{12}$ (or $c_{21}$), $c_{22}$ in the lattice equations (37) and (34) are also determined by Eqs. (36b), (38a) and (79).

From Ref.Ferrabone et al. (2011) the 2D electronic susceptibility is given by $\chi_{e}^{\prime}=0.815\AA^{3}/s=0.151\AA$. Putting this into Eq. (36b) gives a negative ratio $e_{B}^{\prime}/e_{a}^{\prime}=-0.21$; that is, the positive B ions carry a negative Born charge while the negative N ions move with a positive Born charge. As a result, the polarization induced by ionic vibrations $e_{B}^{\prime}(\mathbf{u}_{1}-\mathbf{u}_{2})/s$ is antiparallel (parallel) to the displacements of the B (N) ions, and meanwhile the electric force exerted on the B (N) ions is antiparallel (parallel) to the macroscopic field. However the negative $e_{B}^{\prime}$ value has no influence on the out-of-plane optical (ZO) phonon frequency as $e_{B}^{\prime}$ enters the frequency expression via $e_{B}^{\prime 2}$ [Eq. (63)].

Using the first-principles calculated values $\chi_{e}^{\prime}=0.151~{}\AA$ and $\chi_{0}^{\prime}=0.164~{}\AA$ (Ref.Ferrabone et al. (2011)), we calculated $\omega_{0}^{\prime}$, $e_{a}^{\prime}$, $e_{B}^{\prime}$, $K^{\prime}$, $K_{e}^{\prime}$ for two values of frequency $\omega_{z}$ taken from previous first-principles calculations, 836 cm*-1* (first-principles perturbation result of Ref.Erba et al. (2013), by the same group of Ref.Ferrabone et al. (2011), and also the direct method result of Ref.Miyamoto et al. (1995)), 800 cm*-1* (DFPT value of Refs.Wirtz et al. (2003); Topsakal et al. (2009)), and for one experimental $\omega_{z}$ of 734 cm*-1* (Ref.Rokuta et al. (1997)) and also for a low frequency $\omega_{z}$=405 cm*-1* (see Table 3). A larger static effective charge $e_{a}^{\prime}$, greater than 1.0$e$, occurs for out-of-plane motion than in-plane motion (compare with in-plane effective charge $e_{a}$ values given in Table 1). Charge transfer is quite complicated in semiconductors where covalent bonds and ionic bonds coexist Phillips (1973). For 2D hBN there is a space allowing a more electron charge distribution outside the layer plane, probably causing a larger out-of-plane charge transfer from B to N. In fact, the calculation shows that a smaller $e_{a}^{\prime}$ equal to the in-plane effective charge of $0.61e$ corresponds to a very low phonon frequency 405 cm*-1* (last row of Table 3) that is unacceptably lower than the first-principles and experimental values. We see a negative Born charge $e_{B}^{\prime}$ (on B ions) one order of magnitude smaller than the in-plane Born charge $e_{B}$, while the force constant due to LFEs $K_{e}^{\prime}$ is one-third the effective spring force constant $K^{\prime}$, similar to the in-plane motion case.

To examine the EP of ions and LFEs on the out-of-plane motion, we look at the lattice-dynamical properties $\chi_{e}^{\prime}$, $\chi_{0}^{\prime}$, $e_{B}^{\prime}$, $\omega_{0}^{\prime}$, $\omega_{z}$ obtained with the RIM and PIM each including or excluding the LFEs, whose expressions are listed in Table 4 and also whose values obtained with a set of three first-principles values, $\chi_{e}^{\prime}=0.151$ $\AA$, $\chi_{0}^{\prime}=0.164$ $\AA$, $\omega_{z}=836$ cm*-1* (Refs.Ferrabone et al. (2011); Erba et al. (2013)), are given just below the corresponding expressions (the microscopic quantities $e_{a}^{\prime}$, $e_{B}^{\prime}$, $K^{\prime}$ used in the models were calculated using the same method as above for Table 3). For the adopted $\chi_{e}^{\prime}$ value 0.151 $\AA$ the atomic polarizability of the unit cell is found to fall in the interval, $1.9329~{}\AA^{3}\leq\alpha_{1}^{\prime}+\alpha_{2}^{\prime}\leq 2.5792~{}\AA^{3}$ (Appendix A). In the PIM without LFEs then the dielectric susceptibility of ML hBN is given by $\chi_{e}^{\prime}=(\alpha_{1}^{\prime}+\alpha_{2}^{\prime})/s$, thus yielding $0.3571~{}\AA\leq\chi_{e}^{\prime}\leq 0.4765~{}\AA$. When the LFEs are accounted for, the dielectric susceptibility $\chi_{e}^{\prime}$ is reduced significantly (by over 55%)—different from the in-plane motion case where the dielectric susceptibility is increased—and the Born charge on the B ions becomes negative with a magnitude that is only 21% of their static effective charge. In fact, there is a very small vibrational contribution (only 0.013 $\AA$) to the static susceptibility $\chi_{0}^{\prime}$ due to this small value of Born charge; in contrast, when LFEs are neglected, a large Born charge equal to the static effective charge leads to an overestimation of the vibrational contribution, which rises to a one order of magnitude higher value (0.196 $\AA$). The out-of-plane phonon frequency $\omega_{z}$ is 836 cm*-1* with the PIM including the LFEs, which is overestimated by over 39% when LFEs are not included. In the RIM there is no EP, i.e., $\alpha_{1}^{\prime}=\alpha_{2}^{\prime}=0$, which causes $\chi_{e}^{\prime}$ to vanish and consequently an extremely large out-of-plane phonon frequency, $\omega_{z}\rightarrow\infty$. In consequence the RIM fails to describe the out-of-plane vibrations.

IV Summary and Conclusions

We have studied long wavelength optical vibrations in ML hBN using two pairs of equations [Eqs. (39a), (39b), (40a) and (40b)] to describe the in-plane and out-of-plane lattice vibrations, respectively. These lattice equations, which have similar forms to Huang’s equations for bulk crystals, are deduced from a microscopic dipole lattice model taking into account the LFEs and EP self-consistently. The 2D Lorentz relation connecting the macroscopic and local fields, and the use of the areal polarization, a macroscopic quantity to describe the 2D dielectric polarization, are fundamental to deducing the lattice equations from the atomic theory. From the lattice equations the expressions for the areal energy density are obtained for the in-plane and out-of-plane lattice vibrations, respectively. The three mutually independent $a$- or $c$-coefficients of the equations are expressed in terms of a set of three quantities such as the 2D electronic and static susceptibilities and the intrinsic oscillator frequency, calculated from first principles, thus making the lattice equations very useful for calculating the lattice dynamical properties.

The lattice equations are solved simultaneously with the equation of electrostatics to deduce the optical modes of the TO, LO and ZO vibrations. The explicit expressions have been obtained for phonon frequency, macroscopic field, and also phonon group velocity and density of states. The frequency expressions are found to describe the dispersion relations of previous first-principles calculations very well: while the TO and ZO modes are dispersionless, the LO modes have dispersion with the frequency increasing with the wavevector. In particular, our LO phonon dispersion relation 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 a finite wavevector due to the long-range macroscopic field. It is also shown that the finite ZO phonon frequency exists owing to the electronic polarization of ions, without which the frequency becomes infinitely large in the rigid ion model.

From the lattice equations the frequency-dependent dielectric susceptibilities are deduced for in-plane and out-of-plane lattice motion. By considering the response of the lattice to a test charge with a charge distribution asymmetric or symmetric with respect to the ML, a 2D lattice dielectric function $\epsilon(k,\omega)$ is derived [expression (87)], which is shown to be due solely to the LO vibrations, and the driven LO and ZO vibrations are further discussed. It is also shown that such a 2D longitudinal DF allows the LO phonon dispersion [expression (49)] to be rederived from $\epsilon(k,\omega)=0$, similar to the case of bulk crystals. The 2D LST relation (93) and frequency–susceptibility relation (81) are obtained for in-plane and out-of-plane motion, respectively, connecting the phonon frequencies to the 2D dielectric functions or susceptibilities, which are very useful for quantifying the effects of the macroscopic field on the phonon frequencies.

The analytical expressions have been applied to study the lattice dynamical properties of ML hBN. For the in-plane motion, two sets of three quantities from two independent first-principles calculations, one set of $\chi_{e}$, $\chi_{0}$, $\omega_{0}$ from Ref.Ferrabone et al. (2011) and the other set of $\chi_{e}$, $e_{B}$, $\omega_{0}$ from Ref.Sohier et al. (2017), are used as parameters and very close results are obtained for the calculated properties such as the LO phonon dispersion relation, the effective spring force constant and the vibrational contribution to the static susceptibility. To evaluate the LFEs the unit-cell atomic polarizability is used, which, given a dielectric susceptibility, is found to be limited in an interval using the 2D Clausius-Mossotti relation. The LFEs and electronic polarization of ions should be included simultaneously, but otherwise neglecting either or both causes large discrepancies to the calculated properties: the phonon frequency at $\Gamma$ is overestimated by 15%-19%, whereas the Born charge, the LO phonon group velocity and the vibrational contribution to the static susceptibility are underestimated by 77%, 96%, 96%, respectively. With no ionic EP or LFEs, the LO modes display very small linear dispersion, almost flat in the long wavelength region, which is distinct from the LO phonon dispersion calculated after accounting for both EP and LFEs. For out-of-plane motion, using $\chi_{e}^{\prime}$, $\chi_{0}^{\prime}$, $\omega_{z}$ calculated mainly from Refs.Ferrabone et al. (2011); Erba et al. (2013) as parameters, the static effective charge is found to be different from that for in-plane motion reflecting anisotropy of the 3D charge density distribution. Further, the positive B ions carry a negative Born charge while the negative N ions move with a positive Born charge. Similar to the in-plane motion case, the RIM can not properly describe the out-of-plane vibrations, which, for instance, yields a four times larger Born charge and infinitely large phonon frequency. When the EP of ions is included, the LFEs have significantly reduced the 2D dielectric susceptibility and the Born charge (by 60% and 80%, respectively), different from the in-plane motion case where the dielectric susceptibility and Born charge are both increased substantially (by two and three times, respectively).

Acknowledgements.

We acknowledge support from the Natural Science Research Funds (Nos. 419080500175 & 419080500260) of Jilin University.

Appendix A The generalized Clausius-Mossotti relations and intervals of the total atomic polarizability in ML hBN

In bulk diatomic crystals the electronic susceptibility $\chi$ relates to the atomic polarizabilities of constituent atoms $\alpha_{1}$ and $\alpha_{2}$ by the Clausius-Mossotti relation Clausius (1879); Mossotti (1850),

[TABLE]

There is no such a simple relation for 2D hBN. Substituting Eqs. (11) and (16) into Eq. (15b), $\chi_{e}$, the in-plane electronic susceptibility of ML hBN, is given by Mikhailov (2013); Dell’Anna and Merano (2016)

[TABLE]

If $Q_{0}=Q_{1}$, then the 2D susceptibility of the dielectric $\chi_{e}$ relates to the total atomic polarizability of the unit cell $\alpha_{1}+\alpha_{2}$ by a simple relation,

[TABLE]

which is similar to the Clausius-Mossotti relation (94) for bulk materials. In fact in ML hBN $Q_{0}\neq Q_{1}$ (Sec. II) so the sum $\alpha_{1}+\alpha_{2}$ can not be determined from the generalized Clausius-Mossotti relation (95); adding the Born charge expression (17) does not help because $e_{B}$ and $\chi_{e}$ are related through Eq. (19) and not independent of each other.

Interestingly, when knowing $\chi_{e}$ the upper and lower bounds of the unit-cell atomic polarizability $\alpha_{1}+\alpha_{2}$ can be determined from expression (95). Let $\alpha=\alpha_{1}+\alpha_{2}$. Then Eq. (95) can be transformed into a quadratic equation in variable $\alpha_{1}$,

[TABLE]

with the discriminant $\Delta_{\alpha}$,

[TABLE]

Eq. (97) has two real roots only when $\Delta_{\alpha}\geq 0$. Therefore it follows from the $\Delta_{\alpha}$ expression (98) that either

[TABLE]

or

[TABLE]

For ML hBN ($Q_{1}>Q_{0}$) the latter inequality leads to only negative $\alpha$ and is dropped. Now one finds that $\alpha_{1}$ and $\alpha_{2}$ are given by

[TABLE]

if $\alpha_{1}$ takes the plus square root then $\alpha_{2}$ takes the minus square root, and vice versa. For $\alpha_{1}$ and $\alpha_{2}$ to be positive, one further requires that the above square root is not greater than $\alpha$ and obtains

[TABLE]

Combining the inequalities (99) and (102) then restricts the atomic polarizability in the unit cell for in-plane polarization to the following interval,

[TABLE]

When the 2D susceptibility $\chi_{e}=0.85~{}\AA$ Ferrabone et al. (2011) is used, for instance, one finds $1.3084~{}\AA^{3}\leq\alpha_{1}+\alpha_{2}\leq 1.7530~{}\AA^{3}$, which is only a narrow range of 0.44 $\AA^{3}$.

For the vertical polarization (parallel to $\mathbf{e}_{z}$), the electronic susceptibility $\chi_{e}^{\prime}$ can be expressed in terms of the atomic polarizabilities $\alpha_{1}^{\prime}$, $\alpha_{2}^{\prime}$ as (refer to Sec. II above and also Ref.Mikhailov (2013))

[TABLE]

Following the preceding derivation we find that the atomic polarizability of the unit cell falls in the interval,

[TABLE]

Taking the 2D susceptibility $\chi_{e}^{\prime}=0.151~{}\AA$ Ferrabone et al. (2011), one finds $1.9329~{}\AA^{3}\leq\alpha_{1}^{\prime}+\alpha_{2}^{\prime}\leq 2.5792~{}\AA^{3}$. Mathematically, the solution $\alpha_{1}^{\prime}+\alpha_{2}^{\prime}\geq 1/(\frac{1}{\chi_{e}^{\prime}s}-Q_{1}-Q_{0})$ is also permitted, but yielding $\alpha_{1}^{\prime}+\alpha_{2}^{\prime}\geq 7.3~{}\AA^{3}$ for ML hBN, which are much larger than the total polarizability 4.0 $\AA^{3}$ of the free atoms B and N (the polarizabilities of free atoms B and N are 3.038 and 1.097 $\AA^{3}$, respectively Schwerdtfeger (2006)). Using the linear volume-polarizability relationship Brinck et al. (1993), the polarizability of a constituent atom $\alpha_{i}$ in a crystal can be estimated by $\alpha_{i}=(v_{i}/v_{i}^{0})\alpha_{i}^{0}$ Tkatchenko and Scheffler (2009); Hod (2012), where $v_{i}$ is the effective volume of the atom in the crystal, and $v_{i}^{0}$ and $\alpha_{i}^{0}$ are the free-atom volume and polarizability, respectively. Numerical calculations have shown that the polarizabilities of the B and N atoms in bulk h-BN are reduced compared to their free-atom polarizabilities Tkatchenko and Scheffler (2009); Hod (2012). Anyhow the values $\alpha_{1}^{\prime}+\alpha_{2}^{\prime}\geq 7.3~{}\AA^{3}$ are unacceptably large to be considered.

Compared to the their free-atom polarizability 4.0 $\AA^{3}$, the atoms B and N have a significantly reduced total polarizability after forming the 2D crystal hBN. There is a slightly larger atomic polarizability in the vertical direction than in the layer plane, i.e., $\alpha_{1}^{\prime}+\alpha_{2}^{\prime}>\alpha_{1}+\alpha_{2}$ (with a 0.6 $\AA^{3}$ difference in upper bound), as the constituent B and N atoms are close packed on the plane causing the electron cloud to be more easily deformed and polarized along the vertical direction. If LFEs are neglected then the dielectric susceptibilities become $\chi_{e}=(\alpha_{1}+\alpha_{2})/s$, and $\chi_{e}^{\prime}=(\alpha_{1}^{\prime}+\alpha_{2}^{\prime})/s$, leading to $\chi_{e}^{\prime}>\chi_{e}$. Therefore the first-principles result $\chi_{e}>\chi_{e}^{\prime}$ as shown in Ref.Ferrabone et al. (2011) is attributed to the LFEs.

We note that $\chi_{e}$ and $\chi_{e}^{\prime}$ are susceptibilities due to polarizable atoms or ions, and thus they also occur in non-ionic 2D crystals such as graphene. For graphene expressions (95) and (104) for $\chi_{e}$ and $\chi_{e}^{\prime}$ are still applicable but become simpler considering there are two C atoms in a unit cell, namely $\alpha_{1}=\alpha_{2}$ and $\alpha_{1}^{\prime}=\alpha_{2}^{\prime}$.

Appendix B Proof of $a_{12}=a_{21}$ and $c_{12}=c_{21}$ from macroscopic theory and areal energy density in 2D hBN

We first consider in-plane motion. Let the ions have a displacement in the layer plane, $\mathbf{w}(\bm{\rho})=\mathbf{w}_{\bm{\rho}}(\bm{\rho})$, in the presence of a point charge $e_{\nu}$ fixed at ($\bm{\rho}_{0}$,0), and let us find the electric field in the monolayer (in-plane component), $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$. The polarization $\mathbf{P}$ is $\mathbf{P(\mathbf{r})}=[a_{21}\mathbf{w}(\bm{\rho})+a_{22}\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)+c_{22}\mathbf{E}_{z}(\bm{\rho},0)]\delta(z)$ upon using Eqs. (39b) and (40b), and the Poisson equation is given by

[TABLE]

Expanding $\phi(\bm{\rho},z)$, $\mathbf{w}(\bm{\rho})$ and $\delta(\bm{\rho}-\bm{\rho}_{0})$,

[TABLE]

[TABLE]

[TABLE]

$A$ being the sample area, Eq. (106) is then transformed to

[TABLE]

Expanding $\varphi_{\mathbf{k}}(z)$

[TABLE]

and $\delta(z)$ as Eq. (2), one finds $\hat{\varphi}_{\mathbf{k}}(q)$ from Eq. (110),

[TABLE]

We note that from the above expansions one has reality conditions $\varphi_{-\mathbf{k}}(z)=\varphi^{*}_{\mathbf{k}}(z)$, $\mathbf{w}_{-\mathbf{k}}=\mathbf{w}^{*}_{\mathbf{k}}$ and $\hat{\varphi}_{-\mathbf{k}}(-q)=\hat{\varphi}^{*}_{\mathbf{k}}(q)$ considering that quantities such as $\phi(\bm{\rho},z)$, $\mathbf{w}(\bm{\rho})$ are real.

Taking $z=0$ in Eq. (111) and integrating $\hat{\varphi}_{\mathbf{k}}(q)$ over $q$ one finds

[TABLE]

and then the in-plane component of the electric field in the monolayer,

[TABLE]

where the latter sum represents the field due to the point charge (including dielectric effects through $a_{22}$), which is zero at the charge site. Eq. (114) will be used below to calculate the work required to displace the ions or the charge.

We have obtained Eq. (27), namely, $a_{12}=a_{21}$, from the microscopic dipole lattice model, and here we show that from the viewpoint of macroscopic theory it follows from the principle of energy conservation. Place the point charge $e_{\nu}$ at the origin, whilst setting the ions in the configuration $\mathbf{w}(\bm{\rho})=0$, and consider the following cycle Born and Huang (1954): (a) keeping the charge at the origin, displace the ions horizontally into an irrotational configuration $\mathbf{w}(\bm{\rho})=\nabla_{\bm{\rho}}\psi(\bm{\rho})$, according to

[TABLE]

by increasing $\xi$ from 0 to 1, $\psi(\bm{\rho})$ being a nonzero arbitrary scalar field; (b) keeping the ions at $\mathbf{w}(\bm{\rho})=\nabla_{\bm{\rho}}\psi(\bm{\rho})$, displace the charge to $\Delta\bm{\rho}$; (c) fixing the charge at $\Delta\bm{\rho}$, reverse process (a), i.e., by decreasing $\xi$ from 1 to 0 according to Eq. (115); (d) move the charge back to the origin to complete the cycle.

$-[a_{11}\mathbf{w}(\bm{\rho})+a_{12}\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)]\cdot\Delta\mathbf{w}(\bm{\rho})$ is the work per unit area required to change $\mathbf{w}(\bm{\rho})$ to $\mathbf{w}(\bm{\rho})+\Delta\mathbf{w}(\bm{\rho})$, and total work expended on the ionic system for the configuration change is

[TABLE]

From Eq. (115) $\Delta\mathbf{w}(\bm{\rho})=\Delta\xi\nabla_{\bm{\rho}}\psi(\bm{\rho})$ follows, and then insert expansion $\psi(\bm{\rho})=\sum_{\mathbf{k}}\psi_{\mathbf{k}}e^{i\mathbf{k}\cdot\bm{\rho}}$ into $\mathbf{w}(\bm{\rho})$ [Eq. (115)] and $\Delta\mathbf{w}(\bm{\rho})$. Comparing the former expansion with Eq. (108) one finds $\mathbf{w}_{\mathbf{k}}$ for Eq. (114), $\mathbf{w}_{\mathbf{k}}=\xi i\mathbf{k}\psi_{\mathbf{k}}$. Inserting these expansions of $\mathbf{w}(\bm{\rho})$ and $\Delta\mathbf{w}(\bm{\rho})$ together with Eq. (114) into the above expression (116), and integrating over $\xi$ from 0 to 1, one obtains the work done during process (a),

[TABLE]

The field acting on charge $e_{\nu}$ during process (b) is given by setting $\bm{\rho}=0$ (as $\Delta\bm{\rho}$ is small) and $\mathbf{w}_{\mathbf{k}}=i\mathbf{k}\psi_{\mathbf{k}}$ in the first term of $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$ [Eq. (114)], and the work expended during the process is

[TABLE]

Process (c) is the reverse of process (a) except for the altered position of charge $e_{\nu}$. Therefore after reversing the sign of expression (117) and multiplying its latter summands by a factor $e^{-i\mathbf{k}\cdot\Delta\bm{\rho}}$ we find the work expended during the process,

[TABLE]

During process (d) $\mathbf{w}(\bm{\rho})=0$, and according to Eq. (114) no field acts on charge $e_{\nu}$, and therefore no work is needed to restore the charge to the origin, $W_{d}=0$.

As $W_{a}+W_{b}+W_{c}+W_{d}=0$, one finds

[TABLE]

Upon making $e^{-i\mathbf{k}\cdot\Delta\bm{\rho}}=1-i\mathbf{k}\cdot\Delta\bm{\rho}$, and changing the index of summation $\mathbf{k}$ to $-\mathbf{k}$ for the former summation, one immediately obtains $a_{12}=a_{21}$.

The relation $a_{12}=a_{21}$ allows us to define 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})$ and $\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)$,

[TABLE]

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

Now we consider out-of-plane motion. Let the ions have a displacement perpendicular to the layer plane, $\mathbf{w}(\bm{\rho})=\mathbf{w}_{z}(\bm{\rho})$, in the presence of a point charge $e_{\nu}$ placed at (0,$z_{0}$), and let us find the $z$-component of the field, $\mathbf{E}_{z}(\bm{\rho},z)$. Now the polarization $\mathbf{P}$ is $\mathbf{P(\mathbf{r})}=[c_{21}\mathbf{w}(\bm{\rho})+c_{22}\mathbf{E}_{z}(\bm{\rho},0)+a_{22}\mathbf{E}_{\bm{\rho}}(\bm{\rho},0)]\delta(z)$ upon using Eqs. (39b) and (40b)], and the Poisson equation is given by

[TABLE]

Expanding $\delta(z)$ [Eq. (2)] and the other spatially varying quantities exactly as in the in-plane motion case [Eqs. (107), (108), (109) and (111)], one finds from the Poisson equation $\hat{\varphi}_{\mathbf{k}}(q)$,

[TABLE]

Taking $z=0$ in Eq. (111) and integrating $\hat{\varphi}_{\mathbf{k}}(q)$ over $q$ one finds

[TABLE]

Differentiating Eq. (111) with respect to $z$ and substituting Eq. (123) for $\hat{\varphi}_{\mathbf{k}}(q)$ and then integrating over $q$, one obtains

[TABLE]

To find $\varphi^{\prime}_{\mathbf{k}}(0)$ we approach $\delta(z)$ with $\delta_{\varepsilon}(z)$ where $\varepsilon$ is a small thickness Michel and Verberck (2009); Sohier et al. (2017) and let $\Lambda=\delta_{\varepsilon}(0)$. Thus from Eq. (125) one finds

[TABLE]

Having $\varphi_{\mathbf{k}}(0)$ and $\varphi^{\prime}_{\mathbf{k}}(0)$ [Eqs. (124) and (126)], now insert them into Eq. (125) for $\varphi^{\prime}_{\mathbf{k}}(z)$, from which one obtains $\mathbf{E}_{z}(\bm{\rho},z)$ via $\mathbf{E}_{z}(\bm{\rho},z)=-\mathbf{e}_{z}\sum_{\mathbf{k}}\varphi^{\prime}_{\mathbf{k}}(z)e^{i\mathbf{k}\cdot\bm{\rho}}$,

[TABLE]

Let $z_{0}>0$ in what follows. Taking $z=0$ yields the $z$-component of the field in the ML,

[TABLE]

while taking $z=z_{0}$ and $\bm{\rho}=0$ one finds the $z$-component of the field acting on the point charge $e_{\nu}$,

[TABLE]

where the sum of the terms containing $e_{\nu}$ ($\mathbf{w}_{\mathbf{k}}$) represents the field due to the point charge (ionic displacements), denoted by $\mathbf{E}_{z,e_{\nu}}(0,z_{0})$ ($\mathbf{E}_{z,w}(0,z_{0})$) for simplicity.

Now we show that $c_{12}=c_{21}$ [Eq. (38b)] follows from the principle of energy conservation. Place the point charge $e_{\nu}$ above the ML at a point $P$ ($\bm{\rho}_{0}=0$,$z_{0}$), while keeping the ions in the configuration $\mathbf{w}(\bm{\rho})=0$, and consider the following cycle: (a) keeping the charge at the point $P$, displace the ions vertically into the configuration $\mathbf{w}(\bm{\rho})=\bm{\psi}(\bm{\rho})=\psi(\bm{\rho})\mathbf{e}_{z}$, according to

[TABLE]

by increasing $\xi$ from 0 to 1; (b) keeping the ions at $\mathbf{w}(\bm{\rho})=\bm{\psi}(\bm{\rho})$, displace the charge vertically by a small $\Delta z$ to point $P^{\prime}$ ($\bm{\rho}_{0}=0$,$z_{0}+\Delta z$); (c) fixing the charge at the point $P^{\prime}$, reverse process (a), i.e., by reducing $\xi$ from 1 to 0 according to Eq. (130); (d) move the charge back to point $P$ to complete the cycle.

The work per unit area required to change $\mathbf{w}(\bm{\rho})$ to $\mathbf{w}(\bm{\rho})+\Delta\mathbf{w}(\bm{\rho})$ is $-[c_{11}\mathbf{w}(\bm{\rho})+c_{12}\mathbf{E}_{z}(\bm{\rho},0)]\cdot\Delta\mathbf{w}(\bm{\rho})$, and thus total work expended on the ionic system for the configuration change is

[TABLE]

From Eq. (130) one has $\Delta\mathbf{w}(\bm{\rho})=\Delta\xi\psi(\bm{\rho})\mathbf{e}_{z}$. Expanding $\mathbf{w}(\bm{\rho})$ and $\psi(\bm{\rho})$ as in-plane motion above, one finds $\mathbf{w}_{\mathbf{k}}$ in Eq. (128) is given by $\mathbf{w}_{\mathbf{k}}=\xi\bm{\psi}_{\mathbf{k}}$. Substituting these expansions in terms of $\psi_{\mathbf{k}}$ together with Eq. (128) into the above expression (131), and integrating over $\xi$ from 0 to 1, one obtains for the work expended during process (a),

[TABLE]

During process (b) the field acting on charge $e_{\nu}$ due to the ionic displacements is $\mathbf{E}_{z,w}(0,z_{0})$, namely, the sum of the terms containing $\mathbf{w}_{\mathbf{k}}$ of Eq. (129), with $\mathbf{w}_{\mathbf{k}}=\bm{\psi}_{\mathbf{k}}$, and the work expended for this field of the ionic displacements is

[TABLE]

The field acting on charge $e_{\nu}$ due to the charge itself is $\mathbf{E}_{z,e_{\nu}}(0,z_{0})$, i.e., the sum of the terms containing $e_{\nu}$ of Eq. (129), and the work done for this field can be simply written as $W_{b,e_{\nu}}=-e_{\nu}E_{z,e_{\nu}}(0,z_{0})\Delta z$.

Process (c) is the reverse of process (a) except for the displacement of charge $e_{\nu}$. Therefore upon reversing the sign of expression (132) and multiplying its latter summands by a factor $e^{-k\Delta z}$, we readily find the work expended during the process,

[TABLE]

During process (d) $\mathbf{w}(\bm{\rho})=0$, and according to Eq. (129) no field is associated with displacement $\mathbf{w}(\bm{\rho})$, and accordingly no work is needed for the displacement part, $W_{d,w}=0$. For the field due to the point charge itself, the work done is given by reversing the sign of $W_{b,e_{\nu}}$, i.e., $W_{d,e_{\nu}}=e_{\nu}E_{z,e_{\nu}}(0,z_{0})\Delta z$.

As $W_{a}+W_{b,w}+W_{b,e_{\nu}}+W_{c}+W_{d,w}+W_{d,e_{\nu}}=0$, one finds

[TABLE]

Using $e^{-k\Delta z}=1-k\Delta z$, and changing the index of summation $\mathbf{k}$ to $-\mathbf{k}$ for the former summation, one readily finds $c_{12}=c_{21}$.

Having $c_{12}=c_{21}$, now an areal energy density associated with the out-of-plane optical vibrations can be introduced,

[TABLE]

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