Unusual plasmonic responses in phosphorene with topological transition: the interplay of strain and doping

In this paper, plasmonic responses of phosphorene in the presence of strain and doping have been systematically investigated. Based on density functional theory, permittivities include both the intraband and interband transitions of electrons have been calculated. Due to the modification of the band structure, significantly higher Drude plasma frequency has been observed along the zigzag direction, other than the armchair direction as in the usual case. The resulting unusual plasmonic responses change their anisotropy, both in the elliptic as well as the hyperbolic regimes. Based on our calculations, positive strain as large as 5% along the zigzag direction can even lead to so-called reversed hyperbolic plasmonic responses. The k-surfaces of the plasmonic modes in extended monolayer have been analytically solved, and it is found that actively switching the topology (between elliptic and hyperbolic regimes) of the plasmonic responses by changing the Fermi level is possible in phosphorene at certain frequencies. In the end, a simple model has been proposed to describe such plasmonic responses in the infrared and the parameters of the model have been listed in tables which can be used directly in calculating the permittivities. Our studies may extend the scope of existing investigations of phosphorene plasmons and lead to band engineering as a way to control plasmons in two-dimensional materials.


Introduction
Several years ago, monolayer black phosphorus (bP) or phosphorene has been rediscovered as an anisotropic narrow-gap two dimensional (2D) material with relatively high mobility [1]. Due to its puckered structure, phosphorene shows strongly polarization-dependent optical properties within the monolayer [1,2]. Electrons in this material possess largely different effective masses along the zigzag and the armchair directions, leading to anisotropic collective excitations known as plasmons [3][4][5][6]. The plasmonic responses of phosphorene can be deduced from the Drude model, in which different Drude weights are applied for the conductivities along two directions [3,7]. Based on such model, plasmonic devices, such as waveguides, absorbers, etc, have been proposed by researchers [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23], showing promising applications in optics. Usually, higher Drude plasma frequency occurs along the armchair direction. Further investigations on bP plasmons show that such material even supports hyperbolic plasmonic responses, where the imaginary parts of the conductivities or the real parts of the permittivities along zigzag and armchair directions have opposite signs, i.e. one is negative and the other is positive [24][25][26][27]. Hyperbolicity in bP may lead to exotic plasmonic phenomena. Usually, in the hyperbolic regime, the permittivity along the armchair/zigzag direction is negative/positive [26]; however, this is not always true if strain and doping are both applied. There could be the situations where higher Drude plasma frequency occurs along the zigzag direction and the permittivity along the zigzag/armchair direction is negative/positive in the hyperbolic regime, i.e. so-called reversed hyperbolicity [28]. The occurrence of such significant changes in plasmonic responses in the elliptic as well as the hyperbolic regimes is due to the coexistence of strain and doping, which changes the shape of the bands and eventually the Drude plasma frequencies along zigzag and armchair directions. The combination of strain and doping often leads to exotic phenomena in 2D materials [29][30][31]. For example, Yuan et al have theoretically discovered the reversible phase transition induced by charge doping and strain in transition metal dichalcogenides [32]; such phenomena induced by strong orbital-lattice-charge coupling can be found only if both strain and doping are applied. In this paper, elliptic-to-hyperbolic topological transition of plasmonic responses with negative permittivity along the zigzag direction in the hyperbolic regime has also been discovered at specific frequency, which can only occur with proper strain and doping.
Materials with hyperbolic plasmonic responses are rare in nature; bP perhaps is the first of this kind of 2D materials investigated by researchers [24][25][26][27]. The occurrence of hyperbolic plasmons in electron doped phosphorene is mainly due to largely different Drude plasma frequencies caused by intraband transition of electrons along the zigzag and the armchair directions [3][4][5][6][7]. Although the interband transition of electrons is relatively weak in phosphorene especially when the polarization of light is along the zigzag direction, the induced permittivity cannot be neglected and is essential to hyperbolic responses. Simply replacing the interband part of the permittivity by a constant is a good approximation for deriving the total permittivity at low frequencies, which shows the existence of a narrow frequency range where hyperbolic plasmons are supported. Apart from the fact that phosphorene is unstable under ambient condition, hyperbolic plasmons have hardly been observed in experiment, which is possibly due to the lack of systematic studies on how such hyperbolicity responds to strain and doping beyond the effective mass approximation. Investigations on the mechanism behind the hyperbolic plasmons and how such plasmonic responses change with strain and doping might be helpful in finding them during the experiment.
In this paper, we have systematically investigated the plasmonic responses in phosphorene under different levels of strain and doping. The permittivities include both the interband and intraband transitions of electrons have been calculated based on density functional theory (DFT). The strains are applied either along the zigzag or the armchair directions up to 10%, modifying the band structures, and the Fermi energy is shifted upwards by as large as 2 eV representing the effect of electron doping. The band structures are stretched. Since the Drude plasma frequencies (denoted as ω D ) only involve states near the Fermi energy, they are also significantly modified. We have extracted the Drude plasma frequencies and compared their values along zigzag and armchair directions. Normally, higher ω D along the armchair direction is observed, which corresponds to the usual plasmonic responses [7]; however, with both strain and doping present in the system, situations with significantly higher ω D along the zigzag direction become possible, leading to what we refer to as unusual plasmonic responses in this paper. The interplay of strain and doping creates some kind of boundary beyond which one can expect the occurrence of such unusual plasmonic responses. The Drude plasma frequency only determines the intraband part of the permittivity and can only explain the change of anisotropy in the elliptic regime, while the hyperbolic plasmonic responses depend on the total permittivity. Observations on how the Drude plasma frequencies change as functions of the strain and doping may help us narrow down the range if one wants to search for, for example, hyperbolic responses under particular strain and doping conditions. What is more, the intraband transition induced part of the permittivity is always negative, which alone cannot explain the hyperbolicity; the interband transition induced part of the permittivity is essential and can be regarded as positive background. Calculations indicate that 5% strain along the zigzag direction can lead to remarkable unusual plasmonic responses even in the hyperbolic regime. For clear demonstration, the in-plane permittivities have been used to derive the 2D conductivities using a simple model, based on which the dispersion relations of the plasmonic modes in extended monolayer have been solved, and the plotted isofrequency contours are used to directly illustrate the topology of the plasmonic responses. At specific frequencies, actively switching between the elliptic and hyperbolic regimes via doping is possible. In the end, a simple model which treats the interband transition part of the permittivities as background constants has been proposed, aiming to derive the total permittivities at low frequencies without deploying the time-consuming DFT calculations. The parameters of the model have been listed in tables in the supplementary material (https://stacks.iop.org/NJP/23/113036/mmedia) and can be directly used in applications. Our study may extend the scope of existing investigations of bP plasmons and lead to band engineering via strain and doping as a way to control plasmons in 2D materials. where the x and y axes are along the zigzag and armchair directions, respectively. Strains are applied along two directions, from 0% to 10%, as indicated by the red and blue double arrows. (b) Bandgaps as well as energy differences between the conduction and valence bands at Γ under different strains along x and y axes. Black and red curves with scatters respectively denote the bandgap under various strains along the x and y axes, while light blue and yellow bars respectively denote the energy differences at Γ. The highest valence bands and the lowest conduction bands under strains from 0% to 10% along the x and y axes are plotted in (c) and (d), respectively. Figure 1(a) illustrates the crystal structure of phosphorene, where a black rectangle denotes the corresponding primitive cell. The zigzag and armchair edges are along the x and y axes, respectively. For simplicity, positive strains up to 10% have been applied either along the zigzag or the armchair direction, as respectively illustrated by the red and blue double arrows in figure 1(a). Without any strain, phosphorene is a semiconductor with a direct gap located at Γ in the Brillouin zone. With small strains, bP is still a direct-bandgap material; direct-to-indirect bandgap transition will happen with strains 10% along the x direction and 7%-10% along the y direction according to our calculations. For completeness, we plot both the actual bandgaps and the energy differences between the conduction band minimum (CBM) and valence band maximum (VBM) at Γ as functions of strains in figure 1(b). Light blue and yellow bars denote the energy differences, and the black and red curves with scatters denote the corresponding actual bandgaps. Before reaching the threshold of direct-to-indirect bandgap transition, these two should have the same value. Under 10% strain along the x direction and 7%-10% strains along the y direction, bP becomes an indirect bandgap material; the bandgap is smaller than the energy difference between CBM and VBM at Γ. One might notice that the bandgap has a peak near 4% strain along the x direction, which is below the threshold of direct-to-indirect bandgap transition and bP is still a direct-bandgap material. This is caused by the inversion of the first and the second conduction bands. Near the Γ point, the first conduction band is dominated by the bonding states of the 3P z orbits in P atoms; while for the second conduction band, it is dominated by anti-bonding states of the 3P y orbits. With positive stains, the distances between the atoms will basically increase, which causes the energy of the first/second conduction bands becomes higher/lower. Eventually near 4% strain along zigzag/x direction, band inversion will happen: the first/second conduction band becomes the second/first one; the curvature near the bottom is thus totally changed. Band inversion explains the peak near 4% before the direct-to-indirect transition. Further increasing the strains beyond 10% could lead to instability of the material. Strains can significantly modify band structures, which should also largely affect the optical properties as well as plasmonic responses after being doped with electrons. Figures 1(c) and (d) show the highest valence bands and the lowest conduction bands under strains along x and y axes from 0% to 10%, respectively. Not only the gap widths, but also the geometry of the bands changes. The curvature near the bottom of the first conduction band at Γ is related to the electron mass, which is clearly modified especially under strains along the zigzag/x direction. With high level of doping, the Fermi energy is no longer near the bottom/top of the conduction/valence bands, thus the effective mass approximation of the bands is no longer valid; one has to directly calculate the permittivities as well as Drude plasma frequencies. Based on the observations on figures 1(c) and (d), strains along the x axis seem to strongly change the shape of the conduction bands, while strains along the y axis have much smaller impacts; this may be attributed to the crystal structure of phosphorene. One might also notice that the gap is significantly underestimated compared with the experimental value which is around 2.0 eV; as long as we are concerned about the plasmons at low frequencies, this would not be a serious issue. The vertical axes in figures 1(c) and (d) are plotted as E − E F , where E F is the Fermi energy, the exact value of which slightly depends on calculation details. In the calculation of the bands, Fermi-Dirac distribution of the electrons with proper broadening (0.0019 Ry, corresponding to 300 K) has been used. The structure of phosphorene is still theoretically stable with strain as large as 10%.

Band structures and electron masses
Although we intend to directly apply DFT other than the effective mass approximation to derive the plasmonic responses, it would be interesting to extract the effective mass as a function of strain, since this is an important parameter for both calculations and experiments. We have extracted such parameter based on the fitting of the band energies and Bloch wave vectors. Given the asymmetry in the band structure between the conduction and valence bands, only electron doping can lead to significant unusual plasmonic responses in bP; thus, only electron mass has been investigated here. The significant modification of the curvature near the bottom of the conduction band at Γ under different strains as shown in figure 1(c) has already clarified this point. Figures 2(a)/2(b) shows the electron masses along the x (denoted as m xx , plotted with black lines and scatters) and y (denoted as m yy , plotted with red lines and scatters) axes as functions of strains along the zigzag/armchair direction. We have done such fitting near the Γ point (plotted using solid lines) as well as the CBM (plotted using dashed lines). Under 10% strain along the x and 7%-10% strains along the y directions, bP becomes an indirect-bandgap material and the effective mass changes significantly, as shown by the shaded regions in figures 2(a) and (b). Due to the limited numerical accuracy in DFT, quadratic fitting in the shaded area can only give us rough estimating of the effective mass. The significant changing of the electron mass near 4% in figure 2(a) is caused by the inversion of the first and the second conduction bands as mentioned above, which can explain the so-called reversed hyperbolic responses in the low electron doping limit with strains beyond 5% along the x/zigzag direction [28].

Permittivities and Drude plasma frequencies
The permittivity of phosphorene consists of two parts: one is from the intraband transition and the other is from the interband transition. The intraband part of the permittivity can be written as ε intra = −ω 2 D / [ω(ω + iγ)], where γ is the band broadening and ω D is the Drude plasma frequency which can be written as where V is the cell volume, f is the Fermi-Dirac function, ϕ nk is the wavefunction with band index n and momentum k,q is the unit vector denoting the direction, v is the velocity operator; the interband part is quite complex and can be written as where the new symbol η is the band broadening [38]. Two expressions of permittivities above use independent particle approximation (IPA) which works quite well for metals since the electron-hole interaction is screened by conduction electrons. This corresponds to phosphorene under electron doping in this paper. Other approaches based on many-body perturbation theory such as GW and BSE are intended to explore the quasiparticle and exciton effects on optical properties, and usually applied in dielectric materials other than metals; moreover, these two methods require a lot of computational effort and unlikely applicable for investigations in this paper. However, we intend to consider the many-body effects on plasmons in bP in the near future. Since the Drude plasma frequency basically determines the plasmonic responses and it only involves states near the Fermi energy, we believe IPA is accurate enough. The losses, for example caused by scattering, are taken into account phenomenonlogically by two parameters γ and η. The interband transition of electrons corresponds to vertical electron-hole excitation, which just means optical losses if one only cares about the collective plasmonic excitation. Studies show that the interband transition part of the permittivity is more like background, which is essential to in-plane hyperbolicity and cannot be omitted if one concerns about the plasmonic responses in the hyperbolic regime. Let us first focus on phosphorene without any strain. The real parts of the total permittivities include both the intraband and interband components are plotted in figure 3(a) under various doping levels, where the numbers in the figure indicate the up-shift of the Fermi energy in eV with respect to the baseline which is 0.5 eV above the energy of the highest occupied state. In this paper, the Fermi energy has been shifted upwards by as large as 2 eV. This might be unreal since, for example, applying backgate voltages usually cannot create such high level of doping. However, we are interested in the unusual plasmonic responses in bP under a little extreme strain and doping conditions; and furthermore, strains can lower the doping threshold needed for the occurrence of such phenomena. Solid and dashed curves in figure 3(a) respectively correspond to the permittivities along the x and y axes. Small peaks on the curves come from the interband transition of electrons. When the Fermi energy is located within the gap, phosphorene is a dielectric material with positive permittivity. With increasing level of electron doping, the Fermi energy will meet the conduction bands; intraband transition of electrons will then happen and the permittivities become negative which indicates that now phosphorene is a metallic material. In the low frequency end, the intraband part of the permittivity dominates which is mainly determined by the Drude plasma frequency denoted as ω D . ω D is a complicated function of the band energies and wavefunctions. Since it differs along the x and y directions, we respectively denote them as ω D,xx and ω D,yy . Further observations on the real parts of the permittivities in figure 3(a) lead to such conclusion: ω D,yy > ω D,xx for low Fermi energies, while ω D,yy < ω D,xx for high Fermi energies. During our calculations, significant unusual plasmonic responses have been found in bP only with electron doping, i.e. with positive Fermi energy; hole doping corresponding to negative Fermi energy commonly cannot lead to such phenomena especially in the low doping limit, given the asymmetry in the band structures between the conduction and valence bands as shown in figures 1(c) and (d). Based on the expression of ω D above, the Drude plasma frequency relates to the velocity of electrons and states near the Fermi energy. We have directly extracted the Drude plasma frequencies along two directions from the DFT calculations. Figure 3 In the followings, we investigate the plasmonic responses in phosphorene with both strain and doping applied. We have calculated the interband and intraband parts of the permittivities seperately under various Fermi energies and strains, and then extracted the Drude plasma frequencies. Strains from 0% to 10% have been applied either along the x or the y axis. Since undoped phosphorene is a semiconductor, the energy of 0.5 eV above the highest occupied state has been set to zero; and the Fermi energy has been shifted upwards with a step as small as 0.05 eV. Figures 4(a)   volume-averaged, which may differ from their actual values. Since Drude plasma frequency is a complicated function of the band energies and wavefunctions, no simple relation can be found especially in the situation where several conduction bands are involved. Based on its expression, ω 2 D q = (4π/V) nk | ϕ nk |q.v |ϕ nk | 2 (−∂f /∂E), it is natural that the Drude plasma frequency is not monotonic in Fermi level because it is related to the states near the Fermi energy and more than one bands are involved. As the Fermi energy increases, the contribution from one particular band may become less, leading to slightly decreasing of the Drude plasma frequency which certainly cannot be explained by the change on the electron mass. In the low-doping limit where the effective mass approximation is valid, the Drude plasma frequency should be monotonic in Fermi energy, then the electron mass can be used to estimate the plasmonic responses in bP. For example, with 5% strain along the x axis, the change on the electron mass as shown in figure 2(a) can clearly explain the occurrence of significant unusual plasmonic response with small Fermi energies. Now, we focus on the difference between the plasma frequencies along two directions: ω D,xx − ω D,yy = Δω D . Usually, Δω D is negative because ω D is higher along the y axis, as indicated by figure 3(b) and also figure 3(a) where most of the solid curves are above the corresponding dashed ones. With relatively high Fermi energy and large strain, Δω D can be positive, leading to unusual plasmonic responses.
We plot Δω D as functions of strain and Fermi energy in figure 5: if condition Δω D > 0.15 eV is satisfied, the corresponding area is filled with yellow otherwise dark blue. Such criterion (0.15 eV) is chosen mainly for illustration purposes, with which one can clearly see the boundaries between usual and unusual plasmonic responses, as shown in figure 5. We chose such value because Δω D > 0.15 eV commonly corresponds to significant unusual plasmonic responses based on our calculations. Figures 5(a) and (b) respectively show the Δω D when the strains are along the x and y directions. Clearly, strains along the x axis can strongly affect the Drude plasma frequencies, leading to unusual plasmonic responses even in the low-doping limit. White dashed line in figure 5(a) roughly denotes the boundary between the so-called usual and unusual plasmonic responses. Drude plasma frequency ω D is directly related to the intraband transition part of the permittivity; Δω D > 0.15 eV means that in the low frequency end the intraband transition induced permittivity along the x direction (denoted as ε intra,xx ) is significantly below the one along the y direction (denoted as ε intra,yy ), i.e. Re(ε intra,xx ) < Re(ε intra,yy ). To further confirm this, we plot only the intraband part of the permittivities in figure 6. For simplicity, only the real parts are shown. Solid and dashed curves denote Re ε intra,xx and Re ε intra,yy , respectively. Figures 6(a) and (b) plot Re ε intra,xx together with Re ε intra,yy under 5% and 10% strains along the x axis, respectively; and figures 6(c) and (d) plot the ones under 5% and 10% strains along the y axis, respectively. Indeed, there are cases where the dashed curves are above the corresponding solid ones, indicating ω D,xx > ω D,yy . Intraband transition part of the permittivity is always negative (neglecting the vacuum background), as shown in figure 6, which can only explain the plasmonic responses in the elliptic regime. Regarding the interband transition part of the permittivity as positive background, the total permittivity which include both the intraband and interband transitions can have hyperbolic regions. Based on the discussions above, it is possible to create significant  unusual elliptic as well as hyperbolic plasmonic dispersions with proper strain and doping. To prove such judgement, phosphorene without and with 5% strain along the x axis have been selected as examples in the following subsection.  the hyperbolic regime. This is due to the fact that higher Drude plasma frequency occurs along the y direction(the armchair direction), and at low frequencies the blue curves are above the red ones. Since the intraband part of the permittivity is always negative, the interband part of the permittivity can be regarded as positive background, making the total permittivity along one particular direction (along the zigzag direction in the unstrained case) positive within some frequency ranges.

Elliptic and hyperbolic plasmonic responses with topological transition
With 5% strain along the x axis, things are total changed. Significantly higher Drude plasma frequency occurs along the zigzag direction, other than the armchair direction as mentioned in the above paragraph. The elliptic as well as the hyperbolic regimes can also be found based on the signs of the permittivities. Figures 7(c) and (d) are respectively the total permittivities of the strained phosphorene (5% along the x axis) with 0.5 eV and 0.75 eV shifted Fermi energies. The elliptic and hyperbolic regions are also shaded with light red and blue, respectively. Compared with figures 7(a) and (b), the solid red curves representing the permittivities along the y axis are now above the blue ones which represent permittivities along the x axis, meaning that the anisotropy of the plasmonic responses has been reversed. Although the interband transition induced permittivity is small, it cannot be neglected. The hyperbolic plasmonic responses seem to be robust against strains; with large enough strain and doping, the reversed hyperbolicity may occur, where the permittivity along the zigzag/armchair direction is negative/positive, as shown by the region shaded with light blue color in figure 7(d).
In order to derive the dispersion relations of the plasmonic modes in the extended phosphorene monolayer, the total permittivities have been converted into 2D conductivities using a simple model. The real and imaginary parts of the permittivity are denoted as ε 1 and ε 2 , respectively; the conductivity with σ 1 and σ 2 respectively as its real and imaginary parts can be derived from the permittivity as σ 2 = (1 − ε 1 ) ωε 0 t and σ 1 = ωε 0 tε 2 , where ω is the angular frequency, ε 0 is the permittivity of vacuum, and t is the thickness. The real part of the conductivity is always positive in passive materials, corresponding to optical losses; however the imaginary part can be positive or negative depending on whether the material is metallic or dielectric. For simplicity, we neglect the losses and consider only the imaginary parts of the conductivities. The conductivity tensor is diagonal and can be written as σ = diag[σ xx , σ yy ]. Due to the anisotropy, the plasmonic modes consist of both transverse-electric and transverse-magnetic polarized components. After matching the boundary conditions at two sides of the conductive sheet, one can derive the equation for the dispersion relations of plasmons in the monolayer [33]: where β = (β x , β y ) is the in-plane wave vector, k 0 is the free-space wave number, γ = β 2 x + β 2 y − k 2 0 , Z 0 is the impedance of vacuum. Figure 8(a) plots the σ 2 of unstrained phosphorene under various Fermi energies where solid and dashed curves respectively denote the ones along the x and y axes. As the photon energy increases, σ 2 is firstly positive, and then negative; due to the anisotropy, there exist frequency ranges where σ 2 along two directions have opposite signs, corresponding to the hyperbolic plasmonic responses. The isofrequency contours or k-surfaces are plotted in figure 8(b). The topology of the k-surfaces changes from hyperbolic to elliptic as the Fermi energy increases. The photon energy has been selected so that the transition of the topology of the plasmonic responses can be clearly observed. Actively switching the topology of the plasmonic responses between elliptic and hyperbolic regimes via doping is possible in unstrained bP. As for the strained phosphorene (5% along the x axis), the σ 2 is plotted in figure 8(c). Hyperbolic plasmonic responses in 2D materials are results of extreme anisotropy. The permittivity is negative along one direction but positive along the perpendicular one, i.e. the material is metallic only along one particular axis. The isofrequency contour or the k-surface of the plasmonic mode in the extended monolayer is a hyperbola. Due to the anisotropy of the conductivity, these modes cannot be completely separated into transverse-electric and transverse-magnetic polarized components; in most cases, they consist of both polarizations. Using 2D conductivity instead of permittivity, one can analytically solve the dispersion relations. By plotting the isofrequency contours, the topological signature of the plasmonic responses, either elliptic or hyperbolic, can be clearly demonstrated. As for phosphorene, the hyperbolicity is due to the fact that the Drude plasma frequencies caused by intraband transition as well as the positive background permittivities caused by interband transition largely differ along zigzag and armchair directions. There exists a narrow frequency range in the infrared where the conductivities and the permittivities along two directions have different signs. Since the Drude plasma frequency can be drastically modified by changing the Fermi level under specific strain, active switching the topology of the plasmonic responses via doping is certainly possible in phosphorene.

A simple model for the in-plane permittivities
In the end, we point out that at low frequencies the interband transition part of the permittivity can be simply replaced by a constant which can be regarded as background. Having this in mind, the total permittivity could be written as where i = xx or yy, and ε r,i is the constant representing the interband transition induced permittivity, the exact value of which depends on the fitting of the total permittivities from the model to those from DFT. This is the only parameter in the model that needs to be determined. If one only cares about the plasmonic responses in the infrared, such model is sufficient. We have done the fitting for bP with all the strains and Fermi energies in this paper; during the procedure, we have minimized the deviations between data from With these parameters, one can directly calculate the total permittivities using equation (2) without deploying any time-consuming DFT calculations. For simplicity, we show here only two of the fitting results: one is the unstrained phosphorene with 0.8 eV shifted Fermi energy, as shown in figure 9(a); the other is the strained phosphorene (5% along the x axis) with 0.6 eV shifted Fermi energy, as shown in figure 9(b). Solid curves represent the data from the model and open circles the data from the DFT calculations. At low frequencies, small oscillations of the permittivity induced by the interband transition of electrons are not important (as shown by figure 9(a)), and such model is accurate enough. In the low-doping limit, no peaks occur on the curve (as shown by figure 9(b)), such model can very well describe the total permittivities. Based on figure 9, one can see that the model can describe the total permittivities both in the elliptic as well as the hyperbolic regimes.

Methods
Quantum espresso package was used in the DFT calculations [34]. To avoid inter-layer interactions, the total thickness of the cell was 30 angstroms, leaving enough vacuum at both sides of the monolayer. Optimized norm-conserving potentials were used together with Perdew-Burke-Ernzerhof functionals [35,36]. The van der Waals correction was DFT-D2. The structure was fully relaxed before any calculations of the bands and permittivities. The coulomb potential was truncated to represent two-dimensional system [37]. In the calculations of the bands, a 21 × 15 × 1 k-grid was used; the cut-off energies of wavefunctions and densities were 120 and 480 Ry, respectively. Fermi-Dirac function with 0.0019 Ry (about 300 K) broadening was used for electron occupation, simulating room-temperature excitation of electrons on bands. Other thermal effects such as electron-phonon interactions as well as electron-hole pairs were not considered. Bands were reconstructed using the optimal basis (OB) function approach, which were then used to calculate the permittivities [38]. We increased the cut-off energies of the wavefunctions and densities to 150 and 600 Ry, respectively. There were four atoms in the cell with total 40 bands, ten of which were valence bands. After calculations of the density, a uniform coarse grid was used in the generation of OB functions, and the permittivities were calculated on a 40 × 30 × 1 grid. The commutator of the non-local part of the pseudopotential was also included. The threshold for the construction of the optimal bases was 0.005. The broadenings for intraband and interband transitions were both 5 meV which were used to simulate the losses caused by, for example electron-phonon scattering due to thermal excitation, and were just empirical parameters here. Slightly changing these two parameters would not affect the main results in this paper. The smearing for the Fermi energy and Drude plasma frequencies were 25 meV, corresponding to 300 K. Spin-orbit coupling was not considered.

Conclusion
In this paper, we have systematically investigated the plasmonic responses in phosphorene under various strain and doping conditions. Strains can significantly modify the shape of the electron bands, which finally change the Drude plasma frequencies along zigzag and armchair directions. Under specific condition, significantly higher Drude plasma frequency occurs along the zigzag direction other than the armchair direction as in the usual case. The resulting plasmonic responses change their anisotropy; even the so-called reversed hyperbolic responses may be observed. Further investigations indicate that actively switching the topology of the plasmonic responses between elliptic and hyperbolic regimes is possible at specific frequencies. By replacing the interband transition part of the permittivity by a constant, a simple model has been proposed which can give accurate estimations of the permittivities at low frequencies. Our studies may extend the existing investigations of bP plasmons and lead to band engineering, for example via strain and doping, as a way to control plasmons in two-dimensional materials.