Plasmons and screening in monolayer and multilayer black phosphorus

Black phosphorus exhibits a high degree of band anisotropy. However, we found that its in-plane static screening remains relatively isotropic for momenta relevant to elastic long-range scattering processes. On the other hand, the collective electronic excitations in the system exhibit a strong anisotropy. Band non-parabolicity leads to a plasmon frequency which scales as $n^{\beta}$, where $n$ is the carrier concentration, and $\beta<\tfrac{1}{2}$. Screening and charge distribution in the out-of-plane direction are also studied using a non-linear Thomas-Fermi model.

Introduction-Black phosphorus (BP) is one of the thermodynamically more stable phases of phosphorus, at ambient temperature and pressure. It is a layered material, with each layer forming a puckered surface due to sp 3 hybridization. In its bulk crystalline form [1][2][3][4][5], BP is a semiconductor with a direct band gap of about 0.3 eV with measured Hall mobilities in n and p−type samples approaching 10 5 cm 2 /Vs. Recent rediscovery of BP [6][7][8][9][10][11] in its multilayer form revealed highly anisotropic electrical and optical properties.
In this paper we examine the collective electronic excitations of BP, and its electrostatic screening behavior both along the in-and out-of-plane directions. We calculate the dielectric function (q, ω), at finite frequency ω and wavevector q, for monolayer and multilayers of BP within the Random Phase Approximation (RPA), using an effective low-energy Hamiltonian [12]. The inherent anisotropy of screnning is studied, and the out-of-plane screening properties of multi-layer BP flakes are considered within a non-linear Thomas-Fermi model. The 2D plasmon modes, which are obtained from the zeros of the dielectric function, or the electron loss spectra, [−1/ (q, ω)], show a highly anisotropic plasmon dispersion, ω pl (q, ω), and we studied its scaling behavior with doping. Lastly, we discuss the implications of our results for basic electrical and light scattering experiments.
Hamiltonian-BP has an orthorhombic crystal structure consisting of puckered layers. The lattice constant in the out-of-plane direction is about 10.7Å, and the effective layer-to-layer distance is half of this value [4]. In monolayer BP, translational symmetry in the z direction is broken, and its bandstructure has a direct energy gap at the Γ point instead of the Z point in the bulk case. Based on k·p theory and symmetry arguments, the inplane electron dispersion around the Γ point can be described by the following low-energy Hamiltonian [12], where η c,v and ν c,v are related to the effective masses, while γ and β describe the effective couplings between the conduction and valence bands. E c and E v are the energies of the conduction and valence band edges. At present, the energy gap for monolayer BP has not been measured experimentally, but recent ab initio calculation based on the GW method found an energy gap of ∼ 1.5− 2 eV [13,14]. Unlike other layered materials such as graphene and the transition metal dichalcogenides (TMDs), electrons in BP are energetically highly dispersive along the outof-plane direction. Cyclotron resonance experiments on bulk BP [15] found an out-of-plane effective mass considerably smaller than that of TMDs [16]. For multilayer BP, confinement in the out-of-plane z direction leads to multiple subbands. The in-plane dispersion within each subband j can be described by Eq. (1), where E c,v are being replaced with E j c,v . More explicitly, δE j c is given by j 2 2 π 2 /2m cz d 2 + δ c (d), where j labels the subband, d is the thickness of the BP film, and m cz is the electron effective mass along z. Analogous expressions apply also for the hole case. The quantities δ c,v (d) are chosen such that it reproduces the energy gap of the BP film [13], of 2 eV and 0.3 eV in the monolayer and bulk limit respectively. In this work, we adopt an average of experimental [15] and theoretically [15,17] predicted quantization mass i.e. m cz ≈ 0.2 m 0 and m vz ≈ 0.4 m 0 .
The in-plane dispersion is mainly determined by the parameters η c,v , ν c,v and γ. These parameters are chosen such that they yield the known anisotropic effective masses. In the bulk BP limit, we have m cx = m vx = 0.08 m 0 , m cy = 0.7 m 0 and m vy = 1.0 m 0 [4,15], and m cx = m vx ≈ 0.15 m 0 for monolayer BP [17]. Using this knowledge, we arrive at the following parameter set; η c,v = 2 /0.4m 0 , ν c = 2 /1.4m 0 , ν v = 2 /2.0m 0 , and γ = 4a/π eVm. The value of β is taken to be ≈ 2a 2 /π 2 eVm 2 [17], where a ≈ 2.23Å and π/a is the width of the BZ in x direction. tron gas in the RPA can be written as, where v c (q) = e 2 /2 0 q is the 2D Coulomb interaction and κ describes the effective dielectric constant of the medium, which for a common substrate, SiO 2 , is ∼2.5. Π(q, ω) is the 2D polarizability (i.e. the pair bubble diagram) given by, where k = k + q, {s, s } = ±1 denote conduction/valence bands, while {j, j } are the subband indices and g s = 2 is the spin degeneracy. E sjk and Φ sjk are the eigen-energies and eigen-functions after diagonaliz- Finite damping can be modeled with the phenomenological broadening term η. Allowed optical transitions between these quantized subbands occur when ss = ±1 (i.e. intra-and inter-band processes) and j = j . Otherwise the matrix element ... in Eq. (3) vanishes. Screening-In the static limit, we generalize the wellknown analytical form of the polarizability for 2D electron gas (2DEG) [18] to include anisotropy. Since ω = 0, we deal only with intraband processes. In the T → 0 and η → 0 limits, we have [19], where we have made the transformation, where M is the mass tensor with diagonal elements m x and m y , m d is the 2D density-of-states mass given by √ m x m y , and p F = √ 2m d µ/ . After some algebra, we arrive at where g 2D = m d /π 2 is the 2D density-of-states. We make an interesting remark: for q ≤ 2|k F ·q|, we see that Π(q) reduces to the familiar relation for the static polarization of a 2DEG, Π(q) = g 2D . Long range potentials, such as those induced by charged impurities, involve momenta q such that q ≤ 2|k F ·q|, so that screening will be isotropic, at least in the zero temperature and disorder limits.  6), with excellent agreement in the limits of the model. Π(q) has a kink at q = 2|k F ·q|. Fig. 1(b) illustrates how the kink migrates with change in doping. With increasing temperature and disorder, the kink is smoothed out as illustrated in Fig. 1(c)-(d), showing obvious deviation from the analytical model. The otherwise isotropic screening at small momenta now becomes anisotropic. On the other hand, dynamical screening, (q, ω), in BP exhibits strong directional dependence with q. Anisotropic dynamic screening might have important implications to carrier relaxation processes such as scattering with polar optical phonons. In Suppl. Info, we show the calculated real and imaginary part of (q, ω) at finite ω.
Plasmon dispersion-The zeros of the dynamical dielectric function (q, ω) yield the excitation spectrum of the plasmon modes of the electron gas. The loss function, defined as L(q, ω) = − [1/ (q, ω)], quantifies the spectral weight of the plasmon mode, which presents itself as a delta peak in the limit of zero damping. Experimentally, L(q, ω) can be quantified with EELS. In the long wavelength limit, i.e. q k F , these modes are welldescribed by classical Maxwell theory. We consider a BP film sandwiched between two dielectric media 1 and 2 . The bound modes, i.e. plasmons, are characterized by an in-plane wavevector q pointing at an angle θ with respect to x. The dispersion relation for the bound mode can be obtained from the solution to the following equation, In the limits θ = 0, π and σ xx = σ yy , Eq. 7 reduces tō In the non-retarded regime, i.e. q k 0 , hence k zi ≈ iq, we obtain the 'quasi-static' approximation, where κ = ( 1 + 2 )/2. For frequencies up to the midinfrared, the conductivity can be approximated by the Drude model, where D j is the Drude weight and i denotes the subbands.
Within the model Hamiltonian, the in-plane electron effective masses in vicinity to the Γ point are given by the following expressions [12], where ∆ i is the subband energy gap. Similar expressions apply for the hole case. Note that in graphene, D = µe 2 / 2 instead. With Eq. (10) and (11), we have the classical plasmon dispersion along the j = x, y directions, which is ω pl,j (q) = (D j /2π 0 κ)q. Fig. 2(a) plots the RPA loss function L(q, ω) for momentum along the two main crystallographic directions for monolayer BP, with an electron doping of 10 13 cm −2 . The plasmon disperses differently due to their mass anisotropy, where the smaller mass along x leads to higher resonance frequency. Classical plasmon dispersion agrees well with the RPA result in the long wavelength limit. Due to the energy gap of 2 eV for monolayer BP, Landau damping occurs preferentially via intraband processes. This occurs when plasmon enters the SP phase space, whose boundaries are given by, Our calculation suggests that the plasmon along the y direction is damped at mid-infrared frequencies, while the plasmon along x persists up to the near infrared.
The results reported here can be tested by EELS. In addition, plasmon modes in layered materials [20][21][22] can also be probed by Fourier transform infrared (FTIR) light scattering experiments of nanostructures [23,24] or with infrared nano-microscopy techniques [25,26]. For example, nanostructures exhibit prominent resonances in their extinction spectra due to localized plasmons with odd multiple of the momentum q = π/W where W can be the width of nanoribbons, or the diameter of nano-disks. Fig. 2(b) shows L(q, ω) for different angular orientation of q for a momentum corresponding to nanostructures of 100 nm in size. Dashed lines are solutions of Eq. 7. The results suggest polarization sensitive mid-infrared plasmonic resonances in the absorption spectra in BP nanostructures. Fig. 3 studies the scaling of plasmon frequency (along x) with carrier concentration n. For monolayer BP, we obtain the expected scaling relation of ω pl ∝ n 1/2 , as in conventional 2DEGs. However, for thicker samples, we found that ω pl ∝ n β , with β < 1 2 instead. This deviation is due to the strong non-parabolicity caused by interband coupling, particularly when the energy gap of the BP film is γ 2 /η c,v . Hence, non-parabolicity effects are more prominent for thicker films. We also note the general trend of increasing Drude weight with film's thickness due to the decreasing effective masses (see Suppl. Info.).
Screening and charge distribution in multilayers-We complete our study by considering the charge distribution and the electrostatic screening in few-layer BP sheets. For this aim we use a non-linear TF theory, which has been shown to properly account for the screening prop-erties of graphite [27][28][29] and MoS 2 [30]. We start by considering a given charge transfer between the substrate and the BP flakes, whose origin can be due to charge impurities in the substrate or to the action of a gate voltage. This charge transfer leads to a net surface charge density en 0 while a layer below the substrate acquires a charge of −en 0 , see inset of Fig. 4(b). For a BP sample of thickness d, the electrostatic potential V (z) and the carrier distribution n(z) as a function of the distance from the substrate z can be obtained from the energetic balance between kinetic and interlayer capacitance terms, which leads to the non-linear differential equation [30] where f (z) = [en(z)] 2/3 and we have defined and κ ≈ 8.3 are the interlayer separation and dielectric constant, respectively [4]. Using the boundary conditions f (0) = 5 2 β ⊥ en 0 and f (d) = 0, one can obtain the charge density from the solution of the integral equation On the other hand, the potential difference across a BP sample of thickness d can be shown to be given by [30]: where we have defined the dimensionless parameter r d = n 2/3 (d)/n 2/3 (0). The potential difference obtained from the above model is shown, for different carrier concentrations, in Fig. 4(a) for a n-doped sample (see Suppl. Info. for results also on p-doped samples). The screening of charged impurities or the gate potential increases as the thickness of the BP layer grows. The dependence of ∆V (d) on d suggests an intermediate screening behavior between the strong coupling limit of graphene, where the carriers concentrate close to the interface [28], and the weak coupling regime with reduced screening properties that dominates the screening of MoS 2 [30]. Our results suggest that the gate will have negligible effect 10 nm into the bulk of BP, consistent with recent experiments on multilayers BP transistors [9].
We have also calculated n(z) for a sample with a given thickness d but different charge carrier concentrations n 0 as shown in Fig. 4(b). We observe a strong dependence of the screening strength on n 0 , such that stronger screening is achieved for larger n 0 . Interestingly, from those results one could infer a screening length of the order of the inter-layer spacing for σ 0 = 10 13 cm −2 , whereas for lower concentrations, like 10 11 cm −2 , the screening length is one order of magnitude larger.
Conclusions-In conclusion, we have studied the screening properties of BP using a combination of RPA for the dynamic and static in-plane screening, as well as for the dispersion of the collective (plasmon) excitations, and a non-linear TF theory for the inter-layer screening. Whereas we find a relatively isotropic static screening, the band non-parabolicity leads to highly anisotropic plasmons. Most saliently, we find that for multilayer samples, the plasmon resonance scales with doping as n β , where β < 1 2 . Furthermore, the modes dispersing along one of the crystallographic directions are long lived, being Landau damped (i.e. decaying into intra-band electronhole pairs) only for high frequencies, near the infrared. Finally, we find that the charge distribution along the layers and the strength of the electric field screening in BP flakes seem to be between the strong coupling regime characteristic of graphene, and the weak coupling regime of the TMD semiconductors, such as MoS 2 .
Acknowledgements-FG and RR acknowledge support from the Spanish Ministry of Economy (MINECO) through Grant No. FIS2011-23713, the European Research Council Advanced Grant (contract 290846), and the European Commission under the Graphene Flagship, contract CNECT-ICT-604391. R.R. acknowledges financial support from the Juan de la Cierva Program.