Non-local quantum gain facilitates loss compensation and plasmon amplification in graphene hyperbolic metamaterials

TL;DR

This paper develops a quantum non-local conductivity model for graphene hyperbolic metamaterials, enabling loss compensation and plasmon amplification, which enhances their potential for subwavelength imaging and wave manipulation.

Contribution

It introduces a non-local quantum conductivity model combined with exact structure analysis to accurately predict hyperbolic behavior and plasmonic amplification in graphene metamaterials.

Findings

Quantum model predicts plasmonic amplification in graphene sheets.

Loss compensation achieved through chemical imbalance at finite temperature.

Bloch solutions effectively describe electric field distribution in finite structures.

Abstract

Graphene-based hyperbolic metamaterials have been predicted to transport evanescent fields with extraordinarily large vacuum wave-vectors. It is particularly at much higher wave vector values that the commonly employed descriptional models involving structure homogenization and assumptions of an approximatively local graphene conductivity start breaking down. Here, we combine a non-local quantum conductivity model of graphene with an exact mathematical treatment of the periodic structure in order to develop a tool-set for determining the hyperbolic behavior of these graphene-based hyperbolic metamaterials. The quantum conductivity model of graphene facilitates us to predict the plasmonic amplification in graphene sheets of the considered structures. This allows us to reverse the problem of Ohmic and temperature losses, making this simple yet powerful arrangement practically applicable. We analyze the electric field distribution inside of the finite structures, concluding that Bloch boundary solutions can be used to predict their behavior. With the transfer matrix method we show that at finite temperature and collision loss we can compensate for losses, restoring imaging qualities of the finite structure via an introduction of chemical imbalance.

Equations

$${\sigma }_s(q,\omega )=\frac{i\omega e^2\Pi (q,\omega )}{q^2}$$

$$\Pi \left(n\right)$$

$$+\underset{{\Pi }^{\left(h\right)}[\bar{n}]}{\underbrace{\underset{0}{\overset{\infty }{\int }}d\epsilon \left[\frac{\partial \Pi \left(q,\omega \right){|}_{\mu =\epsilon }^{T=0}}{\partial \epsilon }\bar{n}\left(\epsilon \right)\right]}}$$

$$\cos (\text{K}D)=\cos ({\text{k}}_{z,1}{d}_1)-\frac{i{Z}_0{k}_{z,1}{\sigma }_s(q,\omega )}{2{\epsilon }_1{k}_0}\sin ({\text{k}}_{z,1}{d}_1)$$

$$\frac{k_x^2}{{\epsilon }_d}+\frac{{\text{K}}^2}{{\epsilon }_{x-eff}}=k_0^2$$

$${\epsilon }_{x-eff}=\frac{{\text{K}}^{{}^{\prime }2}-{\text{K}}^{{}^{\prime \prime }2}}{k_0^2-\frac{k_x^2}{{\epsilon }_d}}-2i\frac{{\text{K}}^{{}^{\prime }}{\text{K}}^{{}^{\prime \prime }}}{k_0^2-\frac{k_x^2}{{\epsilon }_d}}$$

$$\{\begin{array}{ll}\mathrm{Re}\left[\text{K}D\right]=0,\mathrm{Im}\left[\text{K}D\right]=0& \left(a\right)\\ \mathrm{Re}\left[\text{K}D\right]=\pi ,\mathrm{Im}\left[\text{K}D\right]=0& \left(b\right)\\ \mathrm{Re}\left[\text{K}D\right]=\pi ,\mathrm{Im}\left[\text{K}D\right]=\pi & \left(c\right)\end{array}$$

$$\{\begin{array}{ll}{\sigma }_s^{{}^{\prime \prime }}(q,\omega )=\frac{2{\epsilon }_1}{Z_{0}a}\tanh \left(\frac{k_0{d}_{1}a}{2}\right)& \left(a\right)\\ {\sigma }_s^{{}^{\prime \prime }}(q,\omega )=\frac{2{\epsilon }_1}{Z_{0}a}\left(\coth \left(\frac{k_0{d}_{1}a}{2}\right)+\frac{1}{\sinh \left(\frac{k_0{d}_{1}a}{2}\right)}\right)& \left(b\right)\\ {\sigma }_s^{{}^{\prime \prime }}(q,\omega )=\frac{2{\epsilon }_1}{Z_{0}a}\frac{\cosh \left(k_0{d}_{1}a\right)+\cosh \left(\pi \right)}{\sinh \left(k_0{d}_{1}a\right)}& \left(c\right)\end{array}$$

$$\left|\mathrm{Im}\left[\text{K}D\right]\right|\le \left|\mathrm{Re}\left[\text{K}D\right]\right|$$

$$\Pi [n](q,\omega )=\frac{g}{A}\sum _{s,s^{\prime }=\pm }\sum _{\mathbf{k}}\frac{M_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{s{s}^{\prime }}[n({ϵ}_{\mathbf{k}}^s)-n({ϵ}_{\mathbf{k}+\mathbf{q}}^{s^{\prime }})]}{{ϵ}_{\mathbf{k}}^s-{ϵ}_{\mathbf{k}+\mathbf{q}}^{s^{\prime }}+\hslash \omega +i\times 0}$$

$$\Pi \left(q,\omega \right){|}_{\mu }^{T=0}=\frac{g\mu }{8\pi {\hslash }^2{\nu }_F^2}\tilde{\Pi }\left(\frac{\hslash {\nu }_{F}q}{\mu },\frac{\hslash \omega }{\mu }\right)$$

$$\cos \left(\text{K}D\right)$$

$$+\frac{i{Z}_{0}a{\sigma }_s^{{}^{\prime }}}{2{\epsilon }_1}\sinh \left(k_0{d}_{1}a\right)$$

$$\cos \left(\text{K}D\right)=\cosh \left(k_0{d}_{1}a\right)-\frac{Z_{0}a{\sigma }_s^{{}^{\prime \prime }}}{2{\epsilon }_1}\sinh \left(k_0{d}_{1}a\right)$$

$$\{\begin{array}{ll}\mathrm{Re}\left[\text{K}D\right]=\mathrm{Im}\left[\text{K}D\right]=0& \left(a\right)\\ \mathrm{Re}\left[\text{K}D\right]=\mathrm{Im}\left[\text{K}D\right]=\pi & \left(b\right)\end{array}$$

$$\{\begin{array}{ll}{\sigma }_s^{{}^{\prime \prime }}(q,\omega )=\frac{2{\epsilon }_1}{Z_{0}a}\tanh \left(\frac{k_0{d}_{1}a}{2}\right)& \left(a\right)\\ {\sigma }_s^{{}^{\prime \prime }}(q,\omega )=\frac{2{\epsilon }_1}{Z_{0}a}\frac{\cosh \left(k_0{d}_{1}a\right)+\cosh \left(\pi \right)}{\sinh \left(k_0{d}_{1}a\right)}& \left(b\right)\end{array}$$

$${\sigma }_s^{{}^{\prime \prime }}(q,\omega )=\frac{2{\epsilon }_1}{Z_{0}a}\left(\coth \left(\frac{k_0{d}_{1}a}{2}\right)+\frac{1}{\sinh \left(\frac{k_0{d}_{1}a}{2}\right)}\right)$$

Full text

I Introduction

Graphene-based hyperbolic metamaterials have been predicted to transport evanescent fields with large vacuum wave-vectors (Shung, 1986; Yan et al., 2012; Iorsh et al., 2013; Sreekanth et al., 2013; Othman et al., 2013; Poddubny et al., 2013; Ferrari et al., 2015; Chang et al., 2016; Novoselov et al., 2016). This feature has immense practical importance, as it opens up possibilities of manipulating wavevectors outside of the light cone, that have been inaccessible before (Vorobev and Tyukhtin, 2012; Pendry, 2000). This is a desirable property in fields such as imaging (Belov and Hao, 2006; Salandrino and Engheta, 2006; Wood et al., 2006; Rho et al., 2010; Andryieuski et al., 2012), control of spontaneous emission (Yang et al., 2010; Tumkur et al., 2011; Iorsh et al., 2013), nano-lithography (Xiong et al., 2009; Ferrari et al., 2015), and others (Ferrari et al., 2015; Smolyaninov and Hung, 2011). It has also been reported that due to the properties of the carrier system in graphene, the response of such metamaterials can indeed be dynamically tuned and switched between elliptic and hyperbolic regimes of operation (Iorsh et al., 2013; Othman et al., 2013; Dai et al., 2015).

Concurrently, moving deeper into the sub-wavelength region, the models, currently used to predict properties of such systems, lose their applicability: the medium homogenization approximations fail when the wavelengths are comparable to the unit cell thickness (Orlov et al., 2011, 2014); non-local effects have a significant influence on the graphene properties for the excitations with wavevectors approaching the Dirac cone (Page et al., 2015). At the same time, Ohmic losses, which are an inherent problem to all hyperbolic metamaterials (HMMs) with conductive sheets (Savelev et al., 2013), are an important issue for graphene-based HMMs (GHMMs) too. However, calculation of losses had been limited due to the approximations of the local graphene conductivity model.

In this study we apply a non-local quantum conductivity model, derived from irreducible polarizability of a particle-hole plasma to GHMMs. This facilitates exact in RPA calculations of the carrier distributions and interactions at arbitrary temperatures and levels of doping (Page et al., 2015; Pyatkovskiy, 2009; Wunsch et al., 2006; Hwang and Das Sarma, 2007). Due to the non-locality of the model, we are able to perform exact calculations of the response of graphene to the excitations with in-plane wave-vectors of energies up to the edge of the Dirac cone, which has previously been impossible at the basis of the local conductivity models. It also provides a description of graphene in the photo-inverted state, where plasmons, coupled to the excited charge carriers, will stimulate electron-hole pair recombination with consequential amplification of the triggering plasmons, uncovering a new property of graphene in the context of GHMMs (Page et al., 2015; Pyatkovskiy, 2009; Watanabe et al., 2013; Dubinov et al., 2011).

Further, we couple this conductivity model to a periodic structure description, treating such an arrangement as a photonic crystal, to calculate and identify where the metamaterial exhibits elliptic or hyperbolic properties. Due to the linearity of the Dirac Cone, the proposed band structure is of a general nature, describing simultaneously the behavior of the material under all physical values of excitation frequencies and doping levels for the separate cases of passive and active graphene. This generalized treatment eliminates the necessity to study individual regions of doping and excitation values, as it covers the full parameter space. We propose a scheme, in which we can calculate the effective optical parameters of the structure exactly, avoiding, in particular, homogenization techniques, which have a limited validity. We also provide a set of analytic tools and treatments, which could be used for an efficient characterization of graphene-based HMMs, as well as any HMMs where the conductive layers are a two-dimensional material.

We use the developed methods in order to analyze the behavior of the structure when the graphene sheets are doped, photo-inverted and at realistic temperatures. We show that in such an arrangement the stimulated plasmonic emission occurs in a frequency band of regulated width, as opposed to the line of the single graphene sheet, suspended in air. Further analysis shows that we can compensate and overcompensate Ohmic losses, opening an avenue to removing an important obstacle to a practical realization of hyperbolic metamaterials.

Finally, we calculate the quasi-static electric field distribution in finite GHMM structures, generated by a propagating narrow spatial Gaussian field envelope. In comparison with the loss-less case, we shed light on the degrading effects of collision and temperature losses on imaging of the envelope in the passive structure. Then, putting graphene sheets in the state of optical inversion, we show the length of propagation of the field through the structure being restored due to the plasmonic gain to that of the ideal system with no losses.

II Quantum Conductivity model

In order to calculate exactly the plasmon dispersion in graphene and predict plasmon amplification or attenuation, it is necessary to use a conductivity model which is both non-local and remains valid within the regions of inter-band processes. This requires an exact treatment of large wavevectors and complex frequencies. This is not achieved with the local Drude (Sreekanth et al., 2013) and Kubo (Iorsh et al., 2013; Othman et al., 2013) formulas, most often used in studies of graphene-based HMMs. They become increasingly inaccurate in the parameter space of the inter- and intra-band transitions (Page et al., 2015), limiting their applicability to scenarios of small loss and absence of gain. Therefore, in order to grasp loss compensation and amplification in graphene-based HMM, we employ a non-local quantum conductivity model, derived from the polarization function (Page et al., 2015):

[TABLE]

where $\Pi(q,\omega)$ is the propagator of the electron-hole pairs. For a general 2D electron gas in random phase approximation (RPA) it is given by the Lindhard formula (Stern, 1967; Shung, 1986) - see Appendix Eq. (9). A closed-form analytic expression for gap-less graphene at zero temperature has been derived (Pyatkovskiy, 2009), and is presented in the Appendix A as Eq. (10).

The exact Eq. (10) needs modification for the case of more complicated distribution of electrons $n\left(\varepsilon\right)$ and holes $\bar{n}\left(\varepsilon\right)$ (e.g. non-zero temperature, photo-excitation). Using the fact that the Lindhard formula is linear in the carrier distribution $n\left(\varepsilon\right)$, the polarization function of the arbitrary non-equilibrium graphene sheet can then be expressed in the form (Page et al., 2015; Page, 2016):

[TABLE]

where $\Pi^{\left(e\right)}[n]$ and $\Pi^{\left(h\right)}[\bar{n}]$ are the individual contributions of the electrons and holes. Here, $n\left(\varepsilon\right)$ and $\bar{n}\left(\varepsilon\right)=1-n\left(-\varepsilon\right)$ are arbitrary functions of the distribution of electrons and holes, allowing to incorporate chemical doping, photo-excitation and temperature (Page et al., 2018; Hamm et al., 2016).

The plasmon dispersion of graphene is calculated using the exact in RPA model for zero temperature doped graphene sheet without and with photo-excitation, and is presented as green solid lines in the Figs. 1a and 1b respectively. The excitation frequency on the vertical axis is scaled with the value of chemical imbalance $\bar{\mu}$, and the in-plane wavevector on the horizontal axis is scaled with the excitation frequency. With such scaling the light cone transforms into the vertical line at $k_{x}/k_{0}=1$ and the Dirac cone transforms into the vertical line at $k_{x}/k_{0}=c_{0}/v_{Fermi}=300$. The plots without excitation frequency scaling are shown in the Figs. 6a and 6b in Appendix A.

Chemical doping of the graphene sheet causes a shift of the Fermi level away from the Dirac point. This creates a region in phase space where inter-band processes are not possible due to insufficient excitation energies. This corresponds to the region we denote by II in the Figs. 1a and 1b. The excitations with energies above the Fermi level could decay through Landau damping into the electron-hole pairs with energies above the Fermi level. This corresponds to the region III in the Figs. 1a and 1b. Region IV in both cases corresponds to inter-band excitations (but is not of interest here, as it does not provide a mechanism for the plasmon amplification).

Optical pumping of doped graphene creates a two-component plasma of electrons and holes. This opens up a phase space within region II, where the plasmons couple to the inverted carrier plasma, trigger stimulated recombination and get amplified. We denote this area as the region I in the Fig. 1b.

The relative sizes of regions I, II and III in the active case depend on the value of doping and imbalance between the chemical potentials of electrons and holes (illustrated in the inserts of the Figs. 1a and 1b). For example, the Fermi level is at zero for intrinsic photo-excited graphene ($\mu_{e}=\mu_{h}$), so the region II is absent. For the case of extrinsic graphene with no photo-inversion ($\mu_{e}=-\mu_{h}$), the phase space for an inter-band recombination shrinks to zero and the region I disappears, becoming the case, shown in the Fig. 1a. All other possible relative values of $\mu_{e}$ and $\mu_{h}$ fall between these boundary cases.

For comparison, the Drude model, modified to include the inter-band transitions, has been used in order to calculate the plasmon dispersion in the case of active and passive doped graphene. These are shown in the Figs. 1a and 1b as dashed black curves. As can be seen, they follow the exact plasmon dispersion, calculated from the non-local quantum model in the regions of small in-plane wavevectors, but start deviating for the wavevectors, approaching the region of inter-band transitions.

III Analysis of an infinite structure

III.1 Bloch Modes

Let us first consider an infinite periodic structure, consisting of alternating graphene/dielectric layers and schematically shown in the Fig. 3b. The plasmons, excited on the outer graphene sheet, are purely evanescent in the $z-$direction, but they couple to the plasmon resonances of the subsequent graphene sheets (Das Sarma and Quinn, 1982), acquiring a phase lag upon the transition across the conductive layers (Vinogradov et al., 2010; Orlov et al., 2014). This difference in phase gives rise to a modulation of their amplitude in space. When this modulation satisfies the periodic Bloch boundary conditions, it starts tracing an effective wave, propagating down the structure, with $\operatorname{Re}\left[\text{K}\right]$, being the spatial frequency, and $\operatorname{Im}\left[\text{K}\right]$ - the rate of decay or growth of the wave in space. Due to the convention on forward direction of propagation, adopted in this work, $\operatorname{Im}\left[\text{K}\right]<0$ characterizes waves of spatially decreasing amplitude, and vice versa. The mechanism for an amplification of Bloch waves is based on the coherent amplification of the plasmonic modes due to stimulated electron-hole pair recombination of the photo-inverted graphene sheets, happening in the phase space of the region I in the Fig. 1b.

Applying the transfer-matrix method (Yeh, 2005; Yeh et al., 1977) to an infinite stack of alternating graphene/dielectric layers and imposing periodic Bloch boundary conditions, the dispersion relation of the Bloch modes can be established (Iorsh et al., 2013; Gu and Yeh, 1996; Yeh, 2005)

[TABLE]

where K is the magnitude of the Bloch wave-vector, $D$ - the period of the structure, $z$ - the axis of anisotropy, $\varepsilon_{1}$ - the permittivity of dielectric, $d_{1}$ - its thickness, $k_{0}$ - the wavevector of the incident radiation, $k_{z,1}$ - its $z$-component in the dielectric and $\sigma_{s}(q,\omega)$

• the non-local graphene sheet conductivity. Eq. (3) is an even function with two generally complex roots of opposite signs: $\pm\text{K}D$. As an example, $\operatorname{Re}\left[\text{K}D\right]$ is calculated for the structure with doped passive graphene sheets at zero temperature, and presented in Fig. 2a.

For illustration purposes, we consider an infinite structure with the period of $d_{1}=300\,\mathrm{nm}$ and the permittivity of the dielectric $\varepsilon_{1}=1$. In Fig. 2b we present $\text{K}D$ values (isofrequencies), calculated for three excitation energies $\tilde{\omega}$, marked as dashed lines in Fig. 2a. For $\tilde{\omega}=0.4$, as seen from the bottom panel of Fig. 2b, $\operatorname{Re}\left[\text{K}D\right]=\left[0,\pi\right]$ and $\operatorname{Im}\left[\text{K}D\right]=0$ over the range $k_{x}/k_{0}\in\left[1,\,56\right]$, being characteristic of the hyperbolic operation of the structure (the light cone in this normalized representation is a vertical line at $k_{x}/k_{0}=1$). For these excitation parameters the graphene sheets of the structure are operating in the loss-less region II of the Fig. 1, and thus do not exhibit any inter- or intra-band losses. As a result, the plasmonic modes are loss-less in region $k_{x}/k_{0}\in\left[1,\,56\right]$, as seen from $\operatorname{Im}\left[\text{K}D\right]=0$. For $k_{x}/k_{0}>56$, the structure exhibits loss, as seen from the $\operatorname{Im}\left[\text{K}D\right]<0$. Since excitations occur within the loss-less region II, this damping must originate from the periodic structure itself. The condition $\operatorname{Re}\left[\text{K}D\right]=\left\lvert\operatorname{Im}\left[\text{K}D\right]\right\rvert=\pi$ for excitation at the band edge is marked as a dashed black line in Fig. 2a.

In the middle panel of Fig. 2b the loss-less hyperbolic modes are located within $k_{x}/k_{0}\in\left[45,\,97\right]$, exhibiting band gaps for $k_{x}/k_{0}<45$ and $k_{x}/k_{0}>97$, seen as $\operatorname{Im}\left[\text{K}D\right]<0\cup\operatorname{Re}\left[\text{K}D\right]=0$ and $\left\lvert\operatorname{Im}\left[\text{K}D\right]\right\rvert>\pi\cup\operatorname{Re}\left[\text{K}D\right]=\pi$ respectively. The lower boundary of the hyperbolic modes region can be expressed as $\operatorname{Re}\left[\text{K}D\right]=\operatorname{Im}\left[\text{K}D\right]=0$, and is shown as solid black line in Fig. 2a.

In the top panel of the Fig. 2b we can see absence of a hyperbolic mode, and increasingly larger losses, present both in regions II and III (boundary between them shown as a thin gray line in Fig. 2a).

III.2 Homogenization

The role of thin metallic layers in a general hyperbolic metamaterial is to provide a confining medium for plane plasma oscillations, excited by TM polarized excitations (see Fig. 3, here, we are concerned with TM-polarized waves, as they couple to strongly-bound TM plasmon modes of graphene (Page et al., 2015)). The $x-$component of incident wavevectors will couple to the plasma oscillations and at resonant frequencies provide the necessary negative $\varepsilon_{x-eff}$ component of the dielectric tensor of the effective medium, giving it overall hyperbolic properties.

In case of a metal with finite thickness, both $\varepsilon_{x-eff}$ and $\varepsilon_{z-eff}$ are affected by the negative $\varepsilon$ of the metal, since due to the metal’s finite thickness both in- and out- of plane electric field components couple to the plasma oscillations. A detailed procedure of finding the components of the dielectric tensor of a general layered metamaterial has been described, for example, in (Agranovich and Kravtsov, 1985; Kidwai et al., 2012; Wood et al., 2006).

On the contrary, in the case of the metallic layer thickness shrinking to zero, the plasma is confined to two dimensions (Fig. 3b), and so cannot be polarized in $z-$direction. The $z-$component of the effective permittivity of the structure is thus described exclusively by the dielectric response: $\varepsilon_{z-eff}=\varepsilon_{d}$. The energy conservation for the effective structure will be expressed as:

[TABLE]

Here, $\text{K}=\text{K}^{{}^{\prime}}+i\text{K}^{{}^{\prime\prime}}$ is the wavevector of the Bloch modes of the structure, calculated from Eq. (3). From this, the effective $\varepsilon_{x-eff}$ can be obtained:

[TABLE]

where the quantity $k_{0}^{2}-\frac{k_{x}^{2}}{\varepsilon_{d}}$ is negative outside of the light cone. Analysis of Eq. (5) provides us with a tool-box for an analytic calculation of the regions of hyperbolic behavior and the plasmonic band gap of the structure. Concretely, we can find the regions in the $\omega/k_{x}$ phase space where the following conditions are satisfied:

[TABLE]

From Eq. (5) we can see that $\left[\varepsilon_{x-eff}\right]<0$ for $0\leq\text{K}^{{}^{\prime\prime}2}\leq\text{K}^{{}^{\prime}2}\leq\pi/d_{1}$, with two edge conditions (a) and (c) bounding the phase space of hyperbolic character. The modes are undamped between conditions (a) and (b), and experience spatial damping between conditions (b) and (c). We can derive a corresponding set of analytic expressions of the boundaries in the $\omega/k_{x}$ space, on which the conditions (6) are satisfied. In terms of the imaginary part of graphene’s non-local conductivity $\sigma_{s}^{{}^{\prime\prime}}$ for the loss-less case these boundaries have the following form(details of the derivation may be found in the appendix):

[TABLE]

These conditions are satisfied along the solid, dotted and dashed black lines respectively, as presented in Figs. 2a and 2c. The black lines in Fig. 2c converge to the single sheet plasmon dispersion for large values of the excitation frequency and its in-plane component. Increasing thickness of the dielectric layer weakens the coupling of the plasmonic modes between the graphene sheets, and thus the region of $\operatorname{Re}\left[\varepsilon_{x-eff}\right]<0$ would shrink. These and other properties can be inferred from an analysis of Eqs. (7), which are central to the study of the hyperbolic behavior of the considered plasmonic metamaterials.

III.3 Mode amplification in active HMMs

Plasmonic metamaterials are inherently lossy due to their exploitation of the plasma resonances. This causes the wavevector of the plasmonic modes to have a large imaginary part, corresponding to a decay of the modes in space. In order to overcome this problem, it was proposed to incorporate gain medium in the dielectric layers (Soukoulis and Wegener, 2010; Xiao et al., 2010; Savelev et al., 2013). Here, we propose to explore gain, available from graphene by optically pumping the carrier plasma into a state of inversion. Ideally, the plasmons, coupling to the inverted carriers, would cause their stimulated recombination and undergo amplification. Here we employ the exact in RPA quantum conductivity model in order to describe the response of the inverted carrier plasma in the doped graphene sheets.

In the case of a single graphene sheet, the conditions for excitation of plasmons are satisfied on the line in the $\left(q,\,\omega\right)$ space (Figs. 1a and 1b), while in the case of the multiple coupled graphene sheets, the whole region of plasmonic modes opens up in the $\left(q,\,\omega\right)$ space, as seen from the $\operatorname{Re}\left[\varepsilon_{x-eff}\right]<0$ region in Fig. 2a.

A non-zero imaginary part of the wavevector defines the change in the amplitude of the Bloch wave, as it propagates through the structure. This parameter acquires meaning only in conjunction with $\operatorname{Re}\left[\text{K}D\right]$, such that $\operatorname{Im}\left[\text{K}D\right]<0\,\cup\,\left|\operatorname{Im}\left[\text{K}D\right]\right|<\left|\operatorname{Re}\left[\text{K}D\right]\right|$ corresponds to damped modes, $0<\operatorname{Im}\left[\text{K}D\right]<\operatorname{Re}\left[\text{K}D\right]$ - to amplified modes and $\left|\operatorname{Im}\left[\text{K}D\right]\right|>\operatorname{Re}\left[\text{K}D\right]$ - to unstable modes. Therefore, generally, stable modes have to satisfy

[TABLE]

To illustrate this, we are using the quantum conductivity model together with Eq. (3) to perform a calculation of $\operatorname{Im}\left[\text{K}D\right]$ for an active structure, which is plotted in Fig. 4b. The regions of gain (purple) and loss (blue) of the considered GHMM coincide with the single graphene sheet regions of the inter-band recombination and excitation respectively. The condition $\operatorname{Re}\left[\varepsilon_{x-eff}\right]<0$ is satisfied over the region which is bound by the black curve in the Fig. 4b. As can be seen, its location is independent of the zones of the inter/intra-band processes in the single-sheet graphene, and spans over both amplifying and lossy regions.

The region where $\operatorname{Re}\left[\varepsilon_{x-eff}\right]<0$ in Fig. 4b is shown for a particular thickness of the dielectric layer. Varying the dielectric thickness changes the coupling strength of the plasmonic modes (Das Sarma and Quinn, 1982), thus changing the phase space of the plasmonic region, with the top right corner following the single sheet plasmon dispersion, and approaching the origin with an increase of the dielectric thickness.

The isofrequency contours of representative modes are shown in Fig. 4a for three different values of the excitation frequency $\omega_{0}\hbar/\mu$ as real and imaginary parts of the isofrequency curves. For the frequency $\omega_{0}\hbar/\mu=0.4$ the hyperbolic mode is seen as a continuous curve $\operatorname{Re}\left[\text{K}D\right]\in\left[1,\,\pi\right)$ for the in-plane wavevectors of $k_{x}/k_{0}\in\left[1,\,200\right]$. It is amplified in this range, but amplification is stable only for $k_{x}/k_{0}\in\left[1,\,230\right]$, - within the $\epsilon_{x-eff}\leq 0$ loop in Fig. 4b. Analogous behavior is observed for excitation frequencies in the range $\omega_{0}\hbar/\mu\in\left(0,\,0.5\right)$, where the in-plane components are fully within region I.

For $\omega_{0}\hbar/\mu=0.8$, the modes for $k_{x}/k_{0}\in\left[1,\,110\right]$ are located within the loss-less region II, where $\operatorname{Im}\left[\text{K}D\right]\approx 0$ and $\operatorname{Re}\left[\text{K}D\right]>0$. For $k_{x}/k_{0}\in\left[110,\,240\right]$ the modes enter the region of gain I, thus being amplified, with $\operatorname{Im}\left[\text{K}D\right]>0$.

The last set of modes is demonstratively represented at $\omega_{0}\hbar/\mu=1.2$. There is a band-gap within region II for $k_{x}/k_{0}\in\left[1,\,80\right]$, where $\operatorname{Re}\left[\text{K}D\right]\approx 0$ and $\operatorname{Im}\left[\text{K}D\right]<0$. The over-damped mode exists for $k_{x}/k_{0}>80$ due to the inter-band loss of region III, preventing perfect destructive interference, which would be causing the band-gap.

IV Pulse propagation in a finite stack

IV.1 Loss Compensation

So far we have considered light propagating in infinite periodic stacks. Despite being a good approximation for a description of sufficiently large structures, it is not clear how this model would generalize to structures at the current level manufacturing possibilities, consisting of just several layers. To answer this question we use the transfer matrix method (TMM), which allows us to spatially resolve the electric and magnetic fields, as well as calculate the Poynting vector, inside of the finite structure with the given set of input or output fields. With this we calculate the quasi-static field distribution, resulting from the illumination of the finite stack with a narrow spatial Gaussian intensity envelope. The Gaussian envelope is sufficiently small in real space, thus consisting of $k_{x}$ components that are evanescent in free space, and coupling to the plasmonic modes of the considered heterostructure.

We start with calculations of the forward-traveling component of the electric field distribution in a structure consisting of 14 unit cells, as shown in the Fig. 5a. In this initial scenario the graphene sheets are assumed to be at zero temperature, have a small amount of scattering loss to ensure numerical stability of the calculations, and doped at a level higher than the energy of the exciting filled. Thus the structure is operating in the loss-less region II in Fig. 1a. Therefore, the propagating fields do not undergo any damping or loss, though the overall intensity is decreasing due to the fields being reflected from the subsequent graphene-dielectric boundaries (this can also be seen as increase of intensity of the backward-propagating fields).

Subsequently, we keep the parameters of the structure and stimulation the same, but introduce a finite temperature (resulting in a smearing of the graphene Fermi level, while graphene is still doped to the same level), and add collision loss in the graphene sheets. The calculated field intensity distribution is shown in Fig. 5b. As a result, the depth of propagation of the pulse is greatly reduced as compared to the loss-less case. The propagating pulse also undergoes dispersion, as can be seen from the increased width of the pulse. We notice that the presence of loss introduces an energy flow within the structure, manifested in $S_{z}>0$. This results in degraded imaging performance of the stack.

Finally, we introduce gain to our heterostructure in the form of optical inversion of graphene, entering the model as carrier imbalance of the electrons and holes. The temperature conditions and scattering losses are yet again the same as in the previous calculation. The resultant field distribution is presented in Fig. 5c. It can be seen that the length of propagation of the pulse within the structure has increased as compared with the scenario in the Fig. 5b, as well as showing decreased dispersion (the pulse branches width is decreased as compared to the Fig. 5b). The energy lost due to the temperature effects and scattering is now well compensated. All this demonstrates that with added gain in the form of optical inversion of graphene we can greatly improve pulse propagating properties of the finite stack, even ensuring stable operation with $S_{z}\rightarrow 0$.

IV.2 Using results for infinite structures for the analysis of finite structures

So far we have computationally investigated the propagation of a small Gaussian envelope in a finite stack, as well as studied the behavior of an infinite stack, deriving a set of analytic tools for analysis and prediction of the properties of an infinite structure. An interesting question to answer is would it be possible to use analytic approach used for an infinite structure in order to predict behavior of a finite stack, without having to perform computationally intense simulations of the finite stack?

Let us consider the TMM framework (which provides the basis for the analysis of a finite stack). To propagate a narrow spatial envelope through the stack we decompose it into a finite number of Gaussian spatial components $k^{\prime}_{x}$, propagating them individually through the stack and recombining them at the output end with their Fourier weights. In doing so we impose a condition on the fields at the output layer to be traveling out of the structure only. This means that the obtained field distributions for a particular $k^{\prime}_{x}$ in general are not the eigenmodes of the structure, but are a linear combination of two counter-propagating eigenmodes. The eigenmodes of the finite structure are calculated as eigenmodes of its single unit cell (that is also what we are calculating in Eq. (3)).

Since the eigenmodes of the unit cell of infinite and finite structures are calculated identically, we can use the isofrequency contours of an infinite structure to predict the behavior of a finite one. Therefore, in order to ensure a loss-less propagation of the pulse through the finite structure, the Fourier components of the pulse should be located within the range of $k_{x}/k_{0}$ where $\operatorname{Im}\left[\text{K}D\right]=0$. At the same time, regimes where $\operatorname{Im}\left[\text{K}D\right]>0$, although showing gain, are also detrimental to the pulse propagation, resulting in non-zero Poynting vector. In this regard, $k_{x}/k_{0}$ can be treated as a resolution limit of a finite structure, as it defines the minimum size of the pulse which can propagate, and can be obtained from the isofrequency contours for an infinite structure for the condition where $\operatorname{Im}\left[\text{K}D\right]=0$.

As was shown earlier, the condition of $\operatorname{Im}\left[\text{K}D\right]=0$ is satisfied when the stack is supporting plasmonic modes in the loss-less case, or in case of losses compensated by gain. Therefore, the resolution limit of the structure is defined as the boundary of plasmonic band of operation $\varepsilon_{x-eff}=0$, shown as black contours in Figs. 2 and 4. The analytic calculation of this condition is described in the Section III and in the Appendix.

V Conclusion

We have studied the behavior of graphene-based hyperbolic metamaterials. Employing an RPA model for the graphene conductivity, which accounts for the plasmonic gain and loss, we developed an analytic tool-set for the characterization of the regions of behavior of the infinite structures. This allows us to uncover the possibility of plasmonic amplification in optically excited heterostructures. Considering structures with the finite number of layers, we demonstrate that the theory assuming infinite structures has a good applicability there as well. We also demonstrated the quasi-equilibrium field distributions of the Gaussian envelopes propagating in graphene-based hyperbolic metamaterials, and showed that these structures can be used to compensate occurring graphene losses.

VI Appendix

VI.1 Quantum polarizability model

In general, $\Pi(q,\omega)$ in the Eq. (1) is the irreducible polarizability, quantifying the response of the particle/hole plasma to the external excitations. In RPA, it is approximated by its first order term, describing the non-interacting particle hole plasma, and is expressed via the Lindhard Formula:

[TABLE]

This expression sums over all possible intra-band $\left(\mathbf{k}\right)\rightarrow\left(\mathbf{k}+\mathbf{q}\right)$ and inter-band $s\rightarrow s^{{}^{\prime}}$ transitions, contains spin/valley degeneracy within the constant $g=4$, and is weighted with the square of the transition matrix $M_{\mathbf{k},\mathbf{k}^{\prime}}^{ss^{\prime}}=[1+ss^{\prime}\cos(\theta_{k,k^{\prime}})]/2$.

From Eq. (9) an analytic expression for polarizability of gap-less graphene at zero temperature has been derived in (Page et al., 2015; Pyatkovskiy, 2009), and with some rescaling can be expressed as (Page, 2016):

[TABLE]

where $\tilde{\Pi}\left(\tilde{q},\tilde{\omega}\right)=-4+\tilde{q}^{2}\frac{G^{+}\left(\frac{2+\tilde{\omega}}{\tilde{q}}\right)+G^{-}\left(\frac{2-\tilde{\omega}}{\tilde{q}}\right)}{2\sqrt{\tilde{q}^{2}-\tilde{\omega}^{2}}}$ and $G^{\pm}\left(z\right)=z\sqrt{1-z^{2}}\pm i\operatorname{arccosh}\left(z\right)$. This expression describes the polarizability of doped graphene sheet at equilibrium, with the corresponding chemical potential $\mu$. The particle-hole symmetry is implicit in $\Pi\left(q,\omega\right)|_{\mu}^{T=0}=\Pi\left(q,\omega\right)|_{-\mu}^{T=0}$. The expression (10) is analytic for $\operatorname{Im}\left[\omega\right]>0$, and can be evaluated for $\operatorname{Im}\left[\omega\right]<0$ via the analytic continuation (Page et al., 2015).

The single sheet plasmon dispersion is calculated exactly for the doped equilibrium graphene and is shown in the Figs. 1a and 6a as a green solid line. It begins in the region II, where the excitation energies are below the chemical potential of doping, thus are insufficient to decay into the inter-band transitions ($\gamma_{pl}=0$). At zero doping region II disappears, and all phase space within the Dirac cone is occupied by the region III. After region II the plasmon dispersion passes through the region III with Landau damping; region IV with intra-band excitations and loss-free region V.

Optical pumping of graphene would create a bath of excited electron-hole pairs. This opens up a region I in the phase space of regions II and III. The plasmons within region I would trigger recombination of electron-hole pairs, thus undergoing coherent amplification. For comparison, the single-sheet plasmon dispersion of the extrinsic active graphene is shown in the Fig. 6b.

VI.2 Modes of passive graphene-based HMM

The sheet conductivity $\sigma_{s}(q,\omega)$ in the Eq. (3) is a complex quantity, and can be expressed as $\sigma_{s}(q,\omega)=\sigma_{s}^{{}^{\prime}}+i\sigma_{s}^{{}^{\prime\prime}}$. Then, (3) becomes:

[TABLE]

were a=-i\mbox{\sqrt{-\left(k_{x}/k_{0}\right)^{2}+\varepsilon_{d}}}. The Eq. (10) yields $\sigma_{s}^{{}^{\prime}}=0$ in the region II in the Figs. 1 and 6a, as calculated for the passive graphene at $T=0$. Thus, Eq. (11) becomes:

[TABLE]

Now it is possible to establish the analytic solution for boundaries of stable hyperbolic behavior of the structure, with $\varepsilon_{x-eff}<0$. The hyperbolic modes occupy the phase space where $\operatorname{Re}\left[\text{K}D\right]\in\left(0,\,\pi\right)$, and have to fulfill the condition in Eq. (8) in order to be stable. Combining these, the two boundaries of the modes are:

[TABLE]

Thus, using these rules, the conditions for $\sigma_{s}^{{}^{\prime\prime}}(q,\omega)$ can be derived from the Eq. (12):

[TABLE]

Finding the values of $(q,\omega)$, at which the above equalities are satisfied, will define the region $\varepsilon_{x-eff}<0$.

The same procedure can be used to find the boundary of undamped hyperbolic modes, where $\operatorname{Re}\left[\text{K}D\right]=\pi$, and $\operatorname{Im}\left[\text{K}D\right]=0$:

[TABLE]

The conditions $(a)$ and $(b)$ of the Eqs. (14), bounding the region of $\varepsilon_{x-eff}<0$, are plotted as black solid and black dashed lines respectively in the Figs. 2a and 2c. Within that region, the undamped modes have their $k_{x}/k_{0}$ values bound by the dotted curve, given by the Eq. (15). All three curves tend to the single graphene sheet dispersion curve, shown as a green line. As the distance between the sheets increases ($d_{1}\rightarrow\infty$), the region of $\varepsilon_{x-eff}<0$ becomes smaller, converging to $\sigma_{s}^{{}^{\prime\prime}}(q,\omega)=\frac{2\varepsilon_{1}}{Z_{0}a}$. Expectedly, this coincides with the single-sheet plasmon dispersion.