Anisotropic hybrid excitation modes in monolayer and double-layer phosphorene on polar substrates

TL;DR

This study explores the anisotropic hybrid plasmon-phonon modes in monolayer and double-layer phosphorene on various polar substrates, revealing substrate effects and tunability through layer misalignment and separation.

Contribution

It provides analytical and numerical insights into hybrid mode dispersion relations in phosphorene systems, highlighting substrate influence and tunability mechanisms.

Findings

Strong anisotropy in hybrid modes along light-mass directions

Substrate type significantly affects dispersion relations

Layer misalignment and separation enable mode engineering

Abstract

We investigate the anisotropic hybrid plasmon-SO phonon dispersion relations in monolayer and double-layer phosphorene systems located on the polar substrates, such as SiO2, h-BN and Al2O3. We calculate these hybrid modes with using the dynamical dielectric function in the RPA by considering the electron-electron interaction and long-range electric field generated by the substrate SO phonons via Frohlich interaction. In the long-wavelength limit, we obtain some analytical expressions for the hybrid plasmon-SO phonon dispersion relations which represent the behavior of these modes akin to the modes obtaining from the loss function. Our results indicate a strong anisotropy in plasmon-SO phonon modes, whereas they are stronger along the light-mass direction in our heterostructures. Furthermore, we find that the type of substrate has a significant effect on the dispersion relations of the coupled modes. Also, by tuning the misalignment and separation between layers in double-layer phosphorene on polar substrates, we can engineer the hybrid modes.

Figures

Tables

Table 1. Table 1: Physical parameters of selected polar materials from Ref Mogulkoc  et al. ( 2016a ) .

	SiO2	h-BN	Al2O3
$\hslash {\omega }_{SO}^1$(meV)	60	101	55
$\hslash {\omega }_{SO}^2$(meV)	146	195	94
${ϵ}_0$	3.9	5.1	12.5
${ϵ}_{\infty }$	2.5	4.1	3.2
$\alpha$	0.2	0.132	0.525

Equations

$$H_0=\left[\begin{array}{cc}{E}_c+{\eta }_c{k}_x^2+{\nu }_c{k}_y^2& \gamma k_x+\beta k_y^2\\ \gamma k_x+\beta k_y^2& E_{\nu }-{\eta }_{\nu }{k}_x^2-{\nu }_{\nu }{k}_y^2\end{array}\right],$$

$$H_{el-el}=\frac{1}{2}\sum _{ij}\sum _{\mathbf{kqp}}{v}_{ij}(q)a_{\mathbf{k}+\mathbf{q},i}^{†}{a}_{\mathbf{p}-\mathbf{q},j}^{†}{a}_{\mathbf{p},j}{a}_{\mathbf{k},i},$$

$$H_{el-ph}=\sum _i\sum _{\lambda }\sum _{\mathbf{kq}}[M_0^{\lambda }((q)]a_{\mathbf{k}+\mathbf{q},i}^{†}{a}_{\mathbf{k},i}(b_{\mathbf{q}\lambda }+b_{-\mathbf{q}\lambda }^{†}),$$

$$V_{ij}^{sc}(\omega ,\mathbf{q})=\frac{v_{ij}(q)}{det|{ϵ}_{ij}(\omega ,\mathbf{q})|},$$

$$\displaystyle M_{0}^{\lambda}(q)=\big{[}v(q)\alpha\frac{\omega_{SO}^{\lambda}}{2}e^{-2qz}\big{]}^{1/2}\ ,$$

$$\displaystyle\alpha=\epsilon_{\infty}\big{[}\frac{1}{\epsilon_{\infty}+1}-\frac{1}{\epsilon_{0}+1}\big{]}\ .$$

$${ϵ}_{ij}(\omega ,\mathbf{q})={\delta }_{ij}-U_{ij}(\omega ,q){\Pi }_i(\omega ,\mathbf{q}),$$

$$U_{ij}(\omega ,q)=U_0(\omega ,q)e^{-q{d}_{ij}(1-{\delta }_{ij})},$$

$$v_{ph}(\omega ,q)=\sum _{\lambda }[M_0^{\lambda }(q){]}^2{D}_0^{\lambda }(\omega ),$$

$$D_0^{\lambda }(\omega )=\frac{2{\omega }_{SO}^{\lambda }}{{\omega }^2-({\omega }_{SO}^{\lambda }{)}^2}.$$

$${\Pi }_j(\omega ,\mathbf{q})=\frac{g_s}{\nu }\sum _{\mathbf{k}}\frac{f^j(E_{\mathbf{q}})-f^j(E_{\mathbf{k}+\mathbf{q}})}{E_{\mathbf{q}}-E_{\mathbf{k}+\mathbf{q}}+\hslash \omega +i\eta },$$

$$E_{\mathbf{k}}=\frac{{\hslash }^2}{2}(\frac{k_x^2}{m_x}+\frac{k_y^2}{m_y}),$$

$$\displaystyle\frac{\Pi_{i}(\omega,\mathbf{q})}{g_{2d}}=(\frac{1}{Q_{i}})\Bigg{[}$$

$$-$$

$$\displaystyle R_{i}(\theta)=\bigg{(}\frac{\cos^{2}(\theta-\tau_{i})}{m_{x}}+\frac{\sin^{2}(\theta-\tau_{i})}{m_{y}}\bigg{)}\ .$$

$$\frac{{\Pi }_i(\omega ,q,\theta )}{g_{2d}}\approx R_i(\theta )E_F\frac{q^2}{{\omega }^2}.$$

$${ϵ}_i(\omega ,\mathbf{q})=1-\frac{2\pi e^2}{{ϵ}_{\infty }q}{\Pi }_i(\omega ,\mathbf{q})+\sum _{\lambda }\frac{\alpha e^{-2qz}}{1-\alpha e^{-2qz}-{\omega }^2/({\omega }_{so}^{\lambda }{)}^2}.$$

$$\omega_{(+)}^{\lambda}(q,\theta)=\omega^{\lambda}_{so}\big{(}1+\alpha e^{-2qz}\frac{\omega^{2}_{pl}(q,\theta)}{(\omega_{so}^{\lambda})^{2}}\big{)}\ .$$

$$\omega_{(-)}(q,\theta)=\omega_{pl}(q,\theta)\big{(}1-\frac{\alpha}{2}e^{-2qz}\big{)}\ .$$

$$\hslash {\omega }_{(SPE)}^{\pm }(q,\theta )=\frac{R_i(\theta ){\hslash }^2}{2}\pm q{\nu }_F\sqrt{R_i(\theta )m_d}.$$

$${\omega }_{ac}(q,\theta )=2q\sqrt{\frac{n\pi e^2{d}_{12}}{{ϵ}_{\infty }}\frac{R_1(\theta )R_2(\theta )}{R_1(\theta )+R_2(\theta )}}.$$

$$\omega_{op}(q,\theta)=\sqrt{\frac{2n\pi e^{2}q}{\epsilon_{\infty}}\big{(}R_{1}(\theta)+R_{2}(\theta)\big{)}}\ .$$

$${\omega }_{ac(+)}^{\lambda }(q,\theta )={\omega }_{SO}^{\lambda }\sqrt{1+\frac{4n\pi d_{12}{e}^2{q}^2\alpha e^{-2qz}}{({\omega }_{SO}^{\lambda }{)}^2{ϵ}_{\infty }}\frac{R_1(\theta )R_2(\theta )}{R_1(\theta )+R_2(\theta )}}.$$

$${\omega }_{ac(-)}(q,\theta )=2q\sqrt{\frac{n\pi e^2{d}_{12}(1-\alpha e^{-2qz})}{{ϵ}_{\infty }}\frac{R_1(\theta )R_2(\theta )}{R_1(\theta )+R_2(\theta )}}.$$

$$\omega^{\lambda}_{op(+)}(q,\theta)=\omega^{\lambda}_{SO}\sqrt{1+\frac{2n\pi e^{2}q\alpha e^{-2qz}}{(\omega^{\lambda}_{SO})^{2}\epsilon_{\infty}}\big{(}R_{1}(\theta)+R_{2}(\theta)\big{)}}\ .$$

$$\omega_{op(-)}(q,\theta)=\sqrt{\frac{2n\pi e^{2}q(1-\alpha e^{-2qz})}{\epsilon_{\infty}}\big{(}R_{1}(\theta)+R_{2}(\theta)\big{)}}\ .$$

$$\displaystyle\frac{\Pi(\omega,q,\theta)}{g_{2d}}=\frac{1}{2}\bigg{[}-2$$

$$\displaystyle+\sqrt{(1+\frac{2\omega}{k_{F}\nu_{F}Q^{2}})^{2}-\frac{4}{Q^{2}}}\quad\bigg{]}$$

$$\displaystyle\frac{\Pi(\omega,q,\theta)}{g_{2d}}=\frac{1}{2}\bigg{[}-2$$

$$\displaystyle+(1+\frac{2\omega}{k_{F}\nu_{F}Q^{2}})\big{(}1-\frac{2/Q^{2}}{(1+\frac{2\omega}{k_{F}\nu_{F}Q^{2}})^{2}}\big{)}\bigg{]}$$

$$\displaystyle\frac{\Pi(\omega,q\to 0,\theta)}{g_{2d}}\approx\frac{-2}{Q^{2}}\bigg{(}\frac{1}{1-\frac{4\omega^{2}}{(k_{F}\nu_{F}Q^{2})^{2}}}\bigg{)}$$

$$ϵ(\omega ,q,\theta )=1-\frac{{\omega }_{pl}^2(q,\theta )}{{\omega }^2(q,\theta )}+\sum _{\lambda }\frac{\alpha e^{-2qz}}{1-\alpha e^{-2qz}-{\omega }^2(q,\theta )/({\omega }_{so}^{\lambda }{)}^2}$$

Full text

I Introduction

Phosphorene, a monolayer of black phosphorus (BP), has recently attracted special attention among two dimensional materials (2DMs) due to its unique highly anisotropic electronic and optical properties Prishchenko et al. (2017); Xia et al. (2014); Liu and Chang (2015); Aierken et al. (2015); Padilha et al. (2015); Saberi-Pouya et al. (2016); Zhu et al. (2016); Ostahie and Aldea (2016); Sa et al. (2015); Liu and Chang (2015); Constantinescu and Hine (2016); Doh and Choi (2017); Saberi-Pouya et al. (2017). BP is the most stable allotrope of phosphorus at room temperature and pressure. Few-layer phosphorene can be obtained through mechanical exfoliation method akin to graphene Brent et al. (2014); Hanlon et al. (2015). However, unlike graphene the phosphorene layers are not perfectly flat and form a puckered surface due to the sp3 hybridization of 3s and 3p atomic orbitals. Also, the band gap of BP is direct and can be tuned from 0.3 eV to the visible part of the spectrum Rodin et al. (2014a); Çakır et al. (2014). In contrast, graphene is gapless, and the transition metal dichalcogenides (TMDs) have an indirect gap in the bulk phase and only monolayer TMDs have a direct gapWang et al. (2012). Moreover, BP exhibits a strong in-plane anisotropy, which is absent in graphene and TMDs Carvalho et al. (2016).

In comparison to graphene, phosphorene is chemically reactive and tends to form strong bonds with surface of substrates which lead to some structural changes Brent et al. (2014); Shao et al. (2016); Vazifehshenas et al. (2010). Naturally, chemically stable 2DMs, such as graphene and BN, may be used for protecting fragile and low-chemical-stable 2DMs, such as phosphorene Bokdam et al. (2011).

The properties of substrate often drastically alter the transport behavior of the 2D crystal and the overall characteristics of the device. Recently, phosphorene has been transferred to the h-BN substrate Doganov et al. (2015); Moreno-Moreno et al. (2016). The interaction between phosphorene and substrate is considered to play a crucial role in the modulation of the electronic properties of phosphorene-based devicesGao et al. (2016); Doganov et al. (2015). In most currently available 2DMs, a sample lies on the top of a polar substrate such as h-BN, SiO2, SiC or Al2O3 Dean et al. (2010); Chen et al. (2015); Mogulkoc et al. (2016a); Ishigami et al. (2007); Dai et al. (2015); Zarenia et al. (2012); Buscema et al. (2014). In such heterostructures, the polar optical phonon modes of the substrate are localized near the 2DMs-substrate interface and are coupled to the surface optical (SO) phonon modes of the polar substrates through the long-range Fr$\ddot{\mathrm{o}}$hlich interaction. As a result, the SO phonon could be considered as the dominant plasmon-SO phonon coupling source in 2DMs on polar substrates.

The coupling of SO phonon modes to the plasma oscillations of free carriers is known as plasmon-SO phonon coupling. These coupled modes may be observed by the infrared (IR) transmission measurementsYan et al. (2013); Liu and Aydin (2016). The plasmon-SO phonon coupling changes the dips in the IR reflectivity spectra from isolated plasmon and SO phonon frequencies to the normal coupled ones Brar et al. (2014); Stroscio and Dutta (2001).

The coupled plasmon-SO phonon modes are extensively investigated in 2DMs with isotropic band structure such as graphene, silicene and germaneneChen et al. (2008); Perebeinos and Avouris (2010); Lin and Liu (2013); Wang et al. (2015). It is shown that this phenomenon modifies many-body properties such as plasmon modesHwang et al. (2010); Brar et al. (2014); Xiaoguang et al. (1985) and self energyHwang and Das Sarma (2013) and can also change the mobility Kaasbjerg et al. (2012); Ong and Fischetti (2012), effective massJalabert and Das Sarma (1989); Hwang and Das Sarma (2013), scattering rate Perebeinos et al. (2008); Lin and Liu (2013) and inelastic lifetime Ahn et al. (2014) in the carrier transport phenomena. Moreover, it can be considered as a mechanism for tuning the band gap Mogulkoc et al. (2016a); Fratini and Guinea (2008); Wang et al. (2015).

In this paper, we model a system of phosphorene on polar insulator substrates as shown in Fig. 1. We go beyond the assumption of independent electron and SO phonon modes and consider the coupled plasmon-SO phonon oscillations. Here, we theoretically study the coupled plasmon-SO phonon modes in monolayer and double-layer phosphorene on polar substrates using perturbation theory. We start from the expression for dynamical dielectric function of the coupled system in the random phase approximation (RPA) level of theory, and develop a general formalism, which includes the effect of anisotropic energy dispersion and rotationally misaligned double-layer system. In such systems having the interaction between the electrons and substrate SO phonons could yield to phonon-mediated electron-electron interaction which creates a new set of collective modes with highly anisotropic dispersion. We find that in the case of phosphorene in contrast to graphene Hwang et al. (2010); Brar et al. (2014) and TMDs Zhao et al. (2016), the mode coupling effect not only modifies the plasma dispersion relation but also enables us to tune the hybrid plasmon-SO phonon modes in two crystallographic directions. So, the anisotropy is an important feature of the coupled electron-SO phonon oscillations spectrum in phosphorene systems.

The paper is organized as follows. In Sec.II we present the theory to calculate the generalized dielectric function in the presence of electron-SO phonon interaction. From the generalized dielectric function, we calculate the coupled plasmon-SO phonon modes in Sec.III and present our results for the coupled plasmon-SO phonon modes in monolayer and double-layer phosphorene. Finally, we summarize the remarks of this work in Sec.IV.

II Theory

For an isolated phosphorene system, the anisotropic energy bands in the absence of electron-phonon interaction can be obtained from the $\mathbf{k}.\mathbf{p}$ free electron Hamiltonian:

[TABLE]

where Ec (Ev) is the energy of conduction (valence) band edge and the effective couplings between the bands are described by $\gamma$ and $\beta$ parameters. The $\eta_{c/v}$ and $\nu_{c/v}$ are related to the effective masses along $x$ and $y$ directions in the bands. These effective masses can be used to obtain an approximate anisotropic energy dispersion for phosphorene monolayer near the conduction (valence) band minimum (maximum)Rodin et al. (2014a).

In a phosphorene multilayer system with no electron-phonon interaction, the electrons in each layer interact with themselves and also with the electrons in other layers through the following electron-electron interaction Hamiltonian:

[TABLE]

where $v_{ij}(q)=v(q)e^{-qd_{ij}(1-\delta_{ij})}$ represents the diagonal (intralayer with $i=j$) and off-diagonal (interlayer with $i\neq j$) elements of the bare Coulomb potential matrix. Here, we define $v(q)=2\pi e^{2}/q\epsilon_{\infty}$ with $\epsilon_{\infty}$ being the high-frequency dielectric constant and $d_{ij}$ is the distance between i-th and j-th layers. Also, $a_{k,i}$ $(a_{k,i}^{\dagger})$ is the electron annihilation (creation) operator in layer $i$.

If the phophorene layers are supported by polar materials, an additional interaction term involving the electron-SO phonon coupling can be included by using Fr$\ddot{\mathrm{o}}$hlich Hamiltonian Tse and Das Sarma (2007):

[TABLE]

which renormalizes the screened electron-electron potential through the dielectric function matrix, $\epsilon_{ij}(\omega,\mathbf{q})$:

[TABLE]

In Eq. (3), $b_{\mathbf{q}\lambda}$ ($b_{\mathbf{q}\lambda}^{\dagger}$) is the phonon annihilation (creation) operator with wave vector $\mathbf{q}$ and branch index $\lambda$ and $M_{0}^{\lambda}(q)$ is the amplitude of the electron-SO phonon interaction Wang and Mahan (1972); Hwang et al. (2010):

[TABLE]

where $\omega_{SO}^{\lambda}$ is the SO phonon frequency of the $\lambda^{th}$ branch, $z$ is the vertical distance between phosphorene and substrate and we define:

[TABLE]

with $\epsilon_{o}$ being the zero-frequency dielectric constant. Here, we assume that the phonon-phonon interaction is negligible so each mode couples to the electrons, independently.

In order to study the effect of electron-SO phonon coupling on the collective charge-density excitations, one needs to obtain zeros of the determinant of the dynamical dielectric function matrix (the poles of the screened potential). We use the dynamical RPA dielectric function in which the contribution of the electron-SO phonon interaction is taken into accountZhang and Takahashi (1993):

[TABLE]

Here, $U_{ij}(\omega,q)$ is the combined Coulomb and phonon-mediated interactions between $i$th and $j$th layers:

[TABLE]

where $U_{0}(\omega,q)=v_{ph}(\omega,q)+v(q)$ and the SO phonon-mediated electron-electron interaction, $v_{ph}(\omega,q)$, is given by Sarma and Mason (1985):

[TABLE]

$D_{0}^{\lambda}(\omega)$ is the bare propagator for a phonon of branch index $\lambda$:

[TABLE]

Also, $\Pi_{j}(\omega,\mathbf{q})$ represents the noninteracting dynamic polarization function of layer $j$ which is given by the following expression:

[TABLE]

where $f^{j}(E_{\mathbf{q}})$ is the Fermi distribution function of jth layer at energy $E$ corresponding to the wave vector $\mathbf{q}$, $g_{s}=2$ is spin degeneracy and $\eta$ is the broadening parameter which accounts for the disorder in the system. We use the zero-temperature polarization function because this approximation is valid at the typical doping densities $n=5-50\times 10^{12}$ cm*-2* where the corresponding Fermi temperature ($T_{F}\approx 400-800$K) for phosphorene is higher than room temperature. However, the zero-temperature dynamic polarization function for intraband transitions in an anisotropic 2D material can be calculated by making use of the following anisotropic parabolic energy dispersion relation

[TABLE]

in Eq. (11). So one obtains

[TABLE]

We define $\mathbf{Q}=\sqrt{m_{d}/\hat{M}}(\mathbf{q}/k_{F})$ and $\mathbf{K}=\sqrt{m_{d}/\hat{M}}(\mathbf{k}/k_{F})$ where $\hat{M}$ is the mass tensor with diagonal elements $m_{x}$ and $m_{y}$ along $x$ and $y$ direction, $m_{d}=\sqrt{m_{x}m_{y}}$, $g_{2d}=m_{d}/\pi\hbar^{2}$ and $Z_{\pm}=((\hbar\omega+i\eta)/\hbar Qk_{F}\nu_{F})\pm(Q/2)$ with $\nu_{F}=\hbar k_{F}/m_{d}$. By introducing the rotational angle, $\tau_{i}$, as the angle between $x$-axis in the laboratory frame and $x$ direction of the $i$th layer, we can write $Q_{i}(\theta)=q\sqrt{m_{d}R_{i}(\theta)}/k_{F}$ in which the orientation parameter, $R_{i}(\theta)$, defined asSaberi-Pouya et al. (2016):

[TABLE]

In the case of monolayer, we will have $\tau_{i}=0$. The knowledge of the appropriate limiting behavior of the polarization function is important in investigating the excitation spectrum specially the collective excitations of the system. Therefore, we obtain the polarization function in the dynamic long-wavelength limit where the plasmon excitation are important due to their long lifetimes. In this limit, the polarization function of Eq. (13) can be approximated as (see the appendix A)

[TABLE]

In the following section, we present our calculations for the coupled plasmon-SO phonon modes in both monolayer and double-layer phosphorene systems with a number of experimentally chosen polar substrates/spacers.

III RESULTS AND DISCUSSION

III.1 Coupled plasmon SO-phonon modes in monolayer phosphorene

The RPA dielectric function for a monolayer system in which the electrons are coupled to SO phonons of a polar substrate can be obtained through summing over all the bare bubble diagrams as Hwang et al. (2010)

[TABLE]

In the long-wavelength limit $(q\to 0)$ by inserting Eq. (15) into Eq. (16), we get the following coupled collective modes (see appendix B):

[TABLE]

[TABLE]

where $\omega_{pl}(q,\theta)=\sqrt{2\pi ne^{2}R_{i}(\theta)q/\epsilon_{\infty}}$ is the plasmon frequency of uncoupled anisotropic system. We consider a phosphorene monolayer sandwiched by polar materials such as SiO2, h-BN or Al2O3 which their physical parameters are listed in Table 1. In addition, the effective masses in $x$ and $y$ directions of phosphorene monolayer are given as $m_{x}\approx 0.15m_{0}$ and $m_{y}\approx 0.7m_{0}$ where $m_{0}$ is the free electron massLow et al. (2014a).

In Fig. 2, we depict the coupled and uncoupled plasmon modes of phosphorene monolayer on SiO2 substrate calculated for two main crystallographic directions $\theta=0$ $(q\parallel x)$ and $\theta=\pi/2$ $(q\parallel y)$. According to this figure in the case of SiO2 as the substrate, there are three hybrid modes, one plasmon-like, $\omega_{-}(q,\theta)$, and two phonon-like, $\omega^{\lambda}_{+}(q,\theta)$, as were previously predicted theoretically Yan et al. (2013) and shown experimentally Yan et al. (2013) for graphene monolayer. It can be seen that the $\omega_{-}(q,\theta)$ (plasmon-like mode) is lower in energy than the uncoupled plasmon mode, $\omega_{pl}(q,\theta)$, and the lower (upper) phonon-like branch is starts from $\omega_{so}^{1}$ ($\omega_{so}^{2}$) at $q=0$ and increases below (above) the uncoupled plasmon dispersion by increasing $q$ along both directions. The coupling is strong when the frequencies of phonon-like modes deviate considerably from the bare SO phonon modes, $\omega_{so}^{\lambda}$ (see Eq. (17)). As expected in 2DMs with anisotropic band structure, the hybrid plasmon-SO phonon modes have higher resonance frequenciesRodin and Castro Neto (2015); Low et al. (2014b); Jin et al. (2015) and the plasmon-like mods are more affected by the electron-SO phonon coupling in the direction of the lighter mass, $q\parallel x$. In the coupled system, the electrons can be scattered either by the emission of hybrid plasmon-SO phonon modes ($\omega_{+},\omega_{-}$) or by the single particle excitation (SPE). The boundary of SPE continuum represented by green shaded area in Fig. 2 is given by:

[TABLE]

Though all hybrid plasmon-SO phonon excitations are damped at IR frequencies, the plasmon-like mode, $\omega_{-}$, experiences the Landau damping at smaller $q$ with respect to the phonon-like modes. On the other hand, the higher frequency phonon-like mode, $\omega^{2}_{+}$, stays away from $\omega^{2}_{SO}$ and enters the SPE region at large $q$ for both directions. Therefore, this mode shows a strong coupling and can be easily detected.

Moreover, at high densities (Fig. 2(b)) where the Fermi energy exceeds the SO phonon energy i.e. $E_{F}$ $\gg$ $\hbar\omega^{2}_{SO}$, the upper phonon-like branch strongly deviates from the uncoupled SO phonon’s energy while $\omega^{1}_{+}$ branch still remains between the two phonon energies. Thus, the high energy phonon-like mode, $\omega_{+}^{2}$, which is very sensitive to the density and direction, can be considered as a tunable quantity for applications in optical plasmonic devices.

Furthermore, we study the impact of angular orientation of $\mathbf{q}$ on the behavior of hybrid plasmon-SO phonon modes by plotting the loss function, $|\Im(1/\det\epsilon(\omega,q,\theta))|$, for $q$=0.1 nm (see Fig. 3). The results show that the maximum value of loss function occurs around $\theta$= 0 and 180 (along $x$ direction) for all three branches.

In order to explore the effects of specific substrate on the hybrid plasmon-SO phonon modes, we show the loss functions of phosphorene monolayer on (a) SiO2, (b) h-BN and (c) Al2O3 polar substrates along $x$ direction ($\theta=0$) in Fig. 4. One can see that using the Al2O3 (with higher $\alpha$) as a polar substrate results in the strong coupling as evidenced by considerable deviation of $\omega_{-}$ from $\omega_{pl}$ and of $\omega_{+}^{\lambda}$ from $\omega_{SO}$.

It sould be pointed out that while the phonon frequency $\omega_{SO}$ is an important parameter for the phonon-like modes, $\omega_{+}^{\lambda}$, the $\alpha$ parameter mostly affects the plasmon-like spectrum. Hence, the choice of substrate can be used to engineer the plasmon dispersion in phosphorene.

III.2 Coupled plasmon-SO phonon modes in double-layer phosphorene

Here, we consider a double-layer phosphorene with equal electron densities sandwiched by a homogeneous dielectric medium that models the substrate. It is reasonable to calculate the uncoupled plasmon modes before discussing the hybrid modes from the zeros of determinant of dielectric function matrix Eq. (7). In the leading-$q$ approximation (long-wavelength limit), two plasmonic branches are obtained through the following relationsRodin and Castro Neto (2015):

[TABLE]

[TABLE]

These two uncoupled branches are shown in Fig. 5(a) for $\theta=0$ (along $x$ direction) and in Fig. 5(b) for $\theta=\pi/2$ (along $y$ direction) when SiO2 is considered as a substrate/spacer Rudenko and Katsnelson (2014); Rodin and Castro Neto (2015).

It can be seen that the plasmon modes experience a stronger reduction along $y$ direction compared to $x$ direction. For the case of electron-SO phonon coupling in a double-layer phosphorene, we find two acoustic phonon-like modes, $\omega^{\lambda}_{ac(+)}$, with $\lambda$=1,2 and one acoustic plasmon-like mode, $\omega_{ac(-)}$ (see appendix B). The long-wavelength dispersions of the acoustic modes can be written as:

[TABLE]

[TABLE]

In the same vein, we obtain two optical phonon-like modes, $\omega^{\lambda}_{op(+)}$, and one optical plasmon-like mode, $\omega_{op(-)}$ as

[TABLE]

[TABLE]

The behavior of the coupled modes with SiO2 as substrate/spacer is presented in Fig. 5(c) and 5(d) for $\theta=0$ and $\theta=\pi/2$, respectively. The dispersion relations of $\omega_{op/ac(-)}$ modes resemble the uncoupled acoustic and optical plasmon modes with an additional $(1-\alpha e^{-2qz})^{1/2}$ multiplier which makes these coupled modes lie lower than the uncoupled ones.

Similar to the case of phosphorene monolayer, the hybrid plasmon-SO phonon frequencies along $q\parallel x$ are larger than along $q\parallel y$ because the carriers along $y$ direction have larger mass and so get damped faster. One can also see that the uncoupled plasmon modes for $q\parallel y$ lie lower than $\omega_{+}^{2}$ and therefore plasmon-SO phonon coupling along $y$ direction is substantially weaker compared to the $x$ direction.

In order to understand how the orientation parameter, $R_{i}(\theta)$, impacts the behavior of hybrid plasmon-SO phonon modes in the double-layer system, we show three cases in which $\tau_{1}=0$ and $\tau_{2}$ progressively increases from $\pi/4$ to $\pi/3$ and to $\pi/2$ along $x$ direction and for substrate/spacer SiO2 in Fig. (6). One may notice that, in general, with the reduction of $\tau_{2}$, the coupled plasmon-SO phonon modes become significantly stronger because the hybrid modes are larger along the lower mass direction, i.e. $\theta=0$. Furthermore, it is clear that the three acoustic plasmon/phonon-like modes are very sensitive to the rotation of layers and get damped as the angle of rotation is increased. Finally, we address the effect of separation between layers, $d_{12}$, on the coupled modes along $x$ direction in Fig. (7) for the substrate/spacer SiO2. Here, increasing the separation between layers shows a similar effect on hybrid mode frequencies as the reduction of rotation angle $\tau_{2}$ (see Fig. 6). As expected, the hybrid acoustic branches, similar to the uncoupled ones, depend strongly on the separation between layers at long-wavelengths and move to the optical branches by increasing the layers’ spacing. As a result, by adjusting two parameters $\tau_{2}$ and $d_{12}$, the acoustic branches may get strongly damped or a transition to the optical modes may be observed.

IV Conclusion

In summary, we have considered monolayer and double-layer phosphorene systems located on commonly used polar substrates (SiO2, h-BN and Al2O3) and calculated the anisotropic coupled plasmon-SO phonon dispersion relations. In the presented theory, the dynamical dielectric function is calculated within the RPA which includes the many-body electron-electron interaction in phosphorene layer(s) as well as the interaction between electrons and the long-range electric field generated by the substrate SO phonons. More importantly, we have obtained some analytical expressions for the plasmon-SO phonon dispersion relations in the long-wavelength limit which yield the same results of loss function from dielectric function. In the long-wavelength limit, three hybrid plasmon-SO phonon branches are obtained due to two relevant SO phonon modes of polar substrates in monolayer phosphorene. In case of double-layer, these hybrid modes are doubled with three acoustic, $\omega_{ac(\pm)}$, and three optical modes, $\omega_{op(\pm)}$. We have shown that these hybrid modes are considerably stronger along x direction because of lower effective mass of electrons. Moreover, by increasing the electron density, the hybrid excitation modes become stronger at both directions, simultaneously. Among mentioned substrates, we have found that Al2O3 as a polar substrate leads to more pronounced coupling of plasmon-phonon modes in comparison to h-BN and SiO2 for the electron density used here. Therefore, the phonon frequency of a substrate, $\omega_{SO}$, is the most effective parameter in phonon-like modes, $\omega^{\lambda}_{+}$, whereas $\alpha$ parameter changes the plasmon-like mode, $\omega_{-}$. Hence, the choice of substrate can be utilized in order to engineer the plasmon-SO phonon dispersion in phosphorene. We have also investigated the effect of misalignment of two layers on the hybrid modes in double-layer phosphorene and found that the acoustic phonon-like modes are more affected by rotation. In addition, we have observed that by decreasing separation between two layers the acoustic modes move away from the optical modes and quickly damped. As a result, the rotation angle and separation between two layers can be used as a mechanism for tuning the plasmon-SO phonon coupling effects.

Appendix A

The polarization function in long-wavelength limit

The dynamic polarization function given by Eq. (11) has the following expression in the long-wavelength limit at zero-temperature:

[TABLE]

One can rewrite Eq. (26) as:

[TABLE]

and find the approximate the relation:

[TABLE]

Finally, by ignoring unity in the denominator and making use of $Q_{i}(\theta)=q\sqrt{m_{d}R_{i}(\theta)}/k_{F}$ and $E_{F}=\hbar^{2}k_{F}^{2}/2m_{d}$, one obtains Eq. (15).

Appendix B

Electron-SO phonon modes in long-wavelength limit

B.1 Monolayer

Inserting Eq. (15) into Eq. (16), we have:

[TABLE]

From the zeros of dielectric function (here we drop the $q$ and $\theta$ for simplicity):

[TABLE]

the coupled plasmon-SO phonon modes are obtained as:

[TABLE]

By Taylor expanding the right hand side of Eq. (31), we get:

[TABLE]

[TABLE]

Since $\omega_{pl}\to 0$ at the long-wavelength limit, one can safely ignore $\omega_{pl}$ in the denominator of the above equations and obtain Eqs. (17) and (18).

B.2 Double-layer

Here, we calculate the coupled plasmon-SO phonon modes in double-layer systems with anisotropic band structure. In such systems, the rotation of one layer with respect to the other should be considered. In order to determine the dispersion relation of the coupled modes, we need to calculate the zeros of determinant of the total dielectric matrix (Eq. (7)) in the long-wavelength limit:

[TABLE]

where $\Pi^{\prime}(\omega,q,\theta)=\Pi_{i}(\omega,q,\theta)/R_{i}(\theta)$. By doing some algebra, we get the following relation:

[TABLE]

where $\Pi_{+}$ and $\Pi_{-}$ defined as

[TABLE]

and

[TABLE]

The dispersion relation of the couple modes is given by $\Pi^{\prime}(\omega,q,\theta)=\Pi_{\pm}(\omega,q,\theta)$. Using the Eqs. (5-10) in the $\Pi^{\prime}(\omega,q,\theta)=\Pi_{+}(\omega,q,\theta)$ condition, we get:

[TABLE]

where $f_{ac}(\theta)$ is defined as:

[TABLE]

After some algebra, we derive the following relation:

[TABLE]

which can be simplified in the long-wavelength limit as:

[TABLE]

and

[TABLE]

Substituting the relation for $f_{ac}(\theta)$ in above relations, one obtains Eqs. (22) and (23) for the coupled acoustic modes. In the case of the coupled optical modes, we use the $\Pi^{\prime}(q,\omega)=\Pi_{-}(q,\omega)$ condition and get the following equation:

[TABLE]

where $f_{op}(\theta)$ is given by:

[TABLE]

Solving Eq. (43) yields:

[TABLE]

In the long-wavelength limit, we end up with the following expressions for the coupled optical modes:

[TABLE]

and

[TABLE]

Finally, by inserting $f_{op}(\theta)$ into Eqs. (46) and (47) one obtains Eqs. (24) and (25).