Gate tunable light–matter interaction in natural biaxial hyperbolic van der Waals heterostructures

2.1 Biaxial hyperbolic phonon polaritons in α-MoO₃

The α-MoO₃ unit cell has lattice constants of around 0.396 nm, 0.369 nm and 1.385 nm in the x, y and z crystallographic directions [4] as shown in Figure 1(a). The orthorhombic nature of the unit cell results in different sets of longitudinal optical (LO) and transverse optical (TO) phonons in the three crystallographic directions. The three sets of LO and TO phonons constitute three RBs in the mid-infrared spectral region, where α-MoO₃ exhibits natural biaxial hyperbolicity. Spectral region for RB-1, RB-2 and RB-3 of α-MoO₃ lie in the ranges of 545 cm⁻¹–850 cm⁻¹, 820 cm⁻¹–972 cm⁻¹, and 958 cm⁻¹–1006 cm⁻¹, respectively, where real part of dielectric permittivity is negative along y, x and z crystallographic directions, respectively, as shown in Figure 1(b) [46]. Since α-MoO₃ is a biaxial hyperbolic material, the dispersion of PhPs strongly depends on angle between plane of propagation of the polaritons and x-crystallographic direction as shown in Figure 1(a) and represented by ϕ. We investigate the dispersion of PhPs in an air/α-MoO₃/Si heterostructure at ϕ = 0°, 45° and 90° as shown in Figure 2(a)–(c). Here, the thickness of α-MoO₃ is chosen to be 50 nm, and the imaginary part of the Fresnel reflection coefficient for p-polarized light is plotted in frequency-momentum space.

Figure 2:

Dispersion of polaritonic modes in air/graphene/α-MoO₃ (50 nm)/Si heterostructure as a function of ϕ and chemical potential (μ) of graphene. Dotted lines indicate LO phonon frequencies, and solid lines indicate TO phonon frequencies along y-(green), x-(orange) and z-(blue) crystallographic directions of α-MoO₃. The color axis represents imaginary part of unitless quantity Fresnel’s reflection coefficient for p-polarized light.

Figure 2(a) reveals discrete modes in the RB-2 and RB-3 spectral regions for ϕ = 0°, which are attributed to the hyperbolic anisotropy of α-MoO₃ resulting in the dispersion of fundamental and higher-order modes of hyperbolic PhPs. These modes demonstrate positive dispersion (i.e., v _g  > 0) in RB-2 spectral region (Type-2 hyperbolicity) and negative dispersion (i.e., v _g < 0) in RB-3 spectral region (Type-1 hyperbolicity). Negative dispersion in the RB-3 spectral region is associated with the negative real part of dielectric permittivity in the out-of-plane crystallographic direction. We also observe similar PhP modes for ϕ = 90° in the RB-1 and RB-3 regions (Figure 2(c)). For ϕ = 45° we observe Type-1 hyperbolic PhP modes in both RB1 and RB2 as the in-plane (ϵ ₁₁, ϵ ₁₂, ϵ ₂₁, ϵ ₂₂) entries of the rotated permittivity tensor depend on the permittivity along both x and y principal directions. However, these modes do not span the entire frequency range of either RB. Equation (1) suggests that the ϵ ₁₁ component of the rotated tensor is the most relevant in determining the extent of these PhP modes.

2.2 Hybridization of biaxial hyperbolic PhPs and SPPs

We examine the hybridization of SPPs of graphene with hyperbolic PhPs of α-MoO₃ via modulation of chemical potential (μ) in graphene. We consider a heterostructure having the configuration: air/graphene/α-MoO₃/silicon as shown in Figure 1(a) and explore the dispersion of polaritons, shown in Figure 2. The conductivity of graphene for various values of chemical potential is shown in Figure 1(c). For μ = 0.025 and 0.05 eV, we observe smearing of hyperbolic PhP modes in all the RBs of α-MoO₃ (Figure 2(d)–(i)), however, there are no significant changes in the dispersion of hyperbolic PhPs in the plotted range of wave-vectors. This is due to interband transition in graphene upon excitation by a photon with an energy greater than 2μ (around 403 cm⁻¹ and 806 cm⁻¹ for μ = 0.025 and 0.05 eV, respectively). At μ = 0.1 eV in Figure 2(j), a broad mode is observed above the RB-3 spectral region and below the RB-2 spectral region. This mode corresponds to the SPPs of graphene, since as the chemical potential increases, the graphene plasmon dispersion moves closer to the light cone [28]. Due to the inherently lossy nature of the SPP mode, it is broader than the phonon modes. Furthermore, in Figure 2(k) and (l), a variation is observed in the dispersion of the fundamental mode of hyperbolic PhPs in the RB-1 spectral region. This is ascribed to the hybridization of hyperbolic PhPs of α-MoO₃ with SPPs of graphene forming HPPPs. Moreover, for μ = 0.3 and 0.5 eV, the dispersion of hyperbolic PhPs in RB-1 and RB-2 spectral region become steeper compared to the no graphene case, which is a strong indication of hybridization of hyperbolic PhPs of α-MoO₃ with SPPs of graphene (Figure 2(m)–(r)). Additionally, we observe that the dispersion curve of SPPs becomes nearly horizontal as it approaches the TO phonon frequencies of α-MoO₃. This behaviour is due to the avoided-crossing [28, 30], which suggests a coupling between surface plasmons of graphene and TO phonons of α-MoO₃ and further substantiates the formation of HPPPs. Dispersion curves in Figure 2(o) and (r) demonstrate that these hybrid modes can propagate outside the RB-1 and RB-2 spectral region, respectively, providing additional operational bandwidth compared to HPhPs of α-MoO₃ and relatively less losses compared to the SPPs of graphene.

Apart from fundamental modes of hyperbolic PhPs in RB-1 and RB-2, variation in the dispersion of higher-order PhP modes is observed upon integration of graphene. For a particular mode index (k _ρ/k ₀), a blue-shift is observed in the frequencies of higher-order modes with increasing chemical potential of graphene, which we will discuss in the sections below. Similar to the observed trend in RB-1 and RB-2 spectral regions, hyperbolic PhP modes in the RB-3 spectral region are strongly dominated by graphene’s interband transitions at lower μ values (0.025 and 0.05 eV) and graphene’s SPPs at higher μ values. At μ = 0.3 and 0.5 eV, mode repulsion is observed near the out-of-plane TO phonon frequency of α-MoO₃ which indicates the coupling of SPPs with these TO phonons. A discussion of the effect of thickness of α-MoO₃ on dispersion of HPPPs is presented in the Section S2 of supplementary information.

Hybridization of SPPs and hyperbolic PhPs is further investigated via iso-frequency surfaces for these polaritons, plotted as a function of ϕ, μ and k _ρ/k ₀. We consider ω = 700 cm⁻¹, 900 cm⁻¹, 960 cm⁻¹ and 980 cm⁻¹, which are representative of frequencies in the RB-1, RB-2, transition region between RB-2 and RB-3, and RB-3 spectral regions in Figure 3. At ω = 700 cm⁻¹ in Figure 3(a), we observe that the fundamental and higher-order hyperbolic PhP modes of α-MoO₃ appear as hyperbolae aligned in the ϕ = 90° direction (i.e. y-direction). The iso-frequency surfaces of the heterostructure of α-MoO₃ and graphene with μ = 0.1, 0.3, 0.5 eV appear as ellipses squeezed along ϕ = 90° due to the hybridization of hyperbolic PhPs of α-MoO₃ with the isotropic SPPs of graphene along y-direction. Furthermore, we observe a shift towards lower mode index in the first-order mode of hyperbolic PhPs with increasing chemical potential of graphene compared to that of no graphene (see Figure S5 of supplementary information). Figure 3a confirms that the higher order PhP modes are affected by graphene integration with α-MoO₃. For μ = 0.3 and 0.5 eV, the SPP wave vector decreases resulting in smaller ellipses compared to the μ = 0.1 eV case. Similarly, at ω = 900 cm⁻¹ in Figure 3(b), the hyperbola is found aligned along ϕ = 0° (x-direction) and first-order mode inside hyperbola varies with increasing chemical potential. These results confirm the presence of tunable hybridization of fundamental and higher-order PhP modes with SPPs via chemical potential of graphene. Polar plots for other chemical potentials at ω = 700 cm⁻¹ and 900 cm⁻¹ are presented in the supplementary information. At ω = 960 cm⁻¹, we do not observe hyperbolic PhP modes for the no graphene case in Figure 3(c). This is associated with the inherent dielectric response of α-MoO₃ in the transition from RB-2 to RB-3 spectral region. However, graphene SPPs dominate the formation of HPPPs for μ = 0.3 and 0.5 eV as shown by the squeezed ellipse in Figures 3(c) and S6 of supplementary information. Moreover, at ω = 980 cm⁻¹, we observe the transition of hyperbolic contour into an elliptical contour for the no graphene case as shown in Figure 3(d). For μ = 0.1 eV, we observe that the elliptical contour is elongated along x-crystalline direction, which is attributed to mode repulsion between HPPP modes in the RB-2 and RB-3 spectral regions for ϕ = 0°. As a result, we observe higher mode index at 980 cm⁻¹ along x-direction. On the other hand, HPPP modes in the RB-2 spectral region do not exist for ϕ = 90° and consequently, we have a lower mode index at 980 cm⁻¹ along y-direction as compared to the x-crystalline direction (see Figure S6 of supplementary information). Hence the elliptical contour for μ = 0.1 eV is elongated along x-crystal direction. Lastly, at μ = 0.3 and 0.5 eV, we observe a feeble elliptical contour in Figure S4 of supplementary information (which are absent in Figure 3(d) is due to relatively greater strength of contour for no graphene case). This can be attributed to type-1 HPPPs in the RB-3 spectral region. After integration of graphene with μ = 0.3 and 0.5 eV, the HPPP modes in the RB-3 spectral region move upward (or towards higher mode index) and create a spectral gap around 980 cm⁻¹ (see Figures S4 and S6 of supplementary information), resulting in a feeble elliptical contour.

Figure 3:

Polar plots of HPPP modes in 50 nm thin α-MoO₃ and graphene heterostructure as a function of graphene’s chemical potential (0.1, 0.3, 0.5 eV): HPPP fundamental and first-order modes are represented by imaginary part of Fresnel’s reflection coefficient for p-polarized light. Radial and azimuthal directions correspond to mode index (k _ρ/k ₀) and ϕ, respectively. Effect of graphene’s chemical potential on HPPP modes is demonstrated by green, orange and blue color for μ = 0.1, 0.3, and 0.5 eV, respectively. The fringes at the edges of the polar plots are the artefact of the scatter type plotting scheme.

In order to understand the origin of this active tunability and the angle dependence of these hybrid modes, we derive an analytical expression for these dispersion relations under the quasistatic approximation (k _ρ ≫ k ₀). The expression for the dispersion curve then is given by:

Here, σ is the conductivity of graphene, Z 0 = μ 0 / ϵ 0 is the impedance of free space, ϵ _z is the z-direction (out of plane) permittivity of α-MoO₃ and ψ = i ϵ z / ϵ x cos 2 ϕ + ϵ y sin 2 ϕ and n is an integer. ϵ ₁ = 1 and ϵ ₃ = 2.5 are taken as the permittivities of the air and Si substrate, respectively. This expression can be derived by imposing the condition that the total of phase shift due to reflection at the two interfaces of the α-MoO₃ slab plus the propagation phase shift should be an integral multiple of 2π (see supplementary information). The variation of the graphene Fermi level through a gate voltage modifies σ which affects the reflection phase shift at that interface, thereby enabling the tunability of the hybrid modes. It must be noted that all the results presented in this work are full wave calculations without using the quasi-static approximation.

2.3 Anisotropic spontaneous emission via biaxial hybrid polaritons

Next, we explore the anisotropic SER enhancement from an electric dipole placed in the vicinity of graphene/α-MoO₃/Si heterostructure, as shown in Figure 1. To the best of our knowledge, anisotropic SERs for natural biaxial hyperbolic material have not been reported in the literature so far. Since α-MoO₃ exhibits in-plane hyperbolicity, we examine the anisotropic SERs by considering the polarization of the electric dipole along x-, y-, and z-directions. Distance of the electric dipole from the surface of the heterostructure (represented by z ₀) is taken to be 50 and 200 nm, whereas μ is varied as 0.1 eV, 0.3 eV and 0.5 eV, as shown in Figure 4. For comparison, we have also evaluated the anisotropic SERs through α-MoO₃ without graphene and for lower μ values, shown in the Section S4 of supplementary information. Integration of graphene with α-MoO₃ crystals not only provides tunable spontaneous emission but also enhances the SER to a much broader operational bandwidth.

Figure 4:

Anisotropic spontaneous emission rates for dipole placed above air/graphene/α-MoO₃ (50 nm)/Si heterostructure: (a)–(c) represent SER enhancements at different μ values for x-, y-and z-polarized electric dipole, respectively, at z ₀ = 50 nm. (d)–(f) represent SER enhancements at different μ values for x-, y- and z-polarized electric dipole, respectively, at z ₀ = 200 nm. The green, red and blue shaded regions represent the RB-1, RB-2, and RB-3 spectral regions of α-MoO₃, respectively.

At z ₀ = 50 nm, the order of magnitude of SER is found to be more than 10³ inside the RB spectral regions for x, y and z oriented dipoles as shown in Figure 4(a)–(c) respectively. Anisotropic SER spectra are dominated by graphene SPP modes outside the RBs of α-MoO₃. The difference between SER spectra in the RB spectral regions with and without integration of graphene and α-MoO₃ is associated with the formation of HPPPs, which provides a larger operational bandwidth. Variation of SER spectra in the RBs as a function of chemical potential is due to the change in the dispersion of HPPP modes with chemical potential, as discussed above. Spectral dips in SERs spectra are attributed to the spectral gap due to the absence of strongly confined SPP modes and hyperbolic PhP modes in the avoided crossing regions of the dispersion. Within such an avoided crossing, the dipole cannot couple to any mode effectively and hence spontaneous emission is suppressed. Here, we see that the line-widths of spectral dips decrease with increasing chemical potential, which is due to the fact that the graphene plasmon mode index decreases upon increase of chemical potential [28], decreasing the strength of the avoided crossing at larger mode indices just below the RBs of α-MoO₃. Moreover, a blue-shift in the spectral dip is observed from 960 cm⁻¹–980 cm⁻¹ when μ increases from 0.1 to 0.3 eV. As described above with Figure 3(c) and (d), this is due to the shift of the RB-3 HPPPs towards higher mode index as the graphene Fermi level increases – this trend is opposite to that of the HPPPs in the other bands.

At z ₀ = 200 nm, (shown in Figure 4(d)–(f)), we observe a decrease in the SER by one order of magnitude due to the evanescent nature of SPPs and HPPPs, resulting in reduced interaction of polaritonic modes with the electric dipole and hence reduction of the SERs. The dips in the SER spectra are due to the spectral gap of polariton modes, as discussed above. Reduction in the line-width of spectral dips with increasing chemical potential is found as expected. Below RB-1 spectral region, the SER for μ = 0.3 and 0.5 eV is higher in comparison with smaller μ values. Since higher chemical potential increases the imaginary part of the graphene conductivity, it reduces the confinement of SPPs. This can be seen from dispersion of SPPs in Figure 2, where one can observe SPP modes at higher μ values just below the RB-1 spectral region move to smaller k _ρ values. As a result, SPP modes for higher chemical potential have an evanescent tail at a larger distance away from the surface compared to the SPP modes for smaller chemical potential and hence the SER can be enhanced. A mathematical justification can further be provided based on the envelope k ρ 2 e − 2 k ρ z 0 which peaks at k _ρ ∼ 1/z ₀. So depending on whether the mode is resonant at k _ρ ∼ 1/z ₀, one can observe an increase in the SER. The variation of the anisotropic SER spectra as a function of the thickness of α-MoO₃ is discussed in the Section S4 of the supplementary information. Similarly, spontaneous emission spectra for the air/graphene/Si heterostructure are discussed in Section S5 of the supplementary information. These SER spectra substantiate our claim that one can observe a remarkable spontaneous emission enhancement for an electric dipole in the vicinity of α-MoO₃ integrated with graphene over a broad spectral range.

2.4 Application to valley quantum interference

We demonstrate the power of such a tunable anisotropic SER for application in spontaneous valley coherence [16]. Quantum coherence in solid state systems is important not only from a fundamental point of view but is a key ingredient for realizing the promised applications in quantum technologies like quantum computing [47], quantum key distribution [48] and quantum metrology [49]. In the infrared range, gapped Dirac systems [50, 51] such as biased bilayer graphene [52, 53] and similar systems could serve as an important excitonic building block for such technologies. The electronic bandstructure of these systems consists of two degenerate valleys K and K′, in which the optical selection rule is sensitive to the polarization (σ ^±) of the photon. In order to harness this “valley degree of freedom”, it is important to be able to actively control the coherence between excitons in the two valleys. While this coherence has been demonstrated in the stimulated regime [54], where the valley excitons inherit coherence from the exciting laser, a spontaneous route to such coherence generation would greatly simplify the device architecture. Such a spontaneous route to valley coherence has been previously demonstrated using structured metamaterials [44]. Our graphene/MoO₃ heterostructure would allow for easy fabrication and enable electrostatic control of valley coherence, which could be easily implemented on a chip and provide a low-cost, facile fabrication and chip-scale alternatives for building relevant technologies.

Valley coherence can be quantified experimentally using the degree of linear polarization or DoLP of the emission, where DoLP = (I _X  − I _Y )/(I _X + I _Y ), I _X and I _Y being the intensities of the two linearly polarized emissions. It can easily be shown that in the presence of a weak pump, the DoLP is approximately equal to the quantum interference Q, which is defined as [16, 44]:

where γ _x and γ _y are the Purcell factors for x and y polarized dipoles, respectively. The valley quantum interference is proportional to the anisotropy between the Purcell factors for orthogonal dipoles. The geometrical configuration for this experiment would involve a gapped Dirac system such as a bilayer graphene system placed a certain distance away from the monolayer graphene/α-MoO₃ system. The exciton of this gapped Dirac system will be represented by the dipole. The quantum interference for this heterostructure with a 50 nm thick α-MoO₃ layer is plotted as a function of frequency for various values of chemical potential μ of graphene in Figure 5. Since α-MoO₃ is a naturally biaxial material, we observe large magnitude of quantum interference in RB-1 and RB-2 spectral regions, where dipole is aligned along the y- and x-crystallographic axes, respectively, and can couple more strongly to PhP modes as compared to the dipole polarized orthogonal to the metallic direction of the hyperbolic crystal.

Figure 5:

Quantum interference (Q) for air/graphene/MoO₃/Si heterostructure with MoO₃ thickness of 50 nm, plotted as a function of frequency, for two different dipole heights of z ₀ = 50 and 200 nm, respectively. The color bar represents unitless quantity quantum interference (Q).

Integration of graphene at chemical potentials of 0.025 and 0.05 eV suppresses the quantum interference, as for the chosen dipole distances, the emitter does not couple efficiently to the HPPPs as they occur at very large wave vectors (as seen in Figure 2) following the relation k _ρ ∼ 1/z ₀ and the graphene itself is isotropic. However, when the chemical potential is changed to 0.1 eV, 0.3 eV and 0.5 eV, we begin to see significant differences in the quantum interference, due to coupling with the HPPP modes. For z ₀ = 50 nm, we see large positive quantum interference just below 545 cm⁻¹. This is explained by the avoided crossing between the graphene plasmon and phonon polaritons in RB-1 region, when the plane of propagation is along y–z plane (i.e. ϕ = 90°). However, this phenomenon does not occur in the x–z plane (i.e. ϕ = 0°) for frequencies below RB-1. Thus, the y-dipole sees an avoided crossing unlike the x dipole, and we get large positive quantum interference. Furthermore, a shift in the frequency band of large quantum interference is observed with increasing chemical potential of graphene from 0.025 to 0.5 eV inside RB-1 and RB-2 spectral region which is due to active tuning of HPPP modes, as discussed above.

Moreover, for z ₀ = 200 nm, we observe a similar phenomenon at μ = 0.1–0.5 eV. As we increase μ, the SPPs mode moves to a lower k _ρ/k ₀. Since large k _ρ/k ₀ modes are tightly confined, they remain unavailable to dipoles placed far away from heterostructure. Hence, this lower wave-vector modes at large μ are more prominent for a dipole at z ₀ = 200 nm as can be observed in Figure 5(b). We also observe that the frequency corresponding to large negative quantum interference shifts at the upper end of the RB-1 which, as described above, is due to the direction dependent avoided crossing at the lower end of the RB-2 for one emitter polarization as opposed to coupling with the isotropic graphene SPP outside the RB-1 for the other polarization. Thus, we see a large value of negative quantum interference just below RB-2, instead of a nearly zero value after RB-1. We also see that the large negative quantum interference for most of the RB-1 region in the no graphene case is replaced by a small positive quantum interference when 0.3 eV or 0.5 eV graphene is introduced. This can be understood from Figure 3 where we observe that while there are no available modes in x-direction without graphene (hyperbola), the lemniscate shape of the fundamental HPPP mode with μ = 0.3 and 0.5 eV for graphene provides available modes to both dipole orientations. We also see a corresponding reverse phenomenon in RB-2 region in Figure 5, which has a similar explanation.

In summary, we observe that our proposed platform provides a convenient way to achieve electrostatically tunable quantum interference in a valley exciton material placed nearby via formation of anisotropic hybrid plasmon-phonon polariton modes which couple selectively with dipoles oriented along x- and y-direction.

Model for calculation of polariton dispersion

We used the standard expression for complex valued 2D graphene conductivity derived from the local limit of the random phase approximation [30, 75]

Here, the first-term represents the intraband contribution to graphene conductivity and the second term is the interband contribution. e, k _B and ℏ refer to the charge of electron, Boltzmann’s constant and reduced Planck constant, respectively. We considered the value of temperature as T = 300 K and relaxation time τ = 10⁻¹³ s which matches well with experimental results [75]. Here μ refers to the and the chemical potential of graphene chosen (in eV). Ω = 200πcω is the angular frequency in rad/s (where ω is the frequency in cm⁻¹). This conductivity expression has been plotted for various chemical potentials (μ) in the supplementary document (Figure S1). Dielectric permittivity of graphene is derived using the Drude model to be [76]

where ϵ _gr is the relative permittivity of graphene in the in-plane directions, ϵ ₀ is the permittivity of free space, and the thickness of the monolayer graphene (d) is taken 0.34 nm in accordance with experimental results [77]. Since the out-of-plane conductivity of graphene is zero, the relative permittivity of graphene in the out-of-plane direction is unity.

The permittivity of α-MoO₃ in mid-IR spectral region is given by the Lorentz oscillator model. Frequency-dependent dielectric permittivity along each crystal axis is written as

The Table shows the values we used for the parameters ϵ _∞, ω _LO, ω _TO and γ along each principal axis, as proposed in [46].

The respective expressions for permittivities along principal directions are plotted in the supplementary information (Figure S1).

Permittivity tensor for any particular plane of propagation (i.e. for any ϕ) is obtained from ϵ(ϕ, ω) = R(ϕ)ϵ _principal(ω)R ^T (ϕ) where

is the rotation matrix for transformation of coordinate system, and ϵ _principal(ω) is the permittivity tensor in principle axis coordinate system [40]. This leads to the final expression

Fresnel’s reflection coefficient for p-polarized light is evaluated from well-known 4 × 4 transfer matrix method [40].

Calculation of Purcell factor

Purcell factor for an emitter placed at a height z ₀ above a surface in terms of the Green’s Dyadic function can be written as [41], [42], [43]

where μ → is a unit vector along the direction of the dipole, G = G ( r 0 → , r 0 → ) is the scattering part of the Green’s dyadic function at the position of the dipole, and k 0 = ω c . Here, P and P ₀ are the spontaneous emission rates of a dipole placed above the heterostructure and in free space, respectively.

From [41, 43] we can express the Green’s tensor as

where the integral is taken over all space and ‘M’ matrices are given by:

From here, we can recast the integral in polar coordinates and calculate the spontaneous emission rate enhancements.