Anisotropic plasmons, excitons and electron energy loss spectroscopy of phosphorene

TL;DR

This study investigates the anisotropic electron energy loss spectrum of monolayer phosphorene using ab-initio calculations, revealing anisotropic excitonic and plasmon peaks, and how doping affects intraband plasmon behavior.

Contribution

It provides a detailed theoretical analysis of anisotropic plasmons and excitons in phosphorene, including the effects of doping and effective mass ratios, using advanced ab-initio methods.

Findings

Anisotropic excitonic peaks near the bandgap in phosphorene.

Doping introduces intraband plasmon peaks within the bandgap.

Intraband plasmon dispersion follows a sqrt(q) behavior with anisotropic characteristics.

Abstract

In this article, we explore the anisotropic electron energy loss spectrum (EELS) in monolayer phosphorene based on ab-initio time dependent density functional theory calculations. Similar to black phosphorous, the EELS of undoped monolayer phosphorene is characterized by anisotropic excitonic peaks for energies in vicinity of the bandgap, and by interband plasmon peaks for higher energies. On doping, an additional intraband plasmon peak also appears for energies within the bandgap. Similar to other two dimensional systems, the intraband plasmon peak disperses as $\omega_{\rm pl} \propto \sqrt{q}$ in both the zigzag and armchair directions in the long wavelength limit, and deviates for larger wavevectors. The anisotropy of the long wavelength plasmon intraband dispersion is found to be inversely proportional to the square root of the ratio of the effective masses: $\omega_{\rm pl}(q \hat{y})/\omega_{\rm pl}(q \hat{x}) = \sqrt{m_x/m_y}$.

Figures

Equations

$${\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}^0(\mathbf{q},\omega )=\frac{1}{\Omega }\sum _{\mathbf{k}}^{\mathrm{BZ}}\sum _{n,n^{\prime }}\frac{f_{n\mathbf{k}}-f_{n^{\prime }\mathbf{k}+\mathbf{q}}}{\omega +{ϵ}_{n\mathbf{k}}-{ϵ}_{n^{\prime }\mathbf{k}+\mathbf{q}}+i\eta }\times$$

$$⟨{\psi }_{n\mathbf{k}}|e^{-i(\mathbf{q}+\mathbf{G})\cdot \mathbf{r}}|{\psi }_{n^{\prime }\mathbf{k}+\mathbf{q}}{⟩}_{{\Omega }_{\mathrm{cell}}}⟨{\psi }_{n\mathbf{k}}|e^{i(\mathbf{q}+{\mathbf{G}}^{\prime })\cdot {\mathbf{r}}^{\prime }}|{\psi }_{n^{\prime }\mathbf{k}+\mathbf{q}}{⟩}_{{\Omega }_{\mathrm{cell}}}.$$

$${\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q},\omega )$$

$${ϵ}_M^{-1}(\mathbf{q},\omega )=1-v_c(\mathbf{q}){\chi }_{\mathbf{G}=0{\mathbf{G}}^{\prime }=0}(\mathbf{q},\omega )$$

$${\mathcal{E}}_{\mathrm{Loss}}(\mathbf{q},\omega )=-\Im [{ϵ}_M^{-1}(q,\omega )]=v_c(\mathbf{q})\mathbf{\Im }[{\chi }_{\mathbf{G}=0,{\mathbf{G}}^{\prime }=0}(\mathbf{q},\omega )].$$

$$H=(u{k}_y^2+\Delta ){\sigma }_x+v_{\mathrm{f}}{k}_x{\sigma }_y,$$

$$E_s(\mathbf{k})=s\sqrt{v_{\mathrm{f}}^2{k}_x^2+(u{k}_y^2+\Delta {)}^2},$$

$$E_c(\mathbf{k})=E_c+\frac{{\hslash }^2{k}_x^2}{2m_x}+\frac{{\hslash }^2{k}_y^2}{2m_y},$$

$${\omega }_{\mathrm{pl}}(\mathbf{q})={\alpha }_0(E_{\mathrm{f}}-E_c{)}^{1/2}{\left[{\cos }^2{\theta }_q+\frac{m_x}{m_y}{\sin }^2{\theta }_q\right]}^{1/2}\sqrt{q},$$

$$\Pi (\vec{q},\omega )$$

$$\Pi (\mathbf{q}\to 0,\omega )$$

$$I_1$$

$$I_2$$

$$\Pi (\mathbf{q}\to 0,\omega )$$

$${\zeta }_1(E_{\mathrm{f}})$$

$${\zeta }_2(E_{\mathrm{f}})$$

$${\zeta }_3(E_{\mathrm{f}})$$

$$\nu$$

$$G_0(\nu )$$

$$\epsilon (q,\omega )=1-v_q\Pi (q,\omega )=0,$$

$${\omega }_{\mathrm{pl}}$$

$$\frac{{\omega }_{\mathrm{pl}}(q\hat{x})}{{\omega }_{\mathrm{pl}}(q\hat{y})}=\sqrt{\frac{v_{\mathrm{f}}^2\text{Re}{\zeta }_1(E_{\mathrm{f}})}{u[\text{Re}{\zeta }_2(E_{\mathrm{f}})+\text{Re}{\zeta }_3(E_{\mathrm{f}})]}}.$$

Full text

I Introduction

Within the family of 2D materialsGeim and Novoselov (2007); Schwierz (2010); Radisavljevic et al. (2011); et.al (2013); Dubertret et al. (2015); Liu et al. (2014a); Li et al. (2014a); Ghosh et al. (2015); Novoselov et al. (2016); Nahas et al. (2016); Mardanya et al. (2016), phosphorene (few layers of black phosphorous) - a direct bandgap semiconductor with a puckered atomic structure, has a highly anisotropic band-structure Li et al. (2014b); Liu et al. (2014b); Ling et al. (2015). It has a massive Dirac like energy dispersion along the armchair and a parabolic Schrödinger like energy dispersion along the zigzag directionQiao et al. (2014); Rodin et al. (2014a); Rudenko and Katsnelson (2014a); Ezawa (2014a). This results in highly anisotropic single particle electronicQiao et al. (2014), thermalLiu and Chang (2015), many particle excitonicTran et al. (2014) and plasmonic properties Low et al. (2014). In particular, the plasmon dispersion has been recently calculated based on low energy continuum Hamiltonian, and it is found to be highly anisotropicLow et al. (2014); Rodin and Castro Neto (2015) with different doping dependence in different directions, depending on the number of layers.

From an experimental viewpoint, momentum resolved electron energy loss spectroscopyEgerton (2009) (EELS) directly probes the loss function of a material, which in turn is simply the inverse of the imaginary part of the dynamical interacting dielectric constant of a material: ${\cal E}_{\rm Loss}({\bf q},\omega)=-\Im[1/\epsilon_{M}({\bf q},\omega)]$. It has been used extensively in a variety of materials such as grapheneEberlein et al. (2008); Shin et al. (2011); Kinyanjui, M. K. et al. (2012); Tegenkamp et al. (2011); Langer et al. (2010); Liou et al. (2015); et.al (2008); Yan et al. (2011), transition metal dichalcogenides et.al (2011); Johari and Shenoy (2011a), and bulk black phosphorousSchuster et al. (2015), to explore the single particle and collective excitations such as excitons and plasmons. In terms of computational methods, calculations based on effective low energy continuum Wunsch et al. (2006); Sachdeva et al. (2015) as well as the tight-binding Hamiltonian Hill et al. (2009); Lam and Guo (2015); Jin et al. (2015) are very insightful, but they fail to capture many of the experimental aspects of the EELS spectrum accurately; in particular the low energy intraband plasmons at large wave-vectors and the high energy interband plasmons involving transitions across various energy bands. For example, in doped graphene, the low energyWunsch et al. (2006); Hwang and Das Sarma (2007); Agarwal and Vignale (2015) and tight-binding approachHill et al. (2009) fails to capture the plasmon anisotropy at finite wave-vectors in the $\Gamma-K$ and the $\Gamma-M$ direction, which has been observed experimentally Tegenkamp et al. (2011); Langer et al. (2010). However, the loss function and corresponding plasmon dispersion relation is generally well described by density functional theory based ab-initio calculations Andersen and Thygesen (2013); Vacacela Gomez et al. (2016); Johari and Shenoy (2011b); Despoja et al. (2013); Kaltenborn and Schneider (2013); Andersen et al. (2014); Trevisanutto et al. (2008); Yan et al. (2011). In case of phosphorene, the low energy intraband plasmon spectrum has been studied using the continuum approximation Low et al. (2014); Rodin and Castro Neto (2015). Further, the effects of strainLam and Guo (2015) and disorderJin et al. (2015) on the plasmon spectrum have been incorporated using a tight-binding approach, but to the best of our knowledge, there is no ab-initio based study of the anisotropic plasmon dispersion (both low as well as high energy) and the high energy EELS spectrum of monolayer phosphorene. Motivated by this, we present an ab-initio density functional theory (DFT) based study of the anisotropic EELS spectrum, which includes low energy intraband plasmons at finite doping, excitons and interband plasmons in monolayer phosphorene.

We find that the crystal anisotropy of bulk black phosphorous is preserved down to it’s single layer, leading to a highly anisotropic electronic band structure, which results in a direction dependent EELS. In case of finite doping, we find an intraband plasmon mode which lies well below the band gap of phosphorene. Interestingly, it has a $\sqrt{q}$ dispersion for small wave vector in each of the two principal direction (parallel to the armchair and the zigzag edge of monolayer phosphorene), which signifies the two dimensional nature of the plasmon mode. We also find a highly dispersive mode in the EELS, which appears at slightly higher energy than the bandgap of phosphorene; it is identified as the exciton peak. Interestingly, while the exciton peak appears in the armchair direction, it is completely absent along the zigzag edge. The other high energy peaks correspond to different interband transitions, with a very distinct peak appearing for energies close to $5$ eV. We also observe a general trend that with increasing momentum transfer, all the resonant features (excitations) of the EELS spectrum are blue shifted and they gradually loose their strength.

The paper is organized as follows: in Sec. II we discuss the formulation for calculating the interacting density response function and the corresponding EELS spectrum, along with the computational details of our ab-initio study to get the electronic band structure of phosphorene. Next, we describe the calculated EELS spectrum in Sec. III, followed by a detailed discussion focused on anisotropy of intraband plasmons in Sec. IV. Finally, we summarize our findings in Sec. V.

II Theory and Computational details

II.1 Dynamical dielectric and loss function

Our starting point is the non-interacting density-density response function $(\chi^{0}_{\bf{GG^{\prime}}})$ for a periodic lattice, which is obtained from the Adler-Wiser formula given byAdler (1962); Wiser (1963),

[TABLE]

The Kohn-Sham energy eigenvalues $\epsilon_{n\bf{k}}$,the wave function $\psi_{n\bf{k}}$ and the corresponding Fermi-Dirac occupation function $f_{n\bf{k}}$ for the nth band at wave vector $\bf{k}$ are obtained from the ground state calculations performed using the density functional theory.

Within the framework of time dependent density functional theory (TDDFT) the interacting density-density response function can be obtained by solving a Dyson-type equation. Expanding in a plane wave basis (valid for a periodic system), the interacting response function can be expressed as

[TABLE]

where $\bf{G}$ and $\bf{q}$ are the reciprocal lattice vector and the wave vector, respectively and $K_{\bf{G_{1}G_{2}}}$ is the (2D truncated for monolayer phosphorene) Coulomb kernel. The exchange-correlation part of the kernel is neglected within the framework of the random phase approximation (RPA). Using the $\chi_{\bf{G}\bf{G^{\prime}}}$ matrix, the macroscopic dielectric function ($\epsilon_{M}$) is obtained as,

[TABLE]

where $v_{c}(\bf{q})$ is the (2D truncated for monolayer phosphorene) Coulomb Kernel as described in Ref. [Rozzi et al., 2006]. The dynamical loss function, which is directly related to the EELS, is calculated as

[TABLE]

Plasmons (collective density excitations) are the characterized by the zeroes of the real part of the macroscopic dielectric function (the denominator of the density-density response function within RPA).

II.2 Computational details

Electronic band structure calculations are performed using density functional theory, as implemented in the GPAW package Mortensen et al. (2005); et.al (2010); Bahn and Jacobsen (2002). The kinetic energy cut-off for the plane wave basis set is taken to be 500 eV. Initially, all the atomic positions are relaxed (using GGA) until the forces on each atom are less than 0.001 eV/Å. Next, we calculate the single particle states [to be used as input for evaluating the response function in Eq. (1)] using the GLLB-SC functional which includes the GLLB type exchange proposed by Gritsenko et al. and PBEsol correlation, which has been found to improve the bandgap in case of semiconductorsGritsenko et al. (1995). While a k-point grid of $16\times 22\times 6$ ( $16\times 22\times 1$ ) is used for the Brillouin zone integrations of the bulk (monolayer) structure, for the electronic band structure calculation, a much denser k-point grid of 32 $\times$ 44 $\times$ 14 ( 64 $\times$ 88 $\times$ 1) is used for calculating the $q$ dependence of the EELS, giving a momentum resolution of $\sim$0.043 Å*-1* and $\sim$0.0215 Å*-1* respectively for the bulk and the monolayer phosphorene. In case of the latter, we use a vacuum layer of 20 Å  in the direction perpendicular to the phosphorene plane to suppress any interaction between two replica images in the vertical direction.

EELS calculations are performed using the generalized random phase approximation (RPA). The generalized RPA uses the local field factors to add the impact of the exchange and correlation effects to the Hartree field Giuliani and Vignale (2005). Because of it’s long-range nature, the Coulomb potential of one layer can interact with its periodic replicas, which is avoided by taking a 2D truncated Coulomb kernel, following Ref [Rozzi et al., 2006]. We consider up to 50 empty bands to correctly describe all the electronic excitations. For the local field corrections, a cut-off energy of 50 eV is used for G and ${\rm\textbf{G}}^{\prime}$ vectors, which corresponds to 259 (135) plane waves for monolayer phosphorene (bulk black phsophorous). Doping or the change in carrier concentration is achieved by shifting the position of the Fermi energy ($E_{\rm f}$).

II.3 Electronic band structure of phosphorene

Phosphorene has a layered structure with each phosphorus atom covalently bonded with three adjacent atoms, forming a $sp^{3}$ hybridized puckered honeycomb structure [see Fig 1(a)]. As shown in the figure, mutually perpendicular armchair and zigzag directions are aligned along the $x$ and $y$ axis, respectively. The lattice parameters of the bulk unit cell are found to be $a=4.56$ Å  (armchair), $b=3.31$ Å  (zigzag) and $c=11.30$ Å,  which are in good agreement with the values reported in the literatureLiu et al. (2014b). As shown in Fig 1(a), $\Gamma-X$ and $\Gamma-Y$ are the high symmetry directions in the reciprocal lattice, which are aligned along the armchair and zigzag axes, respectively. Based on previous reports of large anisotropy of calculated and measured electronic and optical properties along these particular directions, in this paper we calculate and compare the EELS spectrum and the intraband low energy plasmon dispersion along $\Gamma-X$ and $\Gamma-Y$, respectively.

As reported in the literature, monolayer phosphorene is found to be a direct bandgap semiconductor, the magnitude of the gap being 1.51 eV (0.91 eV), calculated using the HSE06 (GGA-PBE) functional Qiao et al. (2014); Kumar et al. (2016); Ghosh et al. (2016). The GLLB-SC computed electronic bandstructure has a direct bandgap of 1.13 eV and it is shown in Fig 1, along with the corresponding density of states (DOS). The conduction band minima (CBM) and the valance band maxima (VBM) are located at the $\Gamma$-point. Interestingly, the electron dispersion is anisotropic and ‘semi-Dirac’ like around this point [see Fig 1(b)], having a massive Dirac character along the $\Gamma-X$ direction and parabolic with a large effective mass along the $\Gamma-Y$ direction De Sarkar et al. (2016). This highly anisotropic low energy dispersion of phosphrene is the origin of the direction dependent transport and optical properties of phosphorene.

Pristine phosphorene, being a relatively large bandgap semiconductor, has vanishingly small thermally excited charge carriers even at room temperature and consequently intraband plasmons are absent. In order to explore intraband plasmons we consider phosphorene doped via electrostatic doping (controlled by a varying gate voltage). For electrostatic doping, the bands near the CBM (and thus the DOS) remains unaffected but the Fermi energy shifts to the conduction band giving rise to doped (electronic) charge carriers, whose number density can be tuned by controlling the gate voltage.

The effective low energy Hamiltonian and band-structure of monolayer phosphorene has been derived using ${\bf k}\cdot{\bf p}$ methodRodin et al. (2014b) as well as the tight-binding approach Rudenko and Katsnelson (2014b); Rudenko et al. (2015); Ezawa (2014b); Pereira and Katsnelson (2015), both of which yield a qualitatively similar picture De Sarkar et al. (2016). For this manuscript we work with the bare minimum effective low energy Hamiltonian of phosphorene, retaining only the lowest order terms in the wavevectors Ezawa (2015),

[TABLE]

where $\sigma_{i}$’s are the Pauli matrices. The anisotropic energy spectrum is thus given by

[TABLE]

where $s=+1$ $(s=-1)$ corresponds to the conduction (valance) band. Fitting Eq. (6) to our GLLB-SC dispersion in vicinity of the $\Gamma$ point [see Fig. 1(c)] yields $\Delta=0.56$ eV, $u=3.55$ eV/Å2  and $v_{\rm f}=4.75$ eV/Å. Note that, Eq. (6) can further be approximated as an anisotropic parabolic dispersion given by

[TABLE]

where $E_{c}=\Delta$. The values of the anisotropic effective masses for the conduction band are given by $m_{x}\equiv\hbar^{2}\Delta/v_{\rm f}^{2}=0.20m_{e}$ and $m_{y}\equiv\hbar^{2}/(2u)=1.1m_{e}$, with $m_{e}$ being the electrons rest mass, consistent with earlier studies Qiao et al. (2014). Equations (5)-(6) will be used to obtain the low energy and low momentum transfer plasmon dispersion of monolayer phosphorene analytically which in turn will be compared with our ab-initio results.

III Electron energy loss spectrum

Having discussed the ab-initio and the low energy electronic band structure, we now proceed to calculate the EELS spectrum of phosphorene. We first test the methodology by calculating the EELS spectrum of bulk black phosphorous and comparing the same with reported data Schuster et al. (2015); Wu et al. (2015). The computed EELS spectrum of bulk black phosphorous is shown in Fig. 2. As expected, the anisotropic nature of the electronic bandstructure of bulk phosphorous is manifested in the EELS data. The high energy interband plasmon peak, present in the vicinity of 20 eV in both the $\Gamma-X$ and the $\Gamma-Y$ direction, is in good agreement with the recently published experimental EELS results for bulk black phosphorous Wu et al. (2015). The low energy part of the spectrum, representing intraband plasmons and excitons, are shown in the inset of Fig. 2, and a marked anisotropy along the two principal direction is observed. In particular, the low energy exciton peak in black phosphorous ($E\approx 0.7$ eV) is only present in the $\Gamma-X$ (armchair) direction, which is consistent with recently reported measurements Schuster et al. (2015). Also, the intraband plasmons in the armchair direction ($\Gamma-X$) have a higher intensity as opposed to those in the zigzag ($\Gamma-Y$) direction.

After benchmarking the methodology by successfully reproducing the characteristic features of EELS of bulk black phosphorus, next we focus on monolayer phosphorene. Due to obvious reasons related to it’s electronic band structure [see Fig 1 (b)], EELS of monolayer phosphorene is also found to be highly anisotropic. This can clearly be observed in Fig. 3(a) and (b), where the calculated spectrum of a single layer of doped (by taking $E_{\rm f}=50$ meV, measured from the CBM) phosphorene is plotted for momentum transfer ($q$) along the $\Gamma-X$ (armchair) and $\Gamma-Y$ (zigzag) directions, respectively. Note that, while the lowest energy peak due to the intraband plasmons are observed only in case of finite doping, others appearing at higher energy have intrinsic origin related to excitons and interband plasmons and they are present in undoped phosphorene as well. Specific features of the EELS are discussed in detail in the following subsections.

III.1 Intraband plasmons

As shown in Fig. 3 and Fig. 4, the first peak appears for energies less than 0.25 eV. The energy corresponding to the peak is well below the bandgap of the pristine monolayer phosphorene and in addition this peak is absent in case of undoped phosphorene. Thus we interpret this as the peak originating from the intraband plasmon modes. With increasing momentum transfer, the intensity of intraband plasmons decreases and the peak position shifts to a higher energy, in both $\Gamma-X$ and $\Gamma-Y$ direction [see Fig. 3 and Fig. 4]. Similar type of blue shift is observed as the doping is increased. Alike bulk black phsophorous, the intensity of the intraband peak is lower for momentum transfer along the $\Gamma-Y$ direction, than compared to the $\Gamma-X$ direction, which has the maximum intensity among all the EELS peaks at $q\rightarrow 0$ [see Fig. 3]. The anisotropy of intraband plasmon modes and their momentum and doping dependence are discussed in detail in Sec. IV.

III.2 Low energy excitons

The next peak in the EELS is observed in the vicinity of $1.3$ eV, which is gradually shifted towards the higher energy, accompanied by intensity reduction with increasing $q$. Note that, this peak is much more dispersive compared to the intraband plasmon peaks and extends approximately upto 3 eV (see Fig. 3 and Fig. 4). Since the peak energy coincides approximately with the energy gap of monolayer phosphorene and it exists in undoped phosphorene as well and it is almost independent of doping, we identify this as the lowest energy excitonic peak. As shown in Fig. 3 and Fig. 4, the excitonic peak is highly anisotropic in nature, as it appears only in case of momentum transfer along the $\Gamma-X$ direction, while being completely absent along the $\Gamma-Y$. We believe that the strongly anisotropic optical response reported for monolayer black phsophorene, which is transparent to the incident light in the energy range of 1.1-2.5 eV, but only if it is polarized along the zigzag direction and opaque if the light is polarized in the armchair directionTran et al. (2014); Schuster et al. (2015), originates from the anisotropy of the lowest energy electron-hole pair excitation observed in this work. Note that, since the crystal anisotropy observed in monolayer persists for multilayer, as well as bulk black phosphorous, similar kind of anisotropic excitonic response is expected at higher thickness as well.

III.3 Interband high energy plasmons

As we move towards higher energy, the next prominent EELS peak appears in the vicinity of $5$ eV, and similar to the low energy peaks, it’s intensity reduces accompanied by a blue shift with increasing momentum transfer. Though it is also anisotropic in nature, this peak has a higher intensity along the zigzag direction as compared to the armchair direction – unlike the case of intraband plasomons. As expected, this high energy peak is nearly independent of doping, as it corresponds to very high energy interband transitions.

Comparing the EELS intensity of various peaks as a function of $q$, it is clear that in general the EELS intensity is maximum for direct transitions with $q\to 0$ and decays with increasing $q$. However, the low energy intraband plasmon peak decays more rapidly with increasing $q$ (on account of damping by electron-hole excitations) as compared to the high energy peaks associated with interband transitions. For example, in case of $q\parallel\Gamma-X$, the intraband plasmon peak is the most intense one among all the EELS peaks for $(q=0.021~{}$Å*-1*), which almost vanishes at higher $q$, leaving the interband transition peak around 10-14 eV to be the most prominent one. Owing to it’s anisotropic nature, the scenario is different in case of $q\parallel\Gamma-Y$, where the interband transition peak around 5 eV has the highest intensity (at $q=0.021~{}$Å*-1*) among all the EELS peaks and it broadens significantly and shifts to around 10-14 eV at higher value of $q$.

Note that, other than the intraband plasmons, rest of the peaks corresponding to the interband transitions are likely to be affected due to the bandgap underestimation (approximately 40%) of GLLB-SC based electronic band structure calculations. For example, while the GW bandgap is reported to be 1.84 eVGhosh et al. (2016); Qiao et al. (2014) for monolayer black phosphorus, it is found to be 1.13 eV in our calculation. Thus, the EELS peaks corresponding to interband transitions are expected to be blue shifted in an actual experiment. However, the GLLB-SC bandgap estimation of $0.35$ eV for bulk black phosphorous is very close to its experimentally reported value of $0.31-0.35$ eV Ghosh et al. (2016). Thus it turns out that for bulk black phosphorous, EELS peak of 20 eV (as shown in Fig. 2) based on GLLB-SC calculations is consistent with the 20 eV peak observed in recent experimentsWu et al. (2015).

IV Anisotropic intraband plasmons

Having described the prominent features of the EELS spectrum of monolayer phosphorene, we now focus on the anisotropic nature of the low energy intraband plasmons and investigate their $q$ and doping ($E_{\rm f}$) dependence. The low energy loss function along the $\Gamma-X$ and $\Gamma-Y$ directions, is shown in Fig. 4 (a) and (b), respectively. As discussed earlier, the first low energy peak which appears for energies significantly below the bandgap of phosphorene and is only present in case of finite doping, corresponds to intraband plasmons - or collective charge density excitations. As shown in Fig 4, the intra-band plasmon peak has higher intensity along the $\Gamma-X$ direction as compared to the $\Gamma-Y$ direction, similar to the case of bulk phosphorus. Interestingly, the intensity decay and blue shift of the intraband plasmon peak is more rapid in the $\Gamma-X$ direction, than compared to the $\Gamma-Y$ direction.

The momentum dependence of the intraband plasmon peak is further analyzed at different doping (by varying the Fermi energy) in Fig. 5. As shown in the figure, for relatively small momentum transfers, the dispersion follows the universal long wavelength $\sqrt{q}$ behavior which is ubiquitous in two dimensional systems Sachdeva et al. (2015). For higher value of the momentum transfer $q$, a clear deviation from the $\sqrt{q}$ behavior is observed. Eventually the plasmon dissipates into single particle continuum of the conduction band, whose boundary is marked by the maxima of $E^{sp}_{c}({\bf q})=E(\Gamma+{\bf q})-E(\Gamma)$ [see the shaded area in Fig. 5]. As shown in Fig. 5, increasing the doping extends the range of validity of $\sqrt{q}$ behavior. For example, at $E_{\rm f}=50$ meV the $\sqrt{q}$ fit along the $\Gamma-X$ direction holds good upto $q=0.1$Å*-1*, which extends upto $q=0.15$Å*-1* with increased doping ($E_{\rm f}=90$ meV). Similar behaviour of the plasmon dispersion is also seen along the $\Gamma-Y$ direction, albeit the spectral weight of the corresponding plasmon peak is smaller in the $\Gamma-Y$ direction.

Starting from the anisotropic parabolic approximation of the bandstructure of monolayer phosphorene, given by Eq. (7), the low energy plasmon dispersion within the random phase approximation, (for the Coulomb potential $V_{q}=2\pi e^{2}/\epsilon q$) has been obtained in Ref. [Rodin and Castro Neto, 2015]. It is explicitly given by

[TABLE]

where we have defined $\theta_{q}=\tan^{-1}{(q_{y}/q_{x})}$, $\alpha_{0}^{2}=2\pi e^{2}g_{\rm 2d}/(m_{x}\epsilon)$, and $g_{\rm 2d}=\sqrt{m_{x}m_{y}}/(\pi\hbar^{2})$ is the two dimensional density of states for an anisotropic parabolic band system. Eq. (8) directly yields the following: 1) $\sqrt{q}$ dependence of the plasmon dispersion in all directions for small wavevectors, 2) $\sqrt{E_{\rm f}-E_{c}}$ dependence of the low energy (and low wave-vector) plasmon dispersion on the Fermi energy, and 3) the long wave-length anisotropy of the plasmon dispersion $\omega_{\rm pl}(q{\hat{x}})/\omega_{\rm pl}(q{\hat{y}})=\sqrt{m_{y}/m_{x}}$, independent of the doping in the system.

The doping dependence of the intraband plasmon dispersion is shown in Fig. 6, and it clearly shows a reasonable match of the plasmon dispersion with the expected $\omega_{\rm pl}\propto\sqrt{E_{\rm f}-E_{c}}$ dependence of the Fermi energy for small $q$. The plasmon frequency anisotropy ratio is shown in Fig. 7 and it also seems to be more or less consistent with the $\omega_{\rm pl}(q{\hat{x}})/\omega_{\rm pl}(q{\hat{y}})=\sqrt{m_{y}/m_{x}}$ behavior. A more thorough calculation of the long wavelength plasmon dispersion using the semi-Dirac continuum Hamiltonian for phosphorene (see Eq. (5)), is presented in the appendix, and it also yields qualitatively similar results. However, the anisotropy ratio now explicitly depends on the Fermi energy – see Eq. (18). However we note that this $E_{\rm f}$ dependence of the anisotropic ratio of the plasmon frequency is small – as shown in panel b) of Fig. 7.

V Summary and conclusions

In this article we present a thorough study of the anisotropic EELS spectrum of monolayer phosphorene using the TDDFT framework. We find that the anisotropy of the underlaying phosphorene crystal leads to the anisotropy in the band structure and consequently in the EELS spectrum as well – similar to the case of bulk black phoophorous. For finite doping in the system, the lowest energy peak in the EELS corresponds to the intraband plasmon mode (charge density excitations) in the sub electron volt range. At slightly higher energy than the band gap, there is a highly dispersive low energy exciton mode, which is almost independent of the doping and it is absent for the momentum transfer along the $\Gamma-Y$ direction. At even higher energies and completely independent of the doping, there are various plasmon modes arising from the interband transitions with a very distinct peak appearing at $\approx 5$eV for monolayer phosphorene (as per GLLB calculations). With increasing number of phosphorene layers (or thickness) this mode is likely to shift to higher energies eventually merging into a 19 eV peak observed in the bulk black phosphorous. We explore the low energy anisotropic intraband plasmons in detail and compare their behavior to analytical expression of the corresponding plasmon dispersion arising from the effective low energy dispersion using RPA. The anisotropic intraband plasmon modes are found to be highly dispersive in nature, with the large wavelength limit following the $\omega_{\rm pl}\propto\sqrt{q}$ behavior in all directions which is a universal characteristic of plasmons in two dimensions. With increasing doping, the long wavelength plasmon dispersion is found to scale with the Fermi energy as $\omega_{\rm pl}\propto\sqrt{E_{\rm f}-E_{c}}$. The anisotropy of the long wavelength plasmon dispersion is found to be proportional to the ratio of square root of the effective masses: $\omega_{\rm pl}(q\hat{y})/\omega_{\rm pl}(q\hat{x})=\sqrt{m_{x}/m_{y}}$.

Acknowledgements

We acknowledge funding from the following DST (Department of science and technology, government of India) schemes: 1) DST INSPIRE Faculty Award, 2)SERB Fast Track Scheme for Young Scientist 3) Ramanujan Fellowship, and 3) DST Nanomission project. We also thank CC IITK for providing HPC facility.

Appendix A Low energy plasmon dispersion of phosphorene

Here we calculate the low energy plasmon dispersion for a massive semi Dirac system in two dimensions, described by Eq. (5). The polarization (density-density response) function is given by

[TABLE]

where $\omega^{+}=\omega+i\eta$, $E_{k,s}=s\sqrt{v_{\rm f}^{2}k_{x}^{2}+(uk_{y}^{2}+\Delta)^{2}}$, ${\bf k}^{\prime}={\bf k}+{\bf q}$, $f^{ss^{\prime}}({\bf k},{\bf q})$ is the overlap function of spinors and $n_{\rm F}(E_{k,s})$ is the Fermi function. Expanding the polarization function defined in Eq. (A), upto $q_{x}^{2}$ and $q_{y}^{2}$ at zero temperature, leads to

[TABLE]

where $\alpha=\sqrt{(E_{\rm f}-\Delta)/u}$, $\beta=v_{\rm f}^{-1}\sqrt{E_{\rm f}^{2}-(\Delta+uk_{y}^{2})^{2}}$ and $I_{1,2}(k_{x},k_{y})$ are complex functions defined as follows

[TABLE]

Performing the integral in Eq. (A) we obtain,

[TABLE]

where we have defined

[TABLE]

Equation (A) in turn uses the following notation,

[TABLE]

where $\text{E}\left[\phi,-\frac{1}{\nu^{2}}\right]$ and $\text{F}\left[\phi,-\frac{1}{\nu^{2}}\right]$ are incomplete elliptic integral of the first and second kind respectively.

Within RPA, the plasmons modes are given by the zeros of the longitudinal dielectric functions,

[TABLE]

where $v_{q}$ is the Fourier transform of the Coulomb potential. Substituting Eq. (13) in Eq. (16) immediately yields the long wavelength plasmon dispersion to be

[TABLE]

where we have used the two dimensional form of $V_{q}=2\pi e^{2}/(\epsilon q)$ and defined $\beta_{0}^{2}=ge^{2}/(2\pi\hbar^{2}\epsilon)$. Figure 8, shows the angular dependence of the long wavelength plasmon frequency, specified by Eq. (17), and it matches reasonably well with the long wavelength expression given in Eq. (8).

The ratio of the anisotropic plasmon dispersion in $x$ and $y$ direction for same value of the wave vector (in the long wavelength regime) is given by

[TABLE]