Dynamically controlling local field enhancement at an epsilon-near-zero/dielectric interface via nonlinearities of an epsilon-near-zero medium

1 Introduction

An epsilon-near-zero (ENZ) material is a natural or engineered (meta)material with the real part of the dielectric permittivity approaching zero. In what follows, we will consider an ENZ medium to be a material that exhibits a dispersion curve with the real part of dielectric permittivity crossing zero, rather than just approaching it from either positive or negative half-planes. The interest in ENZ media has grown substantially during the last decade since the seminal works by Ziolkowski [1] and Silveirinha and Engheta [2]; they theoretically presented unusual modes of propagation and scattering of electromagnetic waves in ENZ materials. In such a material, the wavelength of light is stretched, and spatial and temporal field oscillations decouple, leading to effects such as zero propagation-phase advance and energy squeezing through super coupling—propagation of light through a narrow ENZ-material-based channel, highly directional emission, nonlinearity enhancement and others [2], [, 3]. Recent advancements in this field have given rise to a new promising technological platform dubbed near-zero refractive index photonics [3]. Among naturally occurring materials, transparent conductive oxides (TCOs), such as indium tin oxide (ITO), titanium nitride (TiN), aluminum doped zinc oxide (Al:ZnO), and phononic materials such as silicon carbide (SiC) have been studied as ENZ media. It is worthwhile to note that our analysis is limited to materials that exhibit the ENZ property due to the presence of screened electron plasma, such as TCOs. In the case of phononic materials, the ENZ property is determined by the effective polarization due to the ionic lattice motion. Even though the response due to phonon resonances will obey equations very similar to those for electronic responses [4], these materials are outside the scope of this study. Here, we focus on one particular property of ENZ materials, namely local field enhancement inside an ENZ medium due to the continuity of the normal component of the displacement vector at the interface with a conventional dielectric material which has a positive value of the real part of the dielectric permittivity. The continuity condition reads as follows:

where subscript 1 refers to the dielectric, subscript 2 to the ENZ medium and nˆ is the boundary unit normal. For linear lossless optical media, Eq. (1) translates to the following equation:

where Eiz is the normal component of the electric field of p-polarized wave and εi is the dielectric constant, which in the case of a lossless medium is real-valued. From Eq. (2), it follows that a zero-valued dielectric constant of material 2 gives rise to a large discontinuity of the normal component of the electric field across the boundary and, consequently, a divergent local electric field inside the ENZ medium. However, if the ENZ medium is lossy—i.e. the imaginary part of ε2(ω) has a nonzero value, ε2″—then the normal component of the electric field inside the ENZ medium attains a finite value (along with a phase shift)

For example, Alam et al. [5] have estimated that the local field of a 310-nm-thick ITO film interfaced with air is only enhanced by 2-fold at 1240 nm, which is the bulk plasmon wavelength of ITO. At the same time, the authors observed a large enhancement of the nonlinear response. This enhancement, as measured by the intensity dependent refractive index, was attributed in part to a simple argument: for a lossless material, one finds that Δn=Δε2√ε and thus the change in refractive index (Δn) for a given change in permittivity (Δε) has a pole under ENZ condition. Thus, it seems that in an ENZ material, even a very low-power optical field would produce a large nonlinear response, and local field enhancement would be one of the sources contributing to this effect. However, it is not entirely clear what is the contribution of local field enhancement to the macroscopic polarization of an ENZ medium, if we take (nonlinear)saturation effect into account. With a few tens of GW/cm² of incident peak intensity, this effect may prove significant and needs a thorough investigation. We believe that this knowledge not only can help interpret experimental results but also can aid in the design of ENZ (meta)materials with huge local field enhancement at a prespecified wavelength for various applications in nonlinear photonics. With the advent of inverse design tools based on machine learning and artificial intelligence, targeted design of metasurfaces with prescribed dispersion of macroscopic quantities, such as transmittance and reflectance, is becoming widely available to researchers [6]. In the inverse design paradigm, the macroscopic parameter space, describing the desired functionality of a new material, is mapped into the microscopic parameter space, describing the constituting engineered meta-atoms. The microscopic parameters can then be fed into our model which, in turn, can predict the dispersion of the local field enhancement. Subsequent iterative feedback cycles of prediction–adjustment–prediction can greatly accelerate the engineering efforts.

In what follows, we analyze the local field enhancement (or the squared-field enhancement factor [SFEF]) generated at the interface between a dielectric and a semiconductor ENZ medium represented by a model of general quantum two-state systems (Figure 1). We assume that our model material possesses second- and third-order nonlinearities and is driven by the optical fields strong enough to be able to generate an appreciable excited state population. It is well known that the real and the imaginary parts of the dynamic response of a two-state system to a monochromatic excitation follow a Lorentzian dependence on the frequency of the excitation. Hence, an ensemble of such systems can be a good approximation to materials with generalized Drude–Lorentz-type dispersion with the number density of carriers, N, corresponding to the number of two-state systems in the ensemble. The optical absorption mechanisms for such materials are either free carrier absorption (FCA) for wide-bandgap semiconductors or interband and intersubband absorption for narrow-bandgap semiconductors. In such a system, phonon-assisted FCA can be significant in both the visible and the infrared range [7]. Although the two-state approximation in many ways is quite crude, it is simple enough that we should be able to at least pinpoint some major trends in the optical response of Drude–Lorentz-type material systems.

Figure 1:

Conceptual sketch of a quantum two-state system comprising lower level 0 and upper level i.

VB = valence band, CB = conduction band, IBA = interband absorption, FCA = free carrier absorption.

2 Model

We consider the interaction of a monochromatic optical field, E=E2exp[−i(ωt+kz)]+c.c., with a quantum two-state system (Figure 1). In SI units, the power series expansion of the macroscopic polarization in the electric field strength up to the third order reads as follows:

where χ⁽ⁿ⁾ is macroscopic (nonlinear) susceptibility tensor. Here, we assume a single scalar susceptibility term which is applicable to isotropic media and may also be applied to a crystal with a cubic symmetry. In an anisotropic medium, such as a crystal with a lower symmetry, one has to consider the full susceptibility tensor. The macroscopic polarization of a perturbed quantum system interacting with a bath is best described by the density matrix formalism. The macroscopic polarization of an ensemble of quantum systems—expressed as the atomic number density N times the expectation value of the dipole operator dˆ—can be written in terms of the density matrix ρ(t) as follows:

where ρ_αβ are the off-diagonal elements of the density matrix. The orders of (nonlinear)susceptibilities are then found as follows [8]:

where rαβ(i) is the orders of a perturbation expansion of the density matrix:

The Liouville–von Neumann equation, governing the temporal evolution of the density matrix, reads in the interaction representation:

where Tr(ρ) = N, V = E⋅d and Γˆ is the so-called relaxation matrix, containing the damping parameters due to both the interactions of the quantum system with vacuum fluctuations (spontaneous decay) and with the bath (dephasing of quantum oscillators). Equation (8) for diagonal and off-diagonal elements of the density matrix can be readily solved in closed forms if we assume quasi-continuous wave (CW) excitation: the optical pulses should be longer than all characteristic relaxation times. This condition can be roughly fulfilled for laser pulses as short as a few nanoseconds:

The perturbation orders of the off-diagonal elements of the density matrix and the excited state population are then readily found:

where Gαβ=Edαβ2ℏ, ω is the excitation frequency and ω_0i is the resonance frequency of the two-state system. Expansion of macroscopic polarization (4) leads to a similar expansion of the displacement vector:

Applying the continuity across the boundary leads to the following equation:

where ε2eff(E2z) is the electric field–dependent effective dielectric constant due to nonlinear responses of the material up to the third order. Equation (12) can be solved numerically, in an iterative fashion, assuming the local field inside the ENZ medium equals the incident field at the first step of the numerical procedure. Note that we did set χ⁽²⁾ to zero in the simulations to reflect the fact that, in Eq. (12), the term in the denominator, that is proportional to χ⁽²⁾, is a fast oscillating term, so the averaging would kill it even if the susceptibility itself is not zero-valued.

3 Results and discussion

Before we solve Eq. (12), we need to choose parameters of our system to find the perturbation orders of the density matrix expansion (10). These important parameters are the permanent dipole moments of states 0 and i; the transition dipole, d_0i; the transition energy (bandgap); diagonal and off-diagonal elements of the relaxation matrix and the number density of free carriers, N. The free carrier number density can vary widely in materials of interest from ∼2.5 × 10²⁶ m⁻³ in doped TCOs such as In:CdO and ITO to ∼2.5 × 10²⁸ m⁻³ in TiN to ∼4.6 × 10³⁰ m⁻³ in gold. In our model, we chose the free carrier number density to be that of TiN.

We set the rate of spontaneous decay (diagonal element of the relaxation matrix) to 10⁹ s⁻¹ which is consistent with typical measured carrier lifetimes of ∼1 ns [9]. Exciton lifetimes can be even longer because of self-trapping effects, but we will use the referenced value, bearing in mind that longer lifetimes will simply result in a lower saturation intensity of the 0 → i transition, without affecting the spectral profiles, at least in the first order of approximation.

The dephasing rates of the carriers are highly dependent on the actual dephasing mechanism. Specifically, the carrier–carrier scattering times have been measured at tens of femtoseconds (GaAs, 3–50 fs) [10], [, 11]; exciton–exciton scattering times and free carrier–exciton scattering times have been reported to be hundreds of femtoseconds at moderate number densities. Exciton–phonon scattering times have been determined to span ∼1 ps (GaAs, GaP) [12], and exciton–impurity scattering times are the longest, measuring a few picoseconds. Obviously, the shortest dephasing time will translate into the largest dephasing rate, which will govern the spectral width of the corresponding off-diagonal density matrix frequency dependences. We chose this number to be 10¹⁴ s⁻¹ in our modeling.

Finally, the dipole moments of the excited state i and of 0 → i transition are chosen to be 500 and 10 Debye, respectively. The value of the transition dipole moment may be severely underestimated, considering the fact that in semiconductors the exciton bound to impurities and defects can possess giant oscillator strength, which can translate to the values of dielectric permittivity up to 10⁵, for example, in MXenes [13]. The value of 10 Debye can be closer to typical values of molecular-type (Frenkel) excitons. At the same time, some crystalline materials or materials containing transition metals can have similar characteristics. The value of the permanent dipole moment value of 500 Debye was chosen in conjunction with the fact that Wannier–Mott excitons, and especially 2D excitons, can have giant permanent dipole moments, up to 10³ Debye, because of the large electron–hole separation [14], [, 15]. The bandgap is set to 1 eV.

The choice of model parameters requires an important caveat: We did not try to model any specific material, but rather a class of materials with similar electronic properties, amenable to the analysis described in Section 2. The rather large pool of existing materials in this one class translates into a correspondingly wide range of material parameters that are primarily determined by the elemental composition. This sets the lower and upper limits, and our parameters are chosen to be somewhere in the middle of this range.

The main results of our modeling are presented in Figure 2(a) and (b). Following the study by Alam et al [5], we define the SFEF as |E2|2/|E1|2, where E₂ is the field strength inside the ENZ layer and E₁ is the incident field. For incident field intensities less than 10 MW/cm², the SFEF is much less than 1 and is spectrally centered at the resonance frequency of the two-state system. At such a low incident intensity, saturation and other nonlinear effects are negligible and the imaginary part of the dielectric constant suppresses the local field at the wavelength corresponding to zero value of the real part (Figure 3(a)). At around 20 MW/cm² of the incident field intensity, the SFEF reaches the value of 0.5 (still less than 1), also showing a blue shift in its peak position. The corresponding imaginary part of the dielectric constant shows a dip (Figure 3(b)) at the frequency of the SFEF peak and the real part flattens out around this peak position, although never reaching zero. Also, the maximum value of the imaginary part decreases in comparison with the lower intensity situation, depicted in Figure 3(a). This can be explained by an onset of saturation of the 0 → i transition. The effect becomes even more prominent at higher intensities. At 50 MW/cm², one can observe 4 enhancement peaks: one is rather broad at around 1.1 eV and three sharp ones at 1.26, 1.30 and 1.34 eV. The corresponding real part of the dielectric constant crosses zero at these frequencies (as shown in the inset in Figure 3(c)), and the imaginary part features a deep and spectrally wide saturation signature.

Figure 2:

Spectral dependence of the squared-field enhancement factor (SFEF) for (a) moderate incident field intensities and (b) higher incident field intensities.

Figure 3:

Spectra of the real and imaginary parts of the effective dielectric constant.

At even higher intensities, one can observe the onset of nonlinear saturation effect at the frequency of two-photon resonance that is around 0.5 eV (Figure 3(d)). The SFEF starts growing consistently at this frequency until it saturates at around 1 GW/cm² of the incident intensity (Figure 2(b)). Moreover, the contribution of the third-order nonlinear term to the effective dielectric constant (as defined by Eq. (12)) becomes visible at ∼0.5 eV (Figure 4). In Figure 4(a), one can observe a slight decrease of the real part of the effective permittivity and an increase of the imaginary part (two-photon absorption) at 100 MW/cm². At even higher incident intensities (∼500 MW/cm²), the contribution of the nonlinear absorption and refraction to the field enhancement around the frequency of the two-photon resonance seems to be on par with the saturated linear term of the effective dielectric constant. Both the real and the imaginary parts of the effective dielectric constant show a decrease in a large spectral interval. This decrease is more prominent for the real part of the effective dielectric constant (Figure 4(b)). Therefore, the mechanism of the field enhancement at two-photon resonance is essentially the same as in the situation when the excitation frequency is close to one-photon resonance: suppression of both real and imaginary parts of the dielectric constant because of the nonlinear saturation effects plus the contribution of the third-order term to the effective dielectric constant.

Figure 4:

Spectra of the real and imaginary parts of the linear term (solid line) and full effective (broken line) dielectric constant as defined by Eq. (12).

We also studied the effect of the carrier density on the field enhancement. The carrier density can be manipulated through doping and, in this context, is not related to the laser-induced population of the conduction band. The results are shown on Figure 5(a–c). At lower carrier densities (between 3 × 10²⁷ m⁻³ and 6 × 10²⁷ m⁻³, Figure 5(a)), we can see a very moderate enhancement peak around the frequency corresponding to resonant excitation and a peak around two-photon resonance at 0.5 eV. It seems that the red-shifted broad features between 3 and 6 eV are not physical for the following reason: the applicability of the two-state model is seriously compromised when the excitation frequency is much larger than the bandgap. In a more realistic many-level system, there would be a resonance with one of the low-lying excited states, with corresponding damping of the field enhancement. Even when the incident intensity is low, this enhancement increases far from the resonance, simply because in the Lorentzian two-state system the imaginary part of the dielectric constant decreases as 1/ω² which is faster than the corresponding decrease of the real part –1/ω. But it does not mean that such an enhancement is physical. In this context, the Rabi frequency can be a logical and natural cutoff metric. In our situation, it roughly measures 10% of the bandgap frequency at the carrier density of 2.5 × 10²⁸ m⁻³ and the incident intensity of 40 MW/cm². Figure 5(b) shows a close-up at the frequency of two-photon resonance: the enhancement decreases as the carrier density increases. The reason for this is that the real part of dielectric constant grows faster with the carrier density than it gets suppressed because of the saturation at the given incident intensity. At higher carrier densities (between 1 × 10²⁸ m⁻³ and 3 × 10²⁸ m⁻³, Figure 5(c)), the nonlinear effects are of course more prominent at the given incident intensity which results at larger values of field enhancement. Again, the real part of the dielectric constant grows faster with the carrier density than it gets flattened and suppressed because of the nonlinear contributions at the given incident intensity.

Figure 5:

Squared-field enhancement factor vs carrier concentration at an incident field intensity of 40 MW/cm².

4 Conclusion

In this study, we showed, using the density matrix formalism, that local field enhancement at the interface between a lossless dielectric (air) and a lossy Lorentz–Drude-type medium with large carrier density can be readily manipulated via adjusting the incident field strength. The mechanism of the enhancement is mainly related to the saturation of absorption (including nonlinear saturation effect), resulting in the previously lossy medium becoming transparent at moderate light intensities. This work also provides, for the first time, an exciting prospect of external field responsive enhancement factor as well as optical loss to dynamically tune the ENZ behavior. Another feature revealed by this work is further tuning of enhancement with the free carrier density, which can be adjusted via doping or by the use of a nanocomposite effective medium. We believe that a proof-of-concept experiment should involve an observation of a nonlinear phenomenon that is not directly related to the ENZ material itself but rather is sensitive to the local field enhancement produced by the ENZ material. In this case, one does not face the problem of extracting the field enhancement from under the overall enhancement of the nonlinear signal of the ENZ material. For example, an aluminum-doped zinc oxide layer deposited on a glass substrate and excited by a wavelength in the vicinity of ENZ may generate a local field strong enough to enhance second-harmonic generation in an add-on layer of ZnO nanoparticles. One will then look for tunability of the enhancement of the second harmonic generation (SHG) signal from ZnO.