1. Introduction

Despite the ongoing interest in nonlinear optical materials and structures, driven by their numerous applications, our progress in engineering high optical nonlinearities has remained rather modest. Sometimes, new ideas emerge, but upon close inspection, it often becomes evident that great expectations lead to limited outcomes. Nevertheless, we continue to make incremental strides towards turning promises into realities.

Broadly speaking, two conceptually different approaches can be distinguished in this field. The first approach involves designing homogeneous (single-phase) materials, where intrinsic nonlinearity can be enhanced through various physical phenomena, including quantum effects. This entails molecular (subnano) scale and quantum-mechanical calculations of the local microscopic fields. The second approach focuses on various nanostructures and involves calculating the mesoscopic and macroscopic fields. Our understanding of the interplay between microscopic, mesoscopic, and macroscopic phenomena could help us to take a step forward for adequate interpreting and exploiting the physics of nonlinearity. Obviously, combining both approaches can yield the best outcomes.

Of particular practical importance is the degenerate third-order nonlinear susceptibility $\chi ^{(3)}(\omega =\omega +\omega -\omega )$, which links the nonlinear refractive index $n_2$ and the nonlinear absorption coefficient $\kappa _2$. In the case of weak nonlinearity, the third-order contribution modifies the refractive index and absorption coefficient as $n=n_0+n_2I$ and $\kappa =\kappa _0+\kappa _2I$, respectively.

Despite the significance of nonlinearity and nonlocality in this context, their relationship has not been adequately addressed. A quarter of a century ago, Sahimi, while studying transport phenomena in heterogeneous media, pointed out that "often non-linearity and non-locality are, in fact, two sides of the same coin," implying their interdependence [1]. Indeed, nonlocality is a fundamental characteristic observed in numerous nonlinear systems. In such systems, the response of the medium at a specific point is not solely determined by the wave intensity at that point alone. Instead, it is influenced by the wave intensity in its surrounding vicinity. This nonlocal behavior can arise from various transport processes, including atom diffusion, heat transfer, or the drift of free charge carriers.

In bulk materials, nonlocality effects are known to shift the main plasmonic resonance to higher frequencies, generate additional absorption resonances, and induce anisotropy [2]. Numerous scholarly works explore the effects of nonlocality on the modulation instability of plane waves and the generation of diverse types of solitons [3,4]. The above effects arise from intrinsic nonlocality or in nonlinear media with a spatially nonlocal nonlinear response.

Another significant body of literature is dedicated to the study of composite or inhomogeneous media. Of particular interest are metamaterials (MMs) - periodic nanostructures with specific microgeometry, which allow for the realization of on-demand optical properties. Notably, recent reviews by Smirnova and Kivshar [5] discussed various multipolar nonlinear effects in resonant plasmonic and all-dielectric photonic structures, and Kolkowski et al. [6] considered nonlinear nonlocal metasurfaces. The authors emphasized that bulk nonlinear polarization mainly arises from higher-order nonlocal magnetic dipole (MD) and electric quadrupole (EQ) interactions with light at the microscopic scale. They also noted that the conventional electric dipole (ED) approximation in the bulk of uniform centrosymmetric media forbids second-order nonlinear effects, but this restriction is lifted when higher-order multipoles come into play.

Wurtz et al. [7] and Neira et al. [8] demonstrated that, under certain conditions, the nonlinear response of plasmonic nanorod MMs can be significantly enhanced by appropriate design of MM geometry. Additionally, Vincenti et al. [9] recently studied a metal - transparent conducting oxide multilayer stack and showed that nonlocality and hot electrons compete to either blueshift or redshift the plasmonic resonances, consequently enhancing or suppressing the nonlinear response.

In this paper, we focus on the broadband enhancement of nonlinearity. The bandwidth plays a crucial role in the functionality and performance of many nonlinear optical devices, as it determines the range of frequencies available for interaction within the material, affects phase matching, and compensates for dispersion effects. For instance, in optical communication, broad bandwidth sources are needed to carry a large amount of information simultaneously, thus improving performance of high-speed data communication networks [10]. Nonlinear optical devices like optical parametric amplifiers and wavelength converters efficiently manipulate and process signals. Laser frequency conversion often requires a broad bandwidth light source to cover the desired range for specific applications. In multiphoton microscopy, a broad bandwidth is desirable to excite various fluorophores and achieve higher-resolution imaging. Broadband optical limiters can respond to a wide frequency range, making them effective against various light sources that might pose a threat to optical systems. Anyway, despite numerous efforts [11–13], it seems that broadening the operational bandwidth of nonlinear optical devices and nanostructures remains to be an unsolved problem which deserves more attention.

The relationship between nonlinearity and nonlocality is too intricate to be addressed in its entirety. As expected, neither model can comprehensively describe the interplay between both phenomena. Thus, in this paper, we focus on the simplest 1d case. Our study is based on the following premises: (i) We consider a periodic multilayer MM, whose unit cell can contain any number of constituents (phases). (ii) The number of constituents coincides with the number of layers within the unit cell. (iii) All constituents are non-magnetic ($\mu _i=1$) and can be described by a scalar complex permittivity $\epsilon _i$ and third-order nonlinear susceptibility $\chi _i^{(3)}$. (iv) The electric field within a unit cell is directed normally to the layers (TM polarization). (v) The thickness of layers is assumed to be at least a few tens of nanometers, thus neglecting the possible modification of their linear and nonlinear optical properties arising from quantum and dielectric confinements of free electrons in the layers [14].

Several key questions can be formulated within the framework of our approach: How to get broadband enhancement of nonlinearity with account for nonlocal effects dealing with 1d microgeometry and realistic material parameters? How to assess the contribution of nonlocal polarization to nonlinearity and how valuable can it be? What is the effect of nonlocality on the bandwidth and nonlinearity enhancement? Can nonlocality give rise to nonlinearity? For a given bandwidth, how to choose an optimal arrangement of the layers?

To address the above issues, in Sec. 2 we first present our main idea to broaden the bandwidth of enhanced nonlinearity and then outline our approach to calculate the effective third-order nonlinear susceptibility for a given 1d microgeometry. Section 3 and Sec. 4 contain the results of our computations for the direct and inverse problem, respectively. Section 5 discusses the results obtained and Sec. 6 concludes the paper with some final remarks.

2. Background

2.1 Broadband enhancement of nonlinearity

With reservation for a not too low linear refractive index and not too high optical field intensity [15], the third-order nonlinear refractive index $n_2$ and absorption $\kappa _2$ of macroscopically homogeneous absorbing media are known to be functions of the permittivity $\epsilon$ and the third-order nonlinear susceptibility $\chi ^{(3)}$ of the form [16]

The above mentioned idea has been extensively discussed, and experimental confirmation of this enhancement has been demonstrated [15,17–19]. However, the ENZ range of homogeneous materials is typically narrow, resulting in a narrow-band enhancement of nonlinearity, thereby limiting its practical applications. Therefore, a plausible solution is to broaden the bandwidth using MMs. The material and geometrical parameters of the metamaterial building blocks (unit cells or metamolecules) can be considered as degrees of freedom, and their appropriate selection could facilitate the design of broadband MMs with enhanced functionalities [20,21].

To treat the effective nonlinear response of composite media, analytical approaches, based on the exact calculations of the electric field inside isolated nonlinear inclusions, embedded in a homogeneous host, have gained acceptance. Typical examples include, e.g., arrays of anisotropic spheres [22] or shape-distributed core-shell ellipsoids [23]. These approaches, however, are accurate only in the dilute limit. A simple, general and rather accurate approach to calculating the nonlinear response of composite materials have been introduced by Stroud and Hui [24]. According to their theorem, an effective nonlinear response of a composite medium $\tilde \chi ^{(3)}$ through first order in $\chi _i$ can be represented as

Generally, to evaluate the effective nonlinear susceptibility, full-wave simulations can be adapted [25–27]. However, for 1d MMs consisting of subwavelength multilayers, the local electric field can be analytically calculated. In the particular case when the applied electric field is directed along the optical axis (TM polarization), Eqs. (3) and (4) take the forms [28–30]

Given the local permittivity $\epsilon$ and third-order nonlinear susceptibility $\chi ^{(3)}$ are continuous functions of coordinate $z$ within the lattice period $d$, Eq. (6) takes the form of

In some cases, $\tilde \epsilon$ and $\tilde \chi ^{(3)}$ can be obtained in a closed form. In particular, it takes place if $\epsilon$ is of the Drude form

Our analysis is based on the premise of a high free carrier density. In this case the free-carrier contribution to the local third-order nonlinear susceptibility tends to saturate, especially at frequencies well below the band gap [32,33]. If so, it can be considered as independent of the free carrier density $N$, i.e., $\chi _i^{(3)}=\chi _0\approx const$. It is then convenient to introduce the enhancement of the nonlinear susceptibility $\mu =\tilde \chi ^{(3)}/\chi _0$, which can be evaluated as

Because the nonlinear susceptibility $\tilde \chi ^{(3)}$ is expected to be enhanced wherever $\epsilon \approx 0$, i.e., at $\omega =\omega _p/\sqrt {\epsilon _\infty }$, it follows from Eq. (6) as well as from Eq. (13) that it can be enhanced within the band $[\omega _{min}/\sqrt {\epsilon _\infty },\ \omega _{max}/\sqrt {\epsilon _\infty }]$. This allows us to introduce the related bandwidth as $\delta =(\omega _{max}-\omega _{min})/\omega _0$.

There are, however, two reservations which should be kept in mind when designing broadband nonlinearity enhancement in a similar way. First, while it is a straightforward task to get $\Re (\epsilon )\approx 0$ at particular frequencies, it is not easy to simultaneously maintain $\Im (\epsilon )$ small. It is obvious that dielectric loss can significantly reduce nonlinearity enhancement. Second, in the above approximation, both linear and nonlinear effective response are independent of the arrangement (relative disposition of the layers within the unit cell). In addition, as the lattice period $d$ rises, the local effective medium approximation is no longer valid. At first glance, to account for the effects of nonlocality, we could substitute into Eq. (8) the nonlocal effective permittivity instead of its local counterpart. Such an approach, however, does not look entirely correct.

Within the framework of the effective medium approximation, after replacing $\chi ^{(3)}$ in Eqs. (1) and (2) with $\tilde \chi ^{(3)}$, one has

While Eqs. (5 )–(9) are valid in the long wavelength limit only, the main results of this paper are obtained using Eq. (4), and involve the accurate calculation of the inhomogeneous electric field (see Supplement 1, Section 2).

2.2 Microgeometry

We deal with 1d multilayered materials whose basic microgeometry is schematically shown in Fig. 1. It is assumed to be infinite in the $x-y$ plane and periodic along the $z-$axis. This simplification allows one to avoid dealing with nonlocal boundary conditions.

 

Fig. 1. The schematic of the unit cell of MM under consideration

Download Full Size | PDF

While we show here a five-layer unit cell, it is not a critical limitation of our paper, as can be seen below. The filling factor of $i$th layer is $f_i$, its thickness is $d_i$, and its permittivity is $\epsilon _i$. The structure under study is assumed to be subwavelength, i.e., $d=\sum d_i < \lambda =2\pi c/\omega$. For simplicity, all the layers are considered to be nonmagnetic, i.e., their magnetic permeability is unity. The electric field is directed along the $z-$axis, i.e., normally to the layers, that corresponds to the TM-polarization, and a wave vector is directed along the $x-$axis.

2.3 Constituent materials

Compared to metal-based metamaterials (MMs), heavily doped semiconductors often exhibit a plasmonic resonance peak in the infrared (IR) range due to their lower carrier density. This distinctive characteristic opens up new avenues for a diverse range of infrared photonic applications that remain less accessible through conventional metal-based plasmonics. For instance, semiconductor-based MMs with broadened absorption bands [21,34,35], all-optical switching [36], or semiconductor nanostructures facilitating quantum engineering techniques [37] deserve special mention.

Moving forward, it becomes imperative to specify both linear and nonlinear material parameters, such as permittivity and third-order nonlinear susceptibility. In principle, we can consider both real and hypothetical constituent materials. Notable examples of semiconductor-based plasmonic MMs encompass heavily-doped silicon and transparent conductive oxides (TCOs), including indium tin oxide (ITO), aluminum-doped zinc oxide (AZO), and gallium-doped zinc oxide (GZO). Among these TCOs, AZO stands out due to its remarkably low losses, rendering it a fundamental material in our investigation. It is crucial to recognize, however, that precise calculations of permittivity and particularly the third-order nonlinear susceptibility of heavily doped semiconductors in the near-IR region present challenges due to the intricate interplay of diverse physical mechanisms. Moreover, available experimental data are sparse and occasionally conflicting. Consequently, both linear and nonlinear material parameters can exhibit variability contingent upon specific experimental conditions. Given this, the constitutive materials that we examine subsequently may be denoted as AZO-like hypothetical materials. Anyway, we consider our subsequent analysis as qualitative rather than quantitative.

There are a few simple yet accurate models for the dielectric response of doped semiconductors in the visible and IR ranges [38]. We use the phenomenological Drude model (10), which is known to be accurate for TCO, especially close to the epsilon-near-zero (ENZ) range [39–41].

The dielectric properties of TCO films can significantly vary based on the fabrication method employed and specific parameters applied throughout the fabrication process [40]. In this context, we consider specific values, namely $\epsilon _\infty = 3.85$, $\gamma = 0.1$ eV, and $\omega _p = 1.747$ eV, for a dopant concentration of 2 wt%. This parameter selection roughly aligns with prior findings by Naik et al. [42] and Tian et al. [43], focusing on AZO samples fabricated through pulsed laser deposition techniques. During the deposition process, the concentration of aluminum can be manipulated by adjusting the number of pulses onto the alumina target. For the sake of simplicity, we maintain a constant free carrier damping rate, $\gamma _p$, independent of the dopant concentration. Additionally, within the practical concentration range of 1.6-2 wt%, the square of the plasma frequency within the $i$th layer, $\omega _{pi}^2$, is presumed to exhibit a rough linear dependency on the dopant concentration.

At this stage, it would be best to specify a model for the intrinsic third-order susceptibility. At present, however, we must admit that physical mechanisms of optical nonlinearity are poorly studied. There are numerous physical mechanisms that could contribute to the third-order nonlinear susceptibility of heavily-doped semiconductors in the near-IR. At photon energies much smaller compared to the band gap, the conduction-band nonparabolicity is considered to be the dominant mechanism of the free-carrier optical nonlinearity [44]. Other contributions result from the energy dependence of the relaxation processes [45].

In heavily-doped semiconductors there is a high concentration of free carriers that can interact with the incident electric field, leading to a nonlinear response [46]. This nonlinearity arises due to the third-order nonlinear susceptibility of the carriers, which is related to the Pauli blocking effect, band filling, and carrier-carrier interactions.

In the case of the free carrier nonlinearity, the Pauli blocking effect and band filling refer to the modification of the carrier distribution in the conduction band due to the presence of a high carrier density. As the carrier density increases, the Pauli exclusion principle comes into play, leading to a modification of the available states for carriers and affecting their response to the electric field. Band filling refers to the filling of these available states, which can further modify the nonlinear response.

Carrier-carrier interactions can lead to energy exchange processes, scattering, and modification of the carrier distribution, affecting the overall nonlinear response. The other example involves phonon-assisted nonlinearities. These nonlinearities can arise due to the anharmonic nature of the lattice vibrations. They can result in a third-order susceptibility that is related to the phonon population and temperature.

Furthermore, defects such as vacancies and interstitials can contribute to the third-order nonlinear susceptibility. This arises due to the nonlinear response of the defect states to the incident electric field.

Finally, mention must be made of multiphoton absorption. It takes place as the absorption of two or more photons occurs simultaneously, leading to the excitation of carriers to higher energy states. These highly energetic carriers, and especially hot electrons, can participate in the nonlinear response of the material. The probability of multiphoton absorption depends on the intensity and wavelength of the incident laser pulses, as well as on the electronic band structure of the material. Hot carriers can undergo various relaxation processes, such as carrier-carrier scattering, phonon interactions, and impact ionization, which can influence the nonlinear response of the material.

In practice, the dominant contribution should depend on several factors, including the specific experimental conditions and the material properties. The free carrier nonlinearity is expected to be one of the dominant contributions in TCO at room temperature and given doping levels. The contribution of the phonon-assisted nonlinearities and defect-related nonlinearities may also be significant, but their relative importance depends on the specific properties of the material and the experimental conditions. Two-photon absorption is typically less important at these wavelengths, but three-photon and four-photon absorption may be substantial.

In ITO and AZO, optical nonlinearity is believed to occur mainly due to hot electrons [17,47]. It results from a modification of the energy distribution of conduction band electrons, which in turn is resulted from laser-induced electron heating. At the same time, in heavily-doped GZO, the hot-electron energy relaxation time is shown to depend on the electron density (doping level) [48].

The precise contribution of each mechanism can be difficult to determine experimentally, as the relative contributions of different mechanisms can depend on the specific material, doping level, and experimental conditions. Overall, a comprehensive understanding of the nonlinear response in heavily-doped semiconductors requires considering all these mechanisms and their interplay. Currently, we lack a detailed understanding of the material properties and the underlying physical mechanisms which are required to accurately model and interpret experimental data.

3. Simulations and results

3.1 Long wavelength limit

In this limit it is assumed that $\Lambda \equiv d\omega /2\pi c \ll 1$ and hence Eqs. (5 )–(9) and (12), (13) are valid. Although Eqs. (7) and (9) can be used only within their range of applicability and have been criticized many times due to their lack of accuracy [49], they can be considered as initial approximations.

It seems useful to show the frequency dependence of $\mid \tilde \epsilon \mid$, as well as the functions $\epsilon _i$, which enter Eq. (7). We plot them in Fig. 2 for the seven-layer unit cell (because all $\Im (\epsilon _i)$ are almost constant within the targeted band, only $\Re (\epsilon _i)$ are shown here). The distribution of $N_i$ is step-like with $N_{i+1}-N_i=const$. The bandwidth, marked by double arrows, is [$\lambda _{min}=2\pi c/\omega _{max}=1433$ nm, $\lambda _{max}=2\pi c/\omega _{min}=1582$ nm], with the central wavelength of $\lambda _0=2\pi c/\omega _0 \approx 1501.5$ nm, that yields $\delta \approx 0.099$.

 

Fig. 2. The absolute value of the effective permittivity (bold line) and the real parts of layers permittivities (dashed lines) for the seven-layer unit cell ($f_i=1/7$) with uniformly distributed $N_i$ and the related bandwidth $\delta =0.099.$

Download Full Size | PDF

It is easy to make sure that $\Re (\tilde \epsilon )$, though not shown here, is a monotonic decreasing function of $\lambda$, which crosses zero at one point. Specifically, $\Re (\tilde \epsilon )=0$ at $\lambda \approx 1506$ nm. As $\Im (\tilde \epsilon )$ is nearly maximum at this point, min$(\mid \tilde \epsilon \mid )$ is shifted to lower wavelengths.

It is also of interest to compare the results, obtained for $\mid \tilde \epsilon \mid$ with the use of step-like distribution of the free-carrier density, Eq. (7), and its continuous distribution, Eqs. (9 )–(12). Our computations, performed for two bandwidths, $\delta =0.066$ and $\delta =0.099$, presented in Fig. 3. Doing so, we note that the continuous distribution of the form (11) is most similar to the step-like distribution ($n$-layer unit cell) with $f_1=f_n=1/2(n-1)$ and $f_2=f_3=\cdots =f_{n-1}=1/(n-1)$, as shown in the inserts of Fig. 3 for the cases of $n=5$ and $n=7$.

 

Fig. 3. The absolute value of the effective permittivity for the five-layer (left panel) and seven-layer (right panel) unit cells (bold lines). The dashed lines display the same for the corresponding continuous distributions of the free-carrier density, shown in the inserts by the dotted lines.

Download Full Size | PDF

The spectra of $\mid \tilde \epsilon (\omega )\mid$ for the continuous distributions of the free-carrier density are almost indistinguishable from those for seven-layer unit cells (step-like distributions) for the $6.6{\% }$ and $9.9 {\% }$ bandwidths, so they are not shown in Fig. 3 (the right panel). As can be seen, while for the $n-$layer unit cell $n$ oscillations can occur in the spectra, the continuous distribution gives rise to the smooth spectra of $\mid \tilde \epsilon (\omega )\mid$. Moreover, the values of $\mid \tilde \epsilon \mid$ are almost constant within the bandwidth.

The main conclusion that can be drawn from Fig. 3 is that we can always reduce the amplitude of the oscillations in the spectra of the linear effective response by increasing the number of the layers within the unit cell. What happens then with the nonlinear effective response? To address this question, in Fig. 4 we show the enhancement of the third-order nonlinear susceptibility $\mu$ for two bandwidths, $\delta =0.066$ and $\delta =0.099$, for the five-layer and seven-layer unit cells. The filling factor of each layer is $f_i=1/5$ in the former case and $f_i=1/7$ in the latter case. Similar to the spectra of $\mid \tilde \epsilon \mid$ shown in Fig. 3, the spectra of $\mu (\omega )$ exhibit $n$ oscillations within targeted bands, but their magnitudes become more pronounced. Additionally, we note that as the value of $\delta$ approaches zero, the multilayered material becomes homogeneous, leading to the trivial outcome of $\mu =1$. In contrast, when the bandwidth $\delta$ increases, the absolute value of $\mu$ also tends to increase.

 

Fig. 4. The enhancement of the third-order nonlinear susceptibility for the five-layer (left panel) and seven-layer (right panel) unit cells for two bandwidths, $\delta =0.066$ and $\delta =0.099$.

Download Full Size | PDF

3.2 Effects of nonlocality

The effects of nonlocality can manifest themselves exclusively through retardation, but also implicitly through various arrangements (layouts) of the layers. In the former case, this only pertains to the dependence of $\mu$ on the unit cell size (lattice period) $d$, while in the latter case it also involves the dependence of $\mu$ on the layer disposition. In both instances, the emergence of nonlocal polarization is driven by higher order multipoles coupled with retardation effects [50–52].

It would be natural to consider the effect of nonlocality in terms of dimensionless parameters $\Lambda =d/\lambda$ or $k_0d=2\pi \Lambda$. However, from the practical point of view, as the targeted band is fixed, we prefer to operate with explicitly specified lattice period $d$. Specifically, at $d=$500, 700, and 1000 nm this yields $\Lambda _0=d/\lambda _0\approx$ 0.33, 0.47, and 0.67, respectively.

For simplicity, in the following discussions, we consider $n=5$ and $\delta =0.066$. To constrain the parameter space, we set all values of $f_i$ to be a constant, specifically $f_i=1/5$. Regarding the layer dispositions, there exist a total of 12 independent permutations of the layers. These permutations can be denoted as ’12345’ (representing the linear, monotonic distribution of $N_i$) and ’12354’, ’12435’, ’12453’, ’12534’, ’12543’, ’13245’, ’13254’, ’13425’, ’12524’, ’14235’, ’14325’ (representing nonmonotonic distributions). The numerical sequences here indicate the order in which the layers are arranged within the unit cell.

Consider first the behavior of the enhancement of the third-order nonlinear susceptibility $\mid \mu \mid$ at various $d$. The spectral profiles of this function at $d=$ 500 nm, 700 nm, and 1000 nm are illustrated in Fig. 5 for two arrangements, ’12345’ and ’14235’. The effective nonlinear third-order susceptibility is determined using Eq. (4), while the mesoscopic electric fields are computed as outlined in Supplement 1. For comparison, Fig. 5 also presents $\mid \mu \mid$ calculated in the long wavelength approximation ($d\to 0$). Remarkably, for the linear arrangement (’12345’), all four curves intersect at three specific wavelengths: $\lambda \approx 1490$ nm, $\lambda \approx 1501$ nm, and $\lambda \approx 1514$ nm. This implies that $\mid \mu \mid$ remains almost independent of $d$ within the range of 1490 to 1514 nm. On the other hand, in the case of the ’14235’ arrangement, the four curves intersect at $\lambda \approx 1488$ nm, $\lambda \approx 1536.5$ nm, and $\lambda \approx 1544.5$ nm.

 

Fig. 5. The absolute value of the enhancement of the third-order nonlinear susceptibility for the five-layer unit cell with the bandwidth $\delta =0.066$. The arrangement is ’12345’ (left panel) and ’14235’ (right panel).

Download Full Size | PDF

The 2d plot of the electric field intensity $\mid E(z)/E_0\mid ^2$ in the layers at $d=$1000 nm, presented in Fig. 6, shows the enhanced magnitude of the electric field at specific wavelengths for particular layers.

 

Fig. 6. The 2d plot of the electric field intensity in the layers for the five-layer unit cells at $d=1000$ nm (arrangement of ’12345’).

Download Full Size | PDF

The effect of layer arrangements on $\mid \mu \mid$ is depicted in Fig. 7. While the presented data is for $d=700$ nm, it is worth noting that the qualitative trends hold for $d=500$ nm and $d=1000$ nm as well. Specifically, the high-frequency peak value experiences a decrease of 3.8% - 6.4% at $d=500$ nm, 7.6% - 13.5% at $d=700$ nm, and 12.8% - 25% at $d=1000$ nm for arrangements of ’14325’ and ’13245’, respectively. Concurrently, the low-frequency peak value increases by 5.4% - 8.3% at $d=500$ nm, 11.1% - 16.8% at $d=700$ nm, and 24.3% - 36.9% at $d=1000$ nm for arrangements of ’13254’ and ’12432’. As can be seen, the linear distribution of $N_i$ (arrangement of ’12345’) yields a value of $\mid \mu \mid$ close to its minimum at the high-frequency edge of the band, while providing an almost maximum value of $\mid \mu \mid$ at its low-frequency edge. The ’14325’ arrangement offers the highest enhancement of nonlinearity at the high-frequency band edge and a substantial, though not maximal, enhancement at the low-frequency edge. Similarly, the ’12435’ arrangement delivers the greatest nonlinearity enhancement at the low-frequency band edge and a notable, yet not the highest, enhancement at the high-frequency edge.

 

Fig. 7. The absolute value of the enhancement of the third-order nonlinear susceptibility for the five-layer unit cells with the bandwidth of $\delta =0.066$ at $d=700$ nm. The bold red curve corresponds to the linear distribution of $N_i$ (arrangement of ’12345’) and the bold blue curve corresponds to the long wavelength approximation, $d \to 0$. In order not to overload this figure, the results for all other 11 arrangements are shown as the dashed black curves.

Download Full Size | PDF

It should be noted that this study does not aim to delve into an exhaustive analysis of the spectra of $\mu (\omega )$ for all possible arrangements. In the case of the five-layer MM, there exist 12 independent arrangements, with two of them possessing symmetry properties. Specifically, both the ’12345’ and ’14325’ arrangements are anti-symmetric. Here, anti-symmetry implies that $\Re [\epsilon (z)-\epsilon _0]=-\Re [\epsilon (-z)-\epsilon _0]$ and, concurrently, $\Im [\epsilon (z)]=\Im [\epsilon (-z)]$. To address the question how anti-symmetry manifests itself in the spectra of $\mu (\omega )$, Fig. 8 displays the (relative) enhancement of nonlinearity, normalized to its value in the long wavelength limit, $\mid \mu (\omega,d)/\mu _0(\omega )\mid$, where $\mu _0(\omega )=\mu (\omega,d\to 0)$, in the vicinity of the band center $\omega _0$. Notably, for the linear ’12345’ arrangement, there is a sub-band of [1490 nm - 1515 nm], where the deviation of the nonlinearity enhancement $\mid \mu (\omega, d)\mid$ from $\mid \mu _0(\omega )\mid$ remains below 1% at $d=$500 nm.

 

Fig. 8. The relative enhancement of the third-order nonlinear susceptibility $\mid \mu /\mu _0\mid$ for the five-layer unit cells, arrangements of ’12345’ and ’14325’, at $d=$500, 700, and 1000 nm. The bandwidth is $\delta =0.066$.

Download Full Size | PDF

It is interesting that the equality $\Im (\mu )=\Im (\mu _0)$ holds at $\lambda \approx 1494$ nm and $\lambda \approx 1508$ nm within this sub-band for all $d$ values, $d=$500, 700, and 1000 nm (not shown here). As a result, within the interval of [1594 nm - 1508 nm] the deviation of $\mid \Im (\mu )-\Im (\mu _0)\mid$ does not exceed 0.03 at $d=$500 nm, while $\Im (\mu _0)\approx 0$.

3.3 Impact of energy loss

Above we employed the standard Drude model with a fixed free carrier damping term $\gamma _p$ to set the permittivity of semiconductor layers. The damping term (electron collision rate) may be expressed as $\gamma _p=e/\mu ^*m^*$, where $e, \mu ^*$ and $m^*$ denote the carrier charge, mobility, and effective mass, respectively. Various scattering mechanisms contribute to this damping term, directly leading to light absorption or energy dissipation. The primary mechanisms are known to include: (i) Ionized impurity scattering: Predominant in heavily doped semiconductors, this mechanism involves the Coulombic interaction between free electrons and ionized impurities. (ii) Acoustic phonon scattering: When free electrons scatter with acoustic phonons, energy exchange occurs, contributing to energy dissipation. (iii) Polar optical phonon scattering: Electrons can scatter with optical phonons, associated with the vibrations of ions in the crystal lattice, leading to an overall reduction in electron mobility.

In many scattering events, the energy transfer is more evident in terms of momentum change rather than energy absorption. Nonetheless, the cumulative effect is an increase in resistance and a decrease in carrier mobility, influencing the overall electrical conductivity of the semiconductor. The damping term can generally depend on free carrier density, temperature, crystallinity, and experimental conditions.

While a detailed exploration of the impact of energy loss on the enhancement of nonlinearity is beyond the scope of our study, it can be crucial for optimizing the performance of nonlinear semiconductor optical devices. So, some remarks on this issue can be made. With this aim, in Fig. 9 we present the enhancement of the third-order nonlinear susceptibility $\mid \mu \mid$ calculated for three different values of the damping term. The baseline value, $\gamma _p=\gamma _p^*$, aligns with our previous calculations, while the other two values vary by 10%, representing a 10% reduction and a 10% increase in damping.

 

Fig. 9. The absolute value of the enhancement of the third-order nonlinear susceptibility for the five-layer unit cell with the bandwidth of $\delta =0.066$ at $d=700$ nm for various values of the damping term $\gamma _p$, arrangement of ’12345’.

Download Full Size | PDF

The results presented in Fig. 9 unequivocally demonstrate that nonlinearity decreases, and the oscillations of $\mid \mu \mid$ become smoother with an increase in $\gamma _p$. While we have showcased our calculations for the ’12345’ arrangement, we have made sure that results for other arrangements are qualitatively similar.

In general, both the real part, $\Re (\mu )$, and the imaginary part, $\Im (\mu )$, contribute to $\mid \mu \mid$. At the same time, analogous to Fig. 4, the contribution of $\Re (\mu )$ is most significant at the band center, while the contribution of $\Im (\mu )$ is most pronounced at the band edges.

4. Inverse problem in nonlinearity design

In the preceding section, our focus was on the direct computation of the enhancement of effective third-order nonlinear susceptibility. We established a versatile numerical technique, made initial assumptions, and applied this method to showcase the impact of nonlocality on the enhancement. Notably, these calculations did not involve any optimization. However, from a practical standpoint, the inverse problem becomes pivotal. This problem revolves around the optimal design of optical nonlinearity–specifically, determining the material and geometrical parameters that yield the desired performance of nonlinear materials and nanostructures.

The design process, in general, can be intricate. It typically involves several steps, including defining design goals based on initial specifications, converting these specifications into an objective function, determining the parameter space, and setting initial values for fitted parameters. The ultimate aim is to minimize the objective (cost) function with respect to the design parameters. Design specifications usually manifest as a set of costs and constraints. The output must adhere to all constraints while minimizing costs.

A comprehensive exploration of the inverse problem goes beyond the scope of this paper. Nevertheless, a few remarks on this issue are warranted. Clearly defining the design goal is imperative, especially when dealing with nonlinearity enhancement, as it is application-dependent. For instance, in applications like self-phase modulation, cross-phase modulation, and Kerr lensing in nonlinear optics, the real part of $\tilde \chi ^{(3)}$ holds significance. In telecommunications, it is crucial for phase modulation of optical signals in various communication systems. The imaginary part of $\tilde \chi ^{(3)}$ is often associated with two-photon absorption, playing a role in applications involving energy transfer processes, such as multiphoton microscopy and photodynamic therapy. The absolute value of $\tilde \chi ^{(3)}$ is vital for applications where the nonlinear refractive index is substantial, as it determines the strength of nonlinear effects.

To illustrate, let us consider the relatively simple task of designing nanostructure that maintains a high and nearly constant real part of the enhancement, $\Re (\mu )\approx \mu _0$, within a specified frequency band [$\omega _{min},\omega _{max}$] (note that if $\Im (\chi _0)$ can be neglected, then $\Re (\tilde \chi ^{(3)}) \propto \Re (\mu )$). This problem can be reduced to minimizing an objective function of the form $\int _{\omega _{min}}^{\omega _{max}}\mid \Re [\mu (\omega )]-\mu _0\mid d\omega$.

Multiparameter optimization often yields multiple solutions due to the presence of multiple local minima of the objective function. In Fig. 10 we present four designed dependencies $\Re (\mu )\ vs\ \lambda$ at $d=700$ nm and three designed dependencies at $d=1000$ nm, where we have set $\mu _0=3.5$, $\lambda _{min}=1455$ nm, and $\lambda _{max}=1555$ nm. As fitting parameters, we have taken the doping levels in layers of the unit cell. The optimization problem has been solved using the interior point algorithm.

 

Fig. 10. The fitted real part of the enhancement of third-order nonlinear susceptibility for the five-layer unit cell with the bandwidth $\delta =0.066$ at $\mu _0=3.5$ for two periods of the unit cell: $d=700$ nm (left panel) and $d=1000$ nm (right panel). The specific values of the doping level in each layer are shown in the square brackets.

Download Full Size | PDF

All solutions, obtained at $d=700$ nm, yield very close values of the objective function. Generally, a new solution can be obtained by repeating the optimization procedure after layer rearrangement, using the previous solution as a starting point.

5. Discussion

Initially, we adopted a uniform selection of $f_i$ values, resulting in a nearly flat dependence of $\mid \tilde \epsilon (\lambda )\mid$. However, in general, $\mid \tilde \epsilon \mid$ exhibits a slight increase with $\lambda$ within the band of interest, as illustrated in Fig. 2. Upon varying $f_i$, we observed the potential for an even flatter $\mid \tilde \epsilon (\omega )\mid$ dependence within the same band. For instance, this takes place at $f_1=0.23, f_2=0.21, f_3=0.2, f_4=0.19$, and $f_5=0.17$, although these outcomes are not presented here.

Both the spectra of $\tilde \epsilon (\omega )$ and $\mu (\omega )$ can exhibit oscillations within the targeted band at frequencies where $\Re (\epsilon _i)=0$. However, these oscillations can be effectively mitigated by judiciously selecting the values of $n$ for a given $\delta$, or $\delta$ for a given $n$ (see Figs. 2, 3, and 4). Notably, in Fig. 2, oscillations in the spectra of $\mid \tilde \epsilon (\lambda )\mid$ are nearly imperceptible at $\delta =0.099$ even with $n=7$. It is obvious that the magnitude of these oscillations cannot be too high as $\mid \epsilon _i\mid$ is bounded from below by a nonzero $\Im (\epsilon _i)$.

As shown in Fig. 4, $\Re (\mu )$ can be maintained almost constant within the targeted band when increasing the number of the layers within the unit cell $n$. Meanwhile, $\Im (\mu )$ displays two peaks of opposite signs at the band edges. This phenomenon gives rise to the high-frequency and low-frequency peaks in the $\mid \mu \mid$ spectra, clearly evident in Fig. 5. Close to the band center, the contribution of $\Im (\mu )$ to the $\mid \mu \mid$ spectra is marginal. By examining Eqs. (1) and (2), it can be deduced that $n_2 \propto \frac {\Re (\tilde \chi ^{(3)})}{\mid \tilde \epsilon \mid }$ and $\kappa _2 \propto -\sqrt {\frac {\mid \tilde \epsilon \mid -\tilde \epsilon '}{\mid \tilde \epsilon \mid +\tilde \epsilon '}}\frac {\Re (\tilde \chi ^{(3)})} {\mid \tilde \epsilon \mid }$ in this range.

Furthermore, several additional insights can be gleaned from Fig. 5. First, the peak positions in the $\mid \mu \mid$ spectra remain relatively unaffected by $d$. As we found out, this applies to the peaks of $\Re (\mu )$ and $\Im (\mu )$ as well. In addition, as observed in Fig. 7, these positions show minimal sensitivity to the arrangements of the layers. Second, nonlocality tends to decrease the magnitude of the high-frequency peak while amplifying the low-frequency peak of the third-order nonlinear susceptibility. In general, this duality in the impact of nonlocality is unsurprising, considering that nonlocal corrections to the effective permittivity can be of different signs [50]. Third, the magnitudes of both peaks are influenced by the layer arrangement. Fourth, the linear layer arrangement (’12345’) results in a $\omega _0$-centered sub-band of a weak $d$-dependence of $\mid \mu \mid$, indicating weak nonlocality impact on $\mid \mu \mid$. We note that asymmetric arrangements can also produce similar sub-bands, though not centered around $\omega _0$.

Of course, generally the equality $\vert \mu (\omega )\vert =\vert \mu _0(\omega )\vert$ does not imply that the influence of nonlocality on nonlinearity enhancement can be disregarded at this point. This is because the real parts of $\mu (\omega )$ and $\mu _0(\omega )$, as well as their imaginary parts remain distinct.

As it becomes clear from Fig. 6, the enhancement of nonlinearity is closely related to the enhancement of the (mesoscopic) electric field. In its turn, the electric field can be enhanced at frequencies where $\Re (\epsilon _i)\approx 0$ (the ENZ range). In particular, at the edges of the targeted band, $\Re (\epsilon _1)\approx 0$ at $\lambda =1555$ nm, while $\Re (\epsilon _5)\approx 0$ at $\lambda =1455$ nm, that provides a 100-nm bandwidth. Between the edges, the ENZ condition holds for the layers ’2’, ’3’, and ’4’, that provides both the enhanced electric field and effective nonlinear susceptibility at specific frequencies within the operational frequency band.

A notable group of arrangements, ’12345’, ’12354’, ’12453’, ’12543’, ’13245’, and ’13254’, stands out in Fig. 7. These arrangements exhibit a reduced nonlinearity enhancement at the high-frequency band edge and a moderate enhancement (being in close proximity) at the low-frequency edge. At the band center, however, these arrangements display a considerable spread in ${\mid \mu \mid }$ values. An identifiable trait shared by these arrangements is the inclusion of the ’45’ layer combination. In our modeling, these layers possess the highest carrier density, thereby providing the most negative permittivity at the low-frequency band edge.

Further examination of Fig. 7 reveals that different arrangements can have varying numbers of crossing points, where $\vert \mu (\omega,d)\vert =\vert \mu _0(\omega )\vert$. Notably, for the ’12435’ arrangement, there is only one such point ($\lambda \approx 1515.5$ nm), whereas for the ’13245’ arrangement there are five such points ($\lambda \approx$ 1487 nm, 1499 nm, 1508 nm, 1535 nm, and 1541.5 nm). Generally, as $n\to \infty$, the number of crossing points can become arbitrarily large.

It is follows from Fig. 8 that sub-bands of weak nonlocality can arise inside the targeted band. This is particularly characteristic of symmetric arrangements, for which $\Im (\mu )\approx 0$ at the band center.

Dealing with the inverse problem of nonlinearity enhancement, we have briefly considered only the particular case of its high and near constant real part to demonstrate the capability of our technique. Anyway, it seems important to address some key points on this matter.

Upon examining Eqs. (5) and (6), it becomes evident that the enhancement of the effective third-order nonlinear susceptibility can occur either through (i) a high $\mid \tilde \epsilon \mid$ or (ii) a low $\mid \epsilon _i\mid$ within a specified bandwidth. It is worth noting that conditions (i) and (ii) are not mutually consistent. While we can set the real parts of all $\epsilon _i$ to zero within the band, this results in a relatively low $\mid \tilde \epsilon \mid$, as depicted in Fig. 2. Additionally, the values of $\mid \epsilon _i\mid$ cannot be extremely low due to the presence of loss.

In essence, efficient broadband nonlinearity enhancement involves a trade-off between a high $\mid \tilde \epsilon \mid$ and a low $\mid \epsilon _i\mid$. This trade-off is evident in the optimization results. In our specific example, at $d=700$ nm only 3 out of 5 zeros of $\Re (\epsilon _i)$ for each solution fall within the band [$\lambda _{min}, \lambda _{max}$], yielding the best outcome (see Fig. 10, left panel). At the same time, at $d=1000$ nm the best outcome has been obtained for the solution characterized by 2 zeros of $\Re (\epsilon _i)$ within the above band. Furthermore, in our example, we arbitrarily set $\mu _0=3.5$. It is crucial to recognize that this choice was arbitrary, but it is essential to keep in mind that a larger $\mu _0$ leads to a larger minimum of the objective function. This implies an increase in the magnitude of oscillations of $\Re (\mu )$ around $\mu _0$.

6. Conclusions

In summary, we have introduced a phenomenological approach that allows for the straightforward computation of the third-order nonlinear susceptibility within multiphase one-dimensional metamaterials, taking into account the influential role of nonlocal effects. Our investigation of these 1d multilayer nanostructures, utilizing realistic material parameters, has substantiated their inherent capacity to enhance nonlinear responses and expand bandwidths. Notably, nonlocality tends to amplify the nonlinear response at the low-frequency band edge while mitigating it at the high-frequency edge. Furthermore, specific layer arrangements give rise to sub-bands where the influence of nonlocality on nonlinearity is diminished. Dealing with the inverse problem of nonlinearity design, we have made sure that nonlocality results in the appearance of additional solutions which need to be properly addressed.

Throughout this study, numerous model parameters were held constant. Exploring a broader parameter space presents a wealth of opportunities for devising innovative nonlinear nanostructures, endowed with enhanced functionalities and performance characteristics. Particularly compelling is the potential to engineer diverse profiles of $\mu (\omega )$ within targeted frequency bands through judicious parameter tuning.

In this light, we regard this work as an additional step towards bridging the gap between theoretical models elucidating the nonlinear response of subwavelength nanostructures and the tangible realm of plasmonic materials, with a distinct emphasis on the expansion of bandwidths. As we hope, the presented results could help researchers and engineers to tailor materials and nanostructures for desired outcomes in various fields, including nonlinear optics, telecommunications, and materials science.