Tuning 2D hyperbolic plasmons in black phosphorus

TL;DR

This paper demonstrates that atomically thin black phosphorus can be engineered as a hyperbolic material across a broad spectrum, with tunable plasmonic properties influenced by gain and loss mechanisms.

Contribution

It introduces the concept of tuning 2D hyperbolic plasmons in black phosphorus through engineering and gain, expanding its potential for plasmonic applications.

Findings

Black phosphorus can operate as a hyperbolic material from visible to ultraviolet spectrum.

Introduction of optical gain creates a new hyperbolic region in the infrared.

Hyperbolic plasmon characteristics depend on gain-loss interplay along crystalline directions.

Abstract

Black phosphorus presents a very anisotropic crystal structure, making it a potential candidate for hyperbolic plasmonics, characterized by a permittivity tensor where one of the principal components is metallic and the other dielectric. Here we demonstrate that atomically thin black phosphorus can be engineered to be a hyperbolic material operating in a broad range of the electromagnetic spectrum from the entire visible spectrum to ultraviolet. With the introduction of an optical gain, a new hyperbolic region emerges in the infrared. The character of this hyperbolic plasmon depends on the interplay between gain and loss along the two crystalline directions.

Figures

Equations

$$\mathcal{H}=\sum _{i\ne j}{t}_{ij}{c}_i^{†}{c}_j+\sum _{i\ne j}{t}_{p,ij}{c}_i^{†}{c}_j,$$

$$t_{ij}({\mathbf{r}}_{ij})=t_{ij}({\mathbf{r}}_{ij}^0)\left(1-{\beta }_{ij}\frac{|{\mathbf{r}}_{ij}-{\mathbf{r}}_{ij}^0|}{|{\mathbf{r}}_{ij}^0|}\right),$$

$$\text{Re}({\sigma }_{\alpha \alpha }(\omega ))=-\frac{g_S}{\Omega \omega }{\int }_{BZ}\text{ Im}[\sum _{i,j}|⟨\mathbf{k}i|J_{\mathbf{k}\alpha }|\mathbf{k}j⟩{|}^2$$

$$\times \frac{f(E_{\mathbf{k}i}-\mu )-f(E_{\mathbf{k}j}-\mu )}{E_{\mathbf{k}i}-E_{\mathbf{k}j}+\omega +i\delta }]d^2\mathbf{k}.$$

$$J_{\mathbf{k}\alpha }=-\frac{ie}{\hslash }\sum _{i,j}{e}^{i({\mathbf{r}}_j-{\mathbf{r}}_i{)}_{\alpha }\cdot \mathbf{k}}{t}_{ij}({\mathbf{r}}_j-{\mathbf{r}}_i{)}_{\alpha }{c}_{\mathbf{k}i}^{†}{c}_{\mathbf{k}j}^{}.$$

$$f(E-\mu )=\frac{1}{e^{(E-\mu )/kT}+1},$$

$$n_F(E)=\theta (E)f(E+\frac{E_g}{2}+\Delta \mu )+\theta (-E)f(E-\frac{E_g}{2}-\Delta \mu ),$$

$$\text{Im}({\sigma }_{\alpha \alpha }(\omega ))=-\frac{2\omega }{\pi }\mathcal{P}{\int }_0^{\infty }\frac{\text{Re}({\sigma }_{\alpha \alpha }({\omega }^{\prime }))}{{\omega }^{\prime 2}-{\omega }^2}d{\omega }^{\prime }.$$

$$\text{Im}({\sigma }_{xx}(\omega ))\times \text{Im}({\sigma }_{yy}(\omega ))<0.$$

$${\sigma }_{F,yy}\sim \frac{(q_F{\Gamma }_{res}/2+\omega -{\omega }_{res}{)}^2}{({\Gamma }_{res}/2{)}^2+(\omega -{\omega }_{res}{)}^2}$$

$$q_{\parallel }^2={\gamma }^2+k_0^2,$$

$$\gamma$$

$$+\sqrt{{\left(\frac{2}{{\eta }_0}+\frac{{\eta }_0}{2}{\sigma }_{xx}{\sigma }_{yy}\right)}^2-4{\sigma }_{\mathbf{qq}}{\sigma }_{\perp \perp }}],$$

$${\sigma }_{\mathbf{qq}}(\chi )={\sigma }_{xx}{\cos }^2\chi +{\sigma }_{yy}{\sin }^2\chi ,$$

$${\sigma }_{\perp \perp }(\chi )={\sigma }_{xx}{\sin }^2\chi +{\sigma }_{yy}{\cos }^2\chi .$$

$$\tan {\chi }_0=\sqrt{\left|\frac{\mathrm{Im}({\sigma }_{xx})}{\mathrm{Im}({\sigma }_{yy})}\right|}.$$

$${\sigma }_{xx}=-{\sigma }_{yy}\tan {\chi }_0.$$

Full text

Introduction. Semiconducting two dimensional (2D) crystals are excellent platforms for tuneable optoelectronics, thanks to their remarkable response to external electrical and mechanical stimuli Avouris et al. (2017); Roldán et al. (2017). In particular, atomically thin black phosphorus Li et al. (2014); Liu et al. (2014); Castellanos-Gomez et al. (2014); Xia et al. (2014) (BP) has shown extraordinary tuneability of its optical and electronic properties by several methods Roldán and Castellanos-Gomez (2017), like electrostatic gating Lin et al. (2016); Peng et al. (2017); Whitney et al. (2016); Deng et al. (2017); Liu et al. (2017), chemical functionalization Kim et al. (2015), quantum confinement (number of layers) Yang et al. (2015), external strain Quereda et al. (2016) or high pressure Rodin et al. (2014); Xiang et al. (2015). This allows the control of light-matter interaction in these materials, in particular the dispersion of collective polaritonic excitations Low et al. (2017).

Apart from being a highly tuneable optoelectronic crystal, the lattice structure of BP is very anisotropic Avouris et al. (2017); Keyes (1953). The in-plane anisotropy implies optical birefringence, of which the extreme limit would be hyperbolicity, where the permittivity tensor has principal components of opposite sign Poddubny et al. (2013); Yermakov et al. (2015); Gomez-Diaz et al. (2015); Nemilentsau et al. (2016). Recently, in-plane hyperbolicity was implemented experimentally in the GHz frequency range using a metallic metasurface Yermakov et al. (2018). Moreover, in-plane hyperbolicity in natural van der Waals material $\alpha$-MoO3 was reported and existence of hyperbolic surface polaritons was experimentally verified Ma et al. (2018); Zheng et al. (2018). The strong anisotropy of BP suggests its potential as a natural hyperbolic material, offering new possibilities for actively manipulating polaritons in 2D, such as directional plasmons, light emitters, superlensing effects, Nemilentsau et al. (2016); Correas-Serrano et al. (2016) etc. In this Letter, we discuss the possiblity of driving BP into the hyperbolic region via electrostatic tuning, strain or layers number. We demonstrate that atomically thin BP can be efficiently tuned to become hyperbolic in a broad spectral range from the entire visible spectrum to the ultraviolet. In addition, the presence of a bandgap in excess of the optical phonon energy lends itself as a possible 2D semiconductor gain medium Shang et al. (2017); Wu et al. (2015). With the introduction of a population inversion, we show that optical gain results in a new hyperbolic region in the infrared. Finally, we study the behavior of plasmons in both of these hyperbolic regions.

Optical Conductivity and Band Model. We describe BP by means of a $p_{z}$-orbitals tight-binding model fitted to ab initio $GW$ methods Rudenko and Katsnelson (2014); Rudenko et al. (2015).

[TABLE]

where $c_{i}^{\dagger}$ ($c_{i}$) creates (annihilates) an electron at site $i$, and ten intra-layer $t_{ij}$ and five inter-layer $t_{p,ij}$ hopping terms are considered in the model. The obtained band structure corresponds to an anisotropic direct band gap semiconductor, with the gap at the $\Gamma$ point of the Brillouin zone.

The model can be straightforwardly extended to incorporate arbitrary electrostatic and strain fields. An electric field is applied by modifying the on-site potentials with $\epsilon_{i}=e\times\Delta U\times z_{i}$, where $e$ is the elementary charge, $\Delta U$ is the bias voltage (with units V/nm) and $z_{i}$ is the $z$-coordinate of site $i$.

On the other hand, the application of external strain leads to a variation of the interatomic bond lengths, which further modifies the hopping terms as Suzuura and Ando (2002)

[TABLE]

where $|{\bf r}_{ij}^{0}|$ is the distance in the equilibrium positions between two atoms $i$ and $j$, $|{\bf r}_{ij}|$ the distance in the presence of strain, and $\beta_{ij}=-d\ln t_{ij}(r)/d\ln(r)|_{r=|{\bf r}_{ij}^{0}|}$ is the dimensionless local electron-phonon coupling. A microscopic estimation of $\beta_{ij}$ can be done based on the direct comparison between ab initio and tight-binding calculations. Here we use $\beta\approx 4.5$ because this value was proved to give a matching between ab initio and tight-binding calculations for the direct-to-indirect bandgap transition under uniaxial strain San-Jose et al. (2016). The mechanical properties of BP are highly anisotropic, with zigzag direction being about four times stiffer than armchair direction Wei and Peng (2014); Quereda et al. (2016). Therefore we use uniaxial strain along the armchair direction for our calculations, accounted for by the strain tensor $\mathbf{\epsilon}_{\mathrm{AC}}=\epsilon_{yy}\mathrm{diag}(\nu^{\mathrm{AC}}_{z},1,-\nu^{\mathrm{AC}}_{z})$, where the Poisson ratios $\nu$ are estimated to be $\nu^{\mathrm{AC}}_{x}\approx 0.2$, and $\nu^{\mathrm{AC}}_{z}\approx 0.2$ Wei and Peng (2014). We notice the importance to consider the out-of-plane Poisson ratio $\nu_{z}$ in our calculations, that accounts for the widening (flattening) of the lattice under compressive (tensile) strain.

The optical conductivity in the zigzag ($\sigma_{xx}$) and armchair ($\sigma_{yy}$) directions is given by the Kubo formula

[TABLE]

Here, $|\mathbf{k}i\rangle$ and $E_{\mathbf{k}i}$ are the eigenstates and eigenenergies for momentum $\mathbf{k}$ and orbital $i$. $g_{S}=2$ is the spin degeneracy, $\Omega$ is the unit cell surface, and $\delta=5$ meV is a small damping parameter. $J_{\mathbf{k}\alpha}$ is the current operator in the $\alpha$-direction

[TABLE]

Moreover, $f(E-\mu)$ is the Fermi-Dirac distribution with Fermi level $\mu$:

[TABLE]

where we set $T=300$ K.

We can replace $f(E)$ with a quasi-equilibrium distribution $n_{F}(E,\Delta\mu)$ to introduce population inversion, which will produce optical gain Low et al. (2018); Page et al. (2015):

[TABLE]

where $E_{g}$ is the band gap.

Using the Kramers-Kronig relations, we can also obtain the imaginary part

[TABLE]

Hyperbolic Regions. We first define the condition for hyperbolicity. We note that the real part of the dielectric permittivity is proportional to $\text{Im}(\sigma)$, a consequence of current continuity. Then, a hyperbolic region appears when

[TABLE]

On the other hand, $\text{Re}(\sigma)$, is directly proportional to the optical absorption of the free-standing 2D layer. For pristine bilayer BP, the optical conductivity components are plotted in Fig. 1. The first thing we observe is that the peculiar puckered structure of BP leads to a strong linear dichroism, i.e., a large difference in optical conductivity for incident polarized light along armchair and zigzag directions Low et al. (2014). For a bilayer sample, its optical absorption revealed two sharp peaks along the armchair direction, due to the two interband excitations indicated in red in the band structure. These resonant-like features, for light polarized along armchair, has also been observed experimentally Zhang et al. (2017). On the contrary, light polarized along zigzag shows a featureless monotonically increasing optical absorption instead. Whereas $\text{Im}(\sigma_{xx})$ is negative throughout the spectrum, $\text{Im}(\sigma_{yy})$ goes from negative to positive around $\omega_{h}=2.8$ eV, which results in a hyperbolic region starting at that frequency.

The sign change in $\text{Im}(\sigma_{yy})$ along the armchair direction is key to the appearance of the hyperbolic region as indicated in Fig. 1. This can be traced to the resonant-like feature in the optical absorption $\text{Re}(\sigma_{yy})$ at $\omega_{res}=2.6$ eV. The spectral shape of $\text{Re}(\sigma_{yy})$ can be described by a Fano resonance curve

[TABLE]

and setting $\sigma_{F,xx}\sim\omega^{n}$. A fit of these curves to the region around the second peak in the optical conductivity of bilayer BP is shown by the dashed lines in Fig. 1. Its Kramers-Kronig pair, which correponds to $\text{Im}(\sigma_{yy})$, reveals a sign change after $\omega_{res}$. Hence, we can attribute the origin of hyperbolicity to the strong and anisotropic resonant-like interband transitions.

If we introduce population inversion (Eq. 6, Fig. 2), however, the situation becomes qualitatively different: a new hyperbolic region appears in the infrared range. Here, we assume that the quasi-Fermi levels are such that $\mu_{h}=−\mu_{e}$, and that the electron and hole baths can be described by a common temperature (see Fig. 2c). Optical pumping Ni et al. (2016); Lui et al. (2010), where electrons and holes are generated in pairs, of a charge neutral system with particle-hole symmetry would fit such scenario. The optical gain causes $\text{Re}(\sigma_{yy})$ to become negative. The spectral window where $\text{Re}(\sigma_{yy})<0$ roughly coincides with $E_{g}<\omega<E_{g}+2\Delta\mu$, where there optical transitions between the population inverted electron and hole bands are allowed. In the region up to $\omega=1.27$ eV, we find that $\text{Im}(\sigma_{yy})>0$. In the zigzag direction, the real part of the conductivity also flips sign between $E_{g}<\omega<E_{g}+2\Delta\mu$. The imaginary part in this direction, on the other hand, remains negative throughout the entire frequency range, causing a new hyperbolic region in the infrared.

Since the origin of the hyperbolicity is related to the strong resonant-like anisotropic interband absorption between the largest conduction and valence subband indices, one expects that the spectral range of hyperbolicity can be tuned with bandstructure engineering. Indeed, the onset of the hyperbolic region $\omega_{h}$ can be tuned with the number of layers, strain, bias and doping (Fig. 3). $\omega_{h}$ goes down for compressive strain, because the peaks in the optical spectrum move to lower frequencies Quereda et al. (2016). For increasing bias in the bilayer case, the hyperbolic onset frequency $\omega_{h}$ first goes up, because the bands corresponding to the excitation causing the second peak move away from one another. Then it goes down as the band gap closes and a new peak appears in between the two existing peaks, because breaking the mirror symmetry in the $z$-direction allows for new hybrid transitions Zhang et al. (2017). For more details about how the optical conductivity is affected by these parameters, we refer to the supplemental material. Moreover, $\omega_{h}$ goes up for an increasing number of layers, because extra layers add peaks to the optical conductivity in the armchair direction, and the hyperbolic region appears after the last peak. Finally, for increased doping, $\omega_{h}$ moves further up, as the first peak becomes less prominent due to Pauli blocking.

Hyperbolic Plasmons. Finally, let us consider characteristics of plasmons that can be supported by an hyperbolic material. We assume that the plasmon propagates at an angle $\chi$ with respect to the $x$-axis, where the plasmon wavector takes the form $\mathbf{q}=q_{x}\mathbf{e}_{x}+q_{y}\mathbf{e}_{y}+i\mathbf{e}_{z}\gamma$, where $q_{x}=q_{\parallel}\cos\chi$, $q_{y}=q_{\parallel}\sin\chi$, $q_{\parallel}=\sqrt{q_{x}^{2}+q_{y}^{2}}$. The dispersion relation for the hyperbolic plasmon takes form

[TABLE]

where

[TABLE]

$k_{0}^{2}=\omega^{2}\mu_{0}\varepsilon_{0}$, and

[TABLE]

The iso-frequency contours ($\omega(q_{x},q_{y})=\mathrm{constant}$), calculated using Eq. (10), are presented in Fig. 4. We consider the cases of BP with gain ($\sigma_{xx}=(0.07-i6.4)\times 10^{-2}\sigma_{0}$, $\sigma_{yy}=(-0.03+i0.60)\sigma_{0}$) and BP under strain ($\sigma_{xx}=(0.20-i0.31)\sigma_{0}$, $\sigma_{yy}=(0.20+i0.54)\sigma_{0}$). It can be seen from Figs. 4a,b, that only in the case of BP with gain the iso-frequency contour resembles a hyperbola with the asymptotes defined as

[TABLE]

For BP under strain, the iso-frequency contour resembles a figure eight shape, even though the hyperbolicity condition (8) is met.

This behavior stems from the fact that the hyperbolic shape of the iso-frequency contour is related to the poles of the denominator in Eq. (Tuning 2D hyperbolic plasmons in black phosphorus), defined by zeros of $\sigma_{\mathbf{q}\mathbf{q}}$ (i.e. when $\sigma_{\mathbf{qq}}\to 0$, then $|\gamma|\to\infty$ and $|q_{\perp}|\to\infty$). In particular, it is straightforward to demonstrate that in the case of a purely imaginary conducitivity tensor (i.e. no losses or gain) the condition $\sigma_{\mathbf{q}\mathbf{q}}=0$ leads to Eq. (14) for hyperbola asymptotes.

If the components of the conductivity tensor are both lossy ($\mathrm{Re}(\sigma_{xx,yy})>0$) or both have gain ($\mathrm{Re}(\sigma_{xx,yy})<0$), the module of the conductivity, $|\sigma_{\mathbf{qq}}|$, is never zero. In fact, for the hyperbola asymptote angle $\chi_{0}$, $\left|\sigma_{\mathbf{qq}}(\chi_{0})\right|=\left|\mathrm{Re}(\sigma_{xx})\right|\cos^{2}\chi_{0}+\left|\mathrm{Re}(\sigma_{yy})\right|\sin^{2}\chi_{0}$. Thus, $\left|\sigma_{\mathbf{qq}}(\chi_{0})\right|$ increases with the increase of $\left|\mathrm{Re}(\sigma_{xx,yy})\right|$ both for the lossy material and the material with gain. This leads to the decrease of $q_{\parallel}(\chi_{0})$ and, eventually, destroys hyperbolicity when either losses or gain are too high. For example, this is the case in BP with strain presented in Figs. 4b,d. The high losses of the material in the hyperbolic regime lead to the hyperbola folding into a figure eight shape iso-frequency contour. Moreover, the plasmons itself are very lossy in this case as is quantified by ratio $\mathrm{Re}q_{\parallel}/\mathrm{Im}q_{\parallel}$ in Fig. 4d.

The case where one of the components of the conductivity tensor is lossy ($\mathrm{Re}(\sigma_{xx})>0$, while the other has gain ($\mathrm{Re}(\sigma_{yy})<0$), requires separate consideration. In this case, $\sigma_{\mathbf{qq}}(\chi_{0})=\left|\mathrm{Re}(\sigma_{xx})\right|\cos^{2}\chi_{0}-\left|\mathrm{Re}(\sigma_{yy})\right|\sin^{2}\chi_{0}$. When both the losses and the gain are small, then $\left|\sigma_{\mathbf{qq}}(\chi_{0})\right|$ is small as well which allows for the iso-frequency contour to preserve the hyperbolic shape, as is the case for BP with gain presented in Fig. 4a,c. This is, however, a rather trivial case which can be observed in pure lossy materials when the losses are small Nemilentsau et al. (2016). A non-trivial property of a material with gain is that $\sigma_{\mathbf{qq}}=0$, when $\tan\chi_{0}=\sqrt{\left|\mathrm{Re}(\sigma_{xx})\right|/\left|\mathrm{Re}(\sigma_{yy})\right|}$. That condition, together with Eq. (14), indicates that the iso-frequency contour preserves its hyperbolic shape for arbitrary large losses and gain, as long as the following holds true,

[TABLE]

This criterion can be verified through its iso-frequencies contours, which we defer to the SI. Hyperbolic metasurfaces also offer the possibility of spontaneous emission rate enhancement through the Purcell effect Purcell (1946); Gjerding et al. (2017). We defer these results to the SI.

Conclusion. In conclusion, we showed in this work how BP can be made hyperbolic across a broad spectral range, and studied the influence of optical gain on the hyperbolic plasmons. The ease of electrical tuning of the plasmons would open up new opportunities for flatland nanophotonics .

Acknowledgements. RR acknowledges financial support from the Spanish MINECO through Ramón y Cajal program, Grants No. RYC-2016-20663 and Grant No. FIS2014-58445-JIN. MIK acknowledges financial support from the European Research Council Advanced Grant program (Contract No. 338957). SY acknowledges financial support from National Key R$\And$D Program of China (Grant No. 2018FYA0305800) and National Science Foundation of China (Grant No. 11774269).