Field theory description of the non-perturbative optical nonlinearity of epsilon-near-zero media

In this paper we introduce a fully non-perturbative approach for the description of the optical nonlinearity of epsilon-near-zero (ENZ) media. In particular, based on the rigorous Feynman path integral method, we develop a dressed Lagrangian field theory for light-matter interactions and discuss its application to dispersive Kerr-like media with order-of-unity light-induced refractive index variations. Specifically, considering the relevant case of Indium Tin Oxide (ITO) nonlinearities, we address the novel regime of non-perturbative refractive index variations in ENZ media and establish that it follows naturally from a scalar field theory with a Born-Infeld (BI) Lagrangian. Moreover, we developed a predctive model that includes the intrinsic saturation effects originating from the light-induced modification of the Drude terms in the linear dispersion of ITO materials. Our results extend the Huttner-Barnett-Bechler electrodynamics model to the case of non-perturbative optical Kerr-like media providing an intrinsically nonlinear, field-theoretic framework for understanding the exceptional nonlinearity of ITO materials beyond traditional perturbation theory.


I. INTRODUCTION
Epsilon-near-zero (ENZ) low-index optical materials exhibit a number of fascinating optical properties including a dramatic enhancement of nonlinear processes such as harmonic generation, frequency mixing and conversion, electro-optical modulation, and light-induced refractive index changes [1][2][3][4][5][6][7][8][9][10] .These effects are boosted in sub-wavelength nanostructures across their ENZ spectral regions and provide unprecedented opportunities for classical and quantum information technologies, optical spectroscopy, and sensing.Specifically, indium tin oxide (ITO) degenerate semiconductors have recently attracted a significant interest due to their extreme optical nonlinearity resulting in order-of-unity light-induced refractive index changes with sub-picosecond response time 5,11 .A number of impressive ENZ-driven optical nonlinear phenomena and devices have been demonstrated based on ITO nanolayers, including enhanced second-harmonic 12 and third-harmonic frequency generation 10 , ultrafast all-optical modulation 11 , highefficiency optical time reversal 9 , as well as time-varying and spatio-temporal photonic devices 7,8,[13][14][15] .Moreover, the linear optical dispersion properties of ENZ materials and structures based on ITO can be tuned widely with a zero permittivity wavelength λ ENZ (i.e., the wavelength at which the real part of the permittivity vanishes) that extends from the near-infrared to mid-infrared spectral range 16,17 , enabling the engineering of efficient nonlinear phenomena across multiple wavelengths.It was also established that a significant enhancement of the nonlinear response of ITO nanolay-ers can be achieved through control of the material microstructure, specifically by exploiting secondary phases with anisotropic texturing that are largely influenced by the fabrication conditions 18 .Recently, Reshef et al. 19 observed that the conventional expression of the intensity-dependent refractive index loses its meaning for nonlinear ITO materials in the ENZ regime because the refractive index cannot be expanded in a perturbation series.In particular, they found that the contributions of susceptibility terms up to seventh-order to refraction are much bigger than the linear term and, consequently, the nonlinear response of ENZ and low-index media can no longer be interpreted as a small perturbation to linear optics.However, as the authors clearly pointed out, the power series expansion of the nonlinear polarization in ENZ materials does not diverge at the ENZ wavelength.The intrinsically non-perturbative response of ITO thin films has been further enhanced by Shubitidze et al. 20 through the coupling with surface-localized plasmon-polariton Tamm states in nonlinear photonic structures.
In this work, we consider the nonlinear optical response of ENZ media from an alternative theoretical perspective and we develop, based on the rigorous Feynman path integral approach, a novel framework for understanding non-perturbative nonlinear refraction effects, such as the ones recently demonstrated in ITO low-index media, based on nonlinear field theory.In particular, our theory shows how non-perturbative refractive responses emerge naturally in dispersive Kerr-like media from an intrinsically nonlinear field theory characterized by the Born-Infeld (BI) Lagrangian.While focused for simplicity on a scalar case, our approach introduces a novel framework, more general than traditional nonlinear optics, for gaining a more fundamental understanding of the exceptional nonlinearity of ENZ materials beyond perturbation theory.Our work is organised as follows: in section II we briefly introduce the non-perturbative nonlinear refraction of ENZ media within a phenomenological model based on the concept of "field-corrected" permittivity and illustrate the significant impact of low-index on the nonlinear response of ITO thin-films in the ENZ regime.In section III we give a brief introduction to the concept and physical meaning of Feynman path integrals and contextualise their use in optics.We also introduce the basic model, i.e., the Huttner-Barnett-Bechler Lagrangian for light-matter interaction in dispersive ENZ media at the microscopic level, and derive the effective dressed electromagnetic Lagrangian, which constitutes the starting point of this work.
In section IV we introduce the nonlinearity and discuss the usual perturbative case, typical of nonlinear optics while in section V we propose a new approach, based on the Born-Infeld Lagrangian, that goes beyond the perturbative limit of ENZ media and includes the distinctive saturation effects resulting from the bleaching of the Drude linear dispersion.Finally, in section VI we draw our conclusions and future perspective.

II. NONLINEAR OPTICAL RESPONSE OF ENZ MEDIA
The objective of this section is to review the standard perturbative treatment of the induced nonlinear polarization of a Kerr-like medium with intensity-dependent refractive index and to introduce the necessary modifications required to describe ENZ low-index materials.We assume for simplicity a centrosymmetric material (i.e., scalar susceptibility) and neglect its magnetic response.Under these assumptions, the total polarization induced by monochromatic light can be written as: where E(ω) is the incident electric field, χ (1) is the linear susceptibility of the medium, and χ (3) is its lowest-order nonlinear susceptibility.From the above expression we introduce the "field-corrected" susceptibility 21 : and the "field-corrected" dielectric function: The linear refractive index is n 0 r (ω) and ε (1) r is the linear contribution to the relative permittivity.In conventional nonlinear materials, susceptibility terms beyond the third-order are negligible in magnitude and the light induced refractive index change ∆n = n − n 0 is very small under the usual assumption ∆n << n 0 .Therefore, in these cases the change in refractive index induced by the incident field can be well-approximated as follows: where we made use of the binomial expansion.The expression above leads to the conventional definition of the intensity dependent refractive index: where: is the nonlinear refractive index of the medium.It has been recently established by Reshef et al. 19 that for ENZ low-index media the conventional expression for the intensity-dependent refractive index in Eq. ( 5) is inapplicable as the nonlinear field dependence of the refractive index cannot be expressed in a perturbation series.Instead, the full non-perturbative expression below must be considered in the study of ENZ nonlinear media 19 : where we introduced the complex nonlinear refractive index ñ(ω) ≡ n + iκ to emphasize that the quantities in Eq. ( 7) are generally complex and c j are the degeneracy for each relevant nonlinear process considered in the summation.Specifically, for ITO materials excited at large pump intensity up to ≈ 250 GW/cm 2 , it was shown that nonlinear susceptibility up to seventh-order should be included in Eq. ( 7) because their contributions far exceed the linear refraction 19 .Closed-form expressions for the real (n) and imaginary (κ) parts of the complex refractive index (and hence of the nonlinear absorption coefficient α = 2k 0 κ, where k 0 = 2π/λ is the free space wave number) can be obtained by inserting the appropriate perturbative terms of the field-corrected dielectric function inside the generally valid relations linking the real and imaginary parts of the complex refractive index, as shown in Ref. 20 .In particular, if we consider low excitation intensities up to ≈ 10 GW/cm 2 , where only third-order susceptibility terms are important, we obtain the following expressions: where: , ε 1 = 0) for a fixed value of the imaginary part ε 2 = 0.25 and (b) effect of varying the imagi- nary part of the linear permittivity (i.e., ε 2 ) when ε 1 = 0.In all cases we considered a real third-order susceptibility χ (3) = 5.2 × 10 −17 m 2 /V 2 as measured in Ref. 20 Here (χ 1 , χ 2 ) denote, respectively, the real and imaginary parts of the susceptibility function 22 .This approach can be naturally generalized to any desired order in the perturbation expansion of the susceptibility, which is important when studying ENZ materials at large pump intensities.
In order to clearly illustrate the importance of ENZ low-index materials for nonlinear optics applications, we show in We considered susceptibility values up to χ (7) that are typical of high-quality ITO films at their ENZ transition wavelengths, as measured in Refs. 19,20.We see that the largest refractive index changes are always obtained when ε 1 = 0 and the linear losses (ε 2 ) are minimized.We also con- sidered in Figures 2 the variation We considered a fixed ε 2 = 0.25 and varied the imaginary part of the complex third-order susceptibility χ (3) 2 .In all cases the real part of the third-order susceptibility was fixed to χ (3) = 5.2 × 10 −17 m 2 /V 2 as measured in Ref. 20 ated to χ (3) 2 at the ENZ condition defined by ε 1 = 0.These saturation trends of the nonlinear losses have been recently observed experimentally in highly nonlinear ITO nanolayers excited around the λ ENZ wavelength 19,20 .
Before proceeding further with the main result of our work, it is worth pointing out that while Eq. ( 7) is a valid model for describing nonperturbative effects in low-index media, and that Eqs.(8a) and (8b) are only valid in the low-intensity regime.Therefore, extrapolating information on the highintensity regime from those equations can lead to erroneous predictions.This is due to the fact that Eqs.(8a) and (8b) only includes the third-order susceptibility and completely neglect the higher-order terms in Eq. ( 7).These terms, as it is pointed out in Ref. 19 , however, play a significant role when the peak intensity of the pump pulse becomes very large.We will discuss in detail how to overcome this problem in Sect.V C.

III. THE DRESSED ELECTROMAGNETIC FIELD LAGRANGIAN
In this section, we introduce a Lagrangian approach to lightmatter interaction, originally proposed by Huttner and Barnett in 1992 23 , then extended to the case of dispersive media by Bechler in 1999 24 , and recently generalised to quantum nonlinear optics 25 .The approach we develop in this work follows these references, and consists of the following steps: (i) we define the model through the Huttner-Barnett-Bechler Lagrangian; (ii) we calculate the dressed (effective) electromagnetic action via Feynman path integrals; (iii) we derive the expression for the displacement vector in frequency domain, and finally (vi) we extract the information on the dielectric function.A summary of the details for the calculation of the effective Lagrangian is given in Appendix B, and we refer the interested reader to Refs. 24,25for a more complete and detailed discussion of the model.To begin with, let us consider the action: S = dt d 3 x L HBB , with L HBB = L em + L matter + L res + L lmi being the Huttner-Barnett-Bechler (HBB) Lagrangian 23,24 , where 26 for the free electromagnetic field.Here, L matter accounts for the dynamics of the matter (polarization) field P(x,t), modelled as an harmonic oscillator field with a characteristic resonance frequency ω 0 and a static po- larizability β , and L res refers instead to the reservoir field, representing the collection of all decay processes occurring in the material and modelled as a collection of harmonic oscillators interacting with the matter field.The explicit expressions of the matter and reservoir Lagrangians are given in Appendix A. The last term, L lmi is the light-matter interaction Lagrangian, describing the interaction between the electromagnetic field and the matter field.Without loss of gener-ality, we consider the interaction in the usual electric dipole approximation, i.e., L lmi = −E(x,t) • P(x,t).
The dressed electromagnetic Lagrangian, containing all the information on the optical properties of the medium, can then be calculated from the action above by eliminating the matter degrees of freedom through Feynman path integration 27,28 as follows: where the integration measure Dϕ indicates path integration over all possible configurations of the field ϕ 27,29 .A pictorial representation of path integrals and their physical meaning is presented in Fig. 3.The procedure to calculate the effective dressed electromagnetic Lagrangian is sketched in Appendix B and its final expression, suppressing the x-dependence for convenience, is given as follows: where Γ(t − τ) is the (tensorial) Green's function of the material 24 , i.e., the optical response of the material.

A. Effective Lagrangian in Fourier Space
To derive the dielectric function of the medium, it is useful to represent the effective action in frequency domain by using Fourier transformation to obtain S e f f = (dω/2π) dx Le f f , where now Le , and Γ(ω) is the Fourier transform of Γ(t − τ), whose explicit expression is 24 : Here, λ (ω) is a spectral function that contains information about the properties of the medium and the spectral coupling to the reservoir 24 , i.e., the linear losses of the material.Choosing an explicit expression for λ (ω) conveniently allows for different models of dispersive materials.From the effective Lagrangian, we can then derive the electric displacement vector as 30 from which we define the linear dielectric function as ε(ω . This result provides the possibility for the model to describe different classes of materials by suitably choosing different values of λ (ω).

B. ENZ Materials
The dielectric permittivity of an ENZ material around its λ ENZ is derived in a similar manner to that of a Drude metal, with the inclusion of losses and a background permittivity ε ∞ , due to the residual matter polarization of the ion cores serving as background for the dynamics of the nearly free electrons of the material 31 .This leads to the following expression for the dielectric function of an ENZ material around the λ ENZ : where γ accounts for the losses, and the sign at the denominator depends on the adopted convention.The second term of the dielectric function above can be readily obtained from our model, by setting λ (ω) = (ω 2 0 ± iγω)/ω 2 .To account for the first term, instead, we need to consider a slightly modified form of the Lagrangian that includes the residual polarization P ∞ (x,t) = [ε 0 (ε ∞ − 1)/2]E(x,t) due to the ion cores.To do so, we introduce this term in the light-matter interaction Lagrangian L lmi = −E(x,t) [P(x,t) + P ∞ (x,t)].Since the term P ∞ (x,t) is constant with respect to the matter field P(x,t), we can immediately bring it out of the path integrals and insert it directly in the expression for the effective Lagrangian.Including this modification into the effective Lagrangian, we get the correct form of the dielectric function for an ENZ material around its ENZ transition, i.e., Eq. ( 13).

IV. NONLINEAR DIELECTRIC FUNCTION: PERTURBATIVE CASE
Nonlinear electromagnetic effects can be introduced in the model by choosing an interaction Lagrangian of the form , where the nonlinear couplings c n can be real or complex, and their explicit expression is determined by the particular kind of nonlinearity.Normally, |c n | ≪ 1 and the interaction is perturbative.This, for example, is the case of standard nonlinear optics 22 .Adding this interaction term into the HBB Lagrangian allows us to write the nonlinear (perturbative) dielectric function as: A. Nonlinear Refraction As an example to familiarise with this formalism, let us derive the standard intensity law for standard Kerr materials, i.e., n = n 0 + n 2 I 22 .The interaction Lagrangian for a Kerr medium is 29,32 : where χ (3) is the (scalar) third-order nonlinear susceptibility.
To make calculations easier, from this point on we consider a scalar electromagnetic field (i.e., we neglect phase matching and the effect of polarisation), and assume that the electric field only depends on one spatial coordinate, i.e., E(x,t) → E(x,t).Then, for convenience, we suppress the spatial dependence, since it does not play any role in this calculation.Following the prescriptions of Sect.III, we first represent the interaction Lagrangian in frequency domain (see Appendix C), then, using Eq. ( 14) we get the nonlinear dielectric function as: If we now introduce the time-averaged field intensity as 22 I = 2c n ε 0 |E(Ω)| 2 and define the refractive index as n = ℜ{ ε(ω, E)}, we get, using the fact that |χ (3) | ≪ 1, the wellknow result 22 n ≃ n 0 + n 2 I + O(I 2 ), where n 0 = ℜ{ √ ε lin } is the linear refractive index, and n 2 is the usual nonlinear refractive index 22 , in accordance with the results obtained in Sect.II

V. BORN-INFELD LAGRANGIAN AND NON-PERTURBATIVE NONLINEARITIES
In this section we propose a new approach to describe nonlinear materials with fully non-perturbative optical responses.The approach developed in Sect.II already provides a partial solution to this problem, by considering the full nonperturbative action of the square root in the definition of the refractive index (see, for example, Eq. ( 7)).However, the calculation of the dielectric function is still carried out with a standard, perturbative approach.Going beyond this requires to rethink the way nonlinearities are treated at the field theory level.The approach used above, for example, fails to describe a non-perturbative nonlinearity, since every nonlinear interaction is inserted ad-hoc in a perturbative manner, and even if one would take into account all the infinite orders of nonlinearities, the interaction Lagrangian cannot be easily resummed into a closed-form, non-perturbative expression.
For these reasons, a paradigm shift in the way we think about optical nonlinearities is necessary.What we are looking for is a generalized electromagnetic theory, which still has Maxwell's equations as its dynamical equations of motion, but that at the same time can automatically account for optical nonlinearities.Remarkably, such a theory was originally proposed in 1934 by Born and Infeld that developed a classical electrodynamics approach to solve the issue of the diverging electron self-energy 33,34 .The Lagrangian they proposed to describe such theory is known nowadays as the Born-Infeld (BI) Lagrangian and has the following expression 35,36 : where b is a positive real constant, with the dimensions of an energy density, whose physical meaning is that of an upper limit to the energy carried by the field and, ultimately, its amplitude.The BI Lagrangian is the archetype of nonlinear electrodynamics in classical field theory (see, for example, the discussions in Refs. 35,36), and has recently experienced a resurgence in various areas of theoretical physics, including cosmology 37 , nuclear physics 38 , and string theory, where the BI Lagrangian governs the electrodynamics on D-branes 39 .
A peculiar property of the BI Lagrangian is that its equations of motion are still Maxwell's equations, where now the electric and magnetic fields are intrinsically nonlinear.For this reason, the usual constitutive relations D = ε 0 E + P and H = B/µ 0 − M need to be replaced with their nonlinear counterparts 36,40 Moreover, in the so-called strong limit b → ∞, the BI Lagrangian converges to the usual Maxwell Lagrangian, and the nonlinear interactions naturally arise as higher order terms in the expansion with respect to b.

A. The Non-perturbative Dielectric Function
Motivated by these arguments, we then propose to use a BIlike model to describe the non-perturbative optical nonlinearities of materials.In order to achieve this goal, we replace the Maxwell Lagrangian L em and the interaction Lagrangian L int in the nonlinear HBB model with the following Lagrangian: where a is a parameter that we introduce to account for the fact that we are describing the dynamics of the electromagnetic field inside a medium and not in vacuum.The second term serves to ensure that the linear dielectric function can be written as ε ∞ + Γ(ω)/ε 0 .Notice, moreover, that since in non- linear optics we are typically only interested on the nonlinearities caused by the electric field, the magnetic part of the BI Lagrangian plays no role in this setting, and can be ignored.By doing so, the effective non-perturbative Lagrangian now reads: + dτ E * (t)Γ(t − τ)E(τ).
To calculate the non-perturbative dielectric function, we need to switch to a frequency-domain representation.However, the presence of the square root term in the BI Lagrangian makes this step particularly complicated, as one would require to find a suitable Fourier representation, which is generally not an easy task.One possible way to tackle this problem is to write the BI Lagrangian as a power series in b, represent each term of the series in frequency domain, and then re-sum the series.If we do so, and follow the procedure highlighted in Appendix D, we get the following result for the nonlinear dielectric function:  7) (red, solid line).As it can be seen, although the two models agree at low intensities, for I 0 20 GW/cm 2 , the two models start to deviate significantly.The black, dashed line in panel (a), representing the perturbative approximation ∆ n = n 2 I 0 , has been inserted for comparison.To plot these curves, the following experimentally measured parameters (see Ref. 20 ) have been used: ε 1 = 0, ε 2 = 0.45, and χ (3) = 7.58 × 10 −17 m 2 /V 2 .With these pa- rameters we get n 0 = 0.4743 and n 2 = 0.9517 cm 2 /GW.The BI model has been plotted using where ε lin is the linear permittivity of the material, and I c = 2cn 0 b 2 is the critical intensity.The physical meaning of the critical intensity follows from the parameter b originally introduced by Born and Infeld to set a limit on the maximum electric field when the magnetic flux B = 0.This feature, which is a general characteristic of the BI electrodynamics, sets a maximum intensity scale to the static electric field in order to avoid the divergence of the electron self-energy (see, for example, Refs. 33,35,36for a more detailed discussion of this point).In our context, however, one could interpret I c as an intrinsic free parameter of the BI model that can be used to fit the experimental data of the particular material considered.A caveat, however, is necessary: I c cannot be chosen completely arbi-trarily, but it must be chosen in such a way that 1 − I/I c > 0 is fulfilled in order to avoid nonphysical solutions potentially arising when the square root denominator in Eq. ( 21) becomes imaginary.This is the main result of this work.Using the BI Lagrangian enables a novel, non-perturbative approach to describe the optical response of strongly nonlinear media.The price to pay, however, is that intrinsically perturbative parameters such as the nonlinear susceptibilities loose their meaning in this broader context, but can be conveniently restored in the perturbative limit of the BI theory by setting b → ∞.

B. Nonlinear Response of ENZ Low-index ITO
We now compare the two models described in this work, i.e., the non-perturbative BI model above, and the one derived in Sect.II.To this aim, we use ITO as the nonlinear ENZ material, for which the parameters appearing in Eq. ( 13) have been extracted from fitting the dielectric function model of ITO described in Ref. 20 around λ ENZ and are ε ∞ = 3.99, ω p = 3.16 × 10 15 Hz, and γ = 0.185 × 10 15 Hz.We take as a figure of merit for the comparison the nonlinear refractive index change ∆ n = n(I) − n(0), and the nonlinear absorption change ∆α = α(I) = α(0).We assume an impinging electromagnetic field oscillating at frequency ω = 1.57× 10 15 Hz, corresponding to the ENZ wavelength of ITO, i.e., λ = 1200 nm.
Figure 4 shows the comparison between the nonperturbative approach developed in Sect.II (blue, solid curve), with the BI model of Sect.V A (red, solid curve).To match the latter to the former, we require that they agree in the small intensity regime, where the BI model predicts ε (np) ≃ ε lin + a I/(2I c ) + O(I 2 ) and the usual permittivity model for Kerr nonlinearities is ε ≃ ε lin + 3χ (3) |E| 2 .This leads to identifying a = 4n 0 n 2 I c .With this fixed, I c remains the only free parameter of the BI model.Figure 4(b) shows the nonlinear absorption calculated from the BI model (red, solid line) compared with the results from Eq. ( 7) (blue, solid line) for the case χ (3) 2 = 0, i.e., assuming that the imaginary part of the third-order susceptibility is zero.The condition χ (3) 2 = 0 is a necessary one for the BI model.Contrary to standard nonlinear optics, where the susceptibility is in general a complex quantity, the free parameter of the BI model, namely the critical intensity I c , is instead an intrinsically real parameter.This limits the possibility for the BI, in its actual form, to correctly reproduce the behaviour of Fig. 2(b), when  χ (3) 2 = 0.As it is discussed below, however, including this feature in the BI model might be possible by introducing extra terms in the Lagrangian, giving rise to an effective field theory with a complex free parameter.

C. Nonlinear Saturation Effects
The model proposed in Eq. ( 7) introduces the susceptibility as a scaling factor between the polarization and the inci-  7) (blue,solid line).Both models agree at low intensities (see inset), but while the non-perturbative model predicts an ever increasing refractive index change, the saturated model of Eq. ( 22) correctly reproduces the physical behaviour at high itensities.The parameters used to plot these curves are the same as in Fig. 4. For the saturated permittivity, we have taken ε sat = 4, compatible with the experimen- tal measurement of ITO presented in Ref. 20 the black, dashed line, representing the perturbative approximation ∆n = n 2 I, has been inserted for comparison.
dent field.Although in principle this effective macroscopic approach is legitimate and the effective susceptibility defined with this model can be used to readily compare the effective and microscopic models, one should be careful in doing so.In fact, for this comparison to be valid, one should ensure that the effective description remains within the bounds of the microscopic process and the experimental conditions upon which such a model is based.In fact, using the effective χ (3) model incorrectly, could lead to erroneous results.
In the case of ENZ materials, for example, it is generally accepted that the intensity-dependent refractive index arises from absorptive effects and free electrons 11,[41][42][43][44] This implies the existence of a finite relaxation time for the nonlinear susceptibility or, equivalently, the intensity-dependent refractive index, which makes the values of the nonlinearity significantly dependent upon the duration of the excitation pulse 45 .As a result, the nonlinear response of ENZ materials will eventually saturate to a value corresponding to the optical response of the material in the absence of the free-electron Drude dynamics.This saturation process has been discussed and experimentally verified under different circumstances and for different ENZ materials 11,18,[46][47][48] .
To correctly take into account this saturation effect starting from Eq. ( 7), one possible strategy is to introduce higher order susceptibilities, which progressively become more relevant as the intensity increases.This is the approach taken, for example, in Ref. 19 .To account for this mechanism in our intrinsically nonperturbative model, described by Eq. ( 21), however, we cannot use this argument, since we do not have acces to the different orders of nonlinear susceptibility.Instead, we add a saturation term in our model, resulting in the following version of our non-perturbative model for the permittivity of ITO; where I 0 is the high-power saturation intensity, and its value can be determined by requiring that at high pump intensities the Drude component of the linear permittivity is completely saturated, and only the Lorentz-Tauc contribution remains to describe the ITO permittivity 18,20 .Calling then δ ε = ε sat − ε lin , where ε sat is the saturated permittivity at high intensity, the saturation intensity I 0 can be fixed by the behaviour of the above equation near infinity as From a phenomenological point of view, the saturation of the nonlinearity of a material at high intensity arises from the fact that the nonlinear coupling constant itself acquires an intensity dependence at high values of I. Since in our nonperturbative model the critical intensity I c plays the role of the nonlinear coupling coefficient, one could insert this saturation mechanism by simply replacing a constant value of I c with an intensity-dependent critical intensity of the form I c (I) = I c /[1 − (I c /I 0 )I], so that for I 0 → ∞ one recovers the unsaturated case (i.e., low intensity), and for I 0 < I c one can use the parameter I 0 to compensate the divergence in Eq. ( 21) and induce saturation.
From a more fundamental field theory perspective, on the other hand, the presence of the extra term in Eq. (22) needs to be justified at the Lagrangian level by adding an extra term in the Born-infeld Lagrangian, leading to consider: where the parameter α is related to I 0 .In Born-Infeld nonlinear electrodynamics, a similar extra term is related to the dual electromagnetic tensor 35 , which is ultimately proportional to The presence of this extra term could also be justified by more complicated models of nonlinear electrodynamics, such as Dirca-Born-Infeld-like Lagrangians used to study gauge fields on D-branes 39 , which allow higher-order terms inside the square root in Eq. ( 17) in terms of coupling of the electromagnetic field to background scalar fields.In Fig. 5 we compare the effective susceptibility model of Eq. ( 7) and the saturated Born-Infeld model given by Eq. (22).As it can be seen, while the models agree very well at low intensities, where the nonlinear response is governed mainly by I c (see inset), at high intensity, the saturated model correctly predicts the saturation of the nonlinear index change, Nonlinear Refractive Index Ref. [19]   Eq. ( 22) FIG. 6. Nonlinear refractive index as a function of the peak intensity caluclated using the saturated (red, solid line) Born-Infeld model, compared with the nonperturbative model with 5 th and 7 th order nonlinearities from Ref. 19 (blue, solid line).The parameters used to plot the saturated Born-Infeld model are I c = 100 GW/cm 2 and ε sat = 1.252 + i0.606, as taken from Ref. 19 .The values for the non- linear susceptibilities for the blue curve have been taken from Table 1 in Ref. 19 .
while the unsaturated one erroneously predicts ever increasing ∆ n values.As an additional proof of the validity of our model, in Fig. 6 we compare the predictions of Eq. ( 22) with the non-perturbative, saturated refractive index model developed by Reshef and co-workers in Ref. 19 , where the summation under the square root symbol in Eq. ( 7) includes also 5 th and 7 th order nonlinear susceptibility contributions.As it can be seen, our saturated BI model reproduces the experimental data from Ref. 19 quite well.Finally, it should be noted that the proposed saturation model does not take into account nonlocal and gradient terms in the nonlinear response of ITO as well as surface nonlinearities, and therefore cannot be directly applied to the study of nonlinear nano-structures where larger refractive index modulations have been experimentally demonstrated 2,20,49 .
In this work we have only considered the impact of electric nonlinearities, neglecting the effect of the magnetic field.Taking this also into account might as well introduce new degrees of freedom, that can be used to tame the divergence.In recent years, several other models of nonlinear electrodynamics, such as NED or ModMax 36 , have been raising some interest in the field theory community.Adapting those models to non-perturbative nonlinear optics might also be interesting, as it could possibly lead to alternative perspectives.Lastly, it is worth noticing that when b → 0, the BI Lagrangian becomes zero and the electromagnetic field in this limit has no dynamics 50 .In this case, a Hamiltonian description of the electromagnetic field is necessary and shows the occurrence of topological objects, such as hopfions 51 and knotted fields 52 .Making the connection between the BI approach and nonlinear ENZ materials could therefore open new horizons for the study of electrodynamics in ENZ media, possibly hinting at the existence of hidden topological properties as well.

VI. CONCLUSION
In this paper we have introduced a fully non-pertubative approach for the description of the nonlinear optical refraction properties of ENZ low-index media and non-perturbative metamaterials, which also includes the intrinsic nonlinear saturation effects originating from the light-induced bleaching of the Drude permittivity.In particular, we used the rigorous Feynman path integral method to develop an effective Lagrangian field theory for light-matter interactions and investigated, within a simple scalar model, the relevant case of Indium Tin Oxide (ITO) nonlinear materials.Remarkably, we established that their non-perturbative refractive index response naturally emerges from an intrinsically nonlinear electrodynamics theory characterized by the Born-Infeld (BI) Lagrangian.At large optical intensities, the developed BI model naturally accounts for all orders of nonlinear susceptibilities of materials, unlike the case of traditional perturbative approaches.Moreover, we have shown that very concept of nonlinear susceptibility coefficients loses its meaning in the non-perturbative regime and must be substituted with the nonlinear dielectric function in Eq. ( 22), directly derived from the BI Lagrangian.A closed-form, two-parameter formula for the non-perturbative dielectric function was derived, containing all orders of nonlinearities of an arbitrary material.Our model goes beyond the state-of-the-art perturbative expansion of the refractive index, and introduces a new paradigm for calculating the dielectric function in a fully non-perturbative manner.These results extend the Huttner-Barnett-Bechler electrodynamics model to the case of non-perturbative photonic ENZ metamaterials and devices providing a nonlinear, fieldtheoretic framework for describing exceptional nonlinearities beyond traditional perturbation theory.

DP DR exp
(B1) where S em = dt d 3 x L em and the quantities F 1,2 are functionals of the matter and reservoir fields, defined as: We are interested in calculating F 2 [P, R] explicitly.Intuitively, since L res only depends on R 2 and Ṙ2 [see Eq. (A2)], we can use Gaussian integrals to compute F 2 [P, R].For continuous fields, they are defined as follows: 27,29 Dφ exp − 1 2 φ , Âφ + ψ, φ = exp 1 2 ψ, Â−1 ψ , (B4) where ψ is an auxiliary field (source term), Â is a differential operator (typically well-behaved and invertible), Â−1 is the propagator associated to the field φ , i.e., the inverse of the differential operator Â, and ψ, φ is a shorthand for ψ, φ = d n x ψ * (x)φ (x), (B5) with n being the dimension of the space in which the fields ψ(x) and φ (x) are defined.
To calculate the integral in Eq. (B3) using Gaussian path integrals, we first need to bring its exponent into a quadratic form similar to that of the integral above.To do so, let us first integrate by part the matter-reservoir terms, so that we shift the time derivative from the reservoir field to the matter field as follows: where in the third line we introduced the • notation with d n x = d 3 x dt dτ dω and we defined: to write this term already in a form compatible with Eq. (B4).Notice, moreover, that we have introduced an extra integration with respect to the variable τ (and a correspondent Dirac delta function) for later convenience.
To transform the rest of L res into a quadratic form like the one in Eq. (B4), we first use the identity (∂ t denotes time derivation) and then integrate by parts once more, so that we are left with one term proportional to R∂ 2 t R and another term proportional to ω 2 R 2 .Putting everything together, and introducing the differential operator Â as: we get, for the exponent of Eq. (B3) We can then use Eq.(B4) to calculate the path integral with respect to the reservoir degrees of freedom appearing in Eq. (B3), obtaining: If we now use the definition of B(τ) given above and define the inverse of the operator Â as the Green's function, i.e., Â−1 = G(t − τ, x) 53 , we get: + dt dτ dx P(t) • G(t − τ, x) • P(τ) .(B13) The calculation of the second path integral, to eliminate the matter degrees of freedom, can be done using the same line of reasoning.The only complication is that the expression of the differential operator needed to define a quadratic form in the matter field P(x,t) is more complicated than Â, as it results in an integro-differential operator, whose Green's function is, in general, very complicated to calculate.The reason for this is the presence, in the result above, of the quadratic form P(t) • G(t − τ) • P(τ), which includes the effects of the reservoir in the Green's function for the matter field.Details on how to calculate explicitly this second integral, to arrive to the final result presented in Eq. ( 11), are given in Refs. 24,25e then calculate the derivative of the interaction Lagrangian with respect to the field E * (ω), using the identity 29 to otain: where to get the last equality we used the definition of the time-averaged intensity of the electromagnetic field 22 .

FIG. 1 .
FIG.1.Nonlinear refractive index change as a function of the peak pump intensity.(a) Effect of varying the real part of the linear permittivity away from the ENZ condition (i.e., ε 1 = 0) for a fixed value of the imaginary part ε 2 = 0.25 and (b) effect of varying the imagi- nary part of the linear permittivity (i.e., ε 2 ) when ε 1 = 0.In all cases we considered a real third-order susceptibility χ (3) = 5.2 × 10 −17 m 2 /V 2 as measured in Ref.20 Figures 1 (a-b) the computed nonlinear refractive index change as a function of the peak pump intensity for different values of the real (a) and imaginary (b) parts of the linear permittivity of the medium.

FIG. 2 .
FIG.2.Nonlinear absorption change as a function of the peak pump intensity.(a) Effect of varying the imaginary part of the linear permittivity (i.e., ε 2 ) when ε 1 = 0. (b) We considered a fixed ε 2 = 0.25 and varied the imaginary part of the complex third-order susceptibility χ

FIG. 3 .
FIG. 3. Pictorial representation of path integrals for fields.The final field configuration ϕ b can be reached from the initial one (ϕ a ) by summing over all the possible field configurations (paths in field space, pictorially represented by the landscape).Each path is assigned a weight function exp [i S/h], where S is the classical action of the field.The classical equations of motions of the field, represented here by a solid, blue line, are those minimising the action.The dark green region around it represents the region where the interference of all the contributing paths linking the initial (ϕ a ) and final (ϕ b ) field configurations is maximum.
FIG. 4. (a) Nonlinear refractive index change and (b) nonlinear absorption as a function of the peak intensity as calculated using the Born-Infeld model (blue, solid line) and the non-perturbative refractive index model in Eq. (7) (red, solid line).As it can be seen, although the two models agree at low intensities, for I 0 20 GW/cm 2 , the two models start to deviate significantly.The black, dashed line in panel (a), representing the perturbative approximation ∆ n = n 2 I 0 , has been inserted for comparison.To plot these curves, the following experimentally measured parameters (see Ref.20 ) have been used: ε 1 = 0, ε 2 = 0.45, and χ (3) = 7.58 × 10 −17 m 2 /V 2 .With these pa- rameters we get n 0 = 0.4743 and n 2 = 0.9517 cm 2 /GW.The BI model has been plotted using I c = 300 GW/cm 2 .

FIG. 5 .
FIG.5.Nonlinear refractive index change as a function of the peak intensity calculated using the saturated (red, solid line) Born Infeld model, compared with the nonperturbative refractive-index model in Eq. (7) (blue,solid line).Both models agree at low intensities (see inset), but while the non-perturbative model predicts an ever increasing refractive index change, the saturated model of Eq. (22) correctly reproduces the physical behaviour at high itensities.The parameters used to plot these curves are the same as in Fig.4.For the saturated permittivity, we have taken ε sat = 4, compatible with the experimen- tal measurement of ITO presented in Ref.20 the black, dashed line, representing the perturbative approximation ∆n = n 2 I, has been inserted for comparison.