Self-organization of frozen light in near-zero-index media with cubic nonlinearity

We theoretically demonstrate the existence of frozen light in near-zero-index media with cubic nonlinearity. Light is stopped to a standstill owing to the divergent wavelength and the vanishing group velocity, effectively rendering, through nonlinearity, a positive-epsilon trapping cavity carved in an otherwise slightly-negative-epsilon medium. By numerically solving Maxwell's equations, we find a soliton-like family of still azimuthal doughnuts, which we further study through an adiabatic perturbative theory that describes soliton evaporation in lossy media or condensation in actively pumped materials. Our results suggest applications in optical data processing and storage, quantum optical memories, and soliton-based lasers without cavities. Additionally, near-zero-index conditions can also be found in the interplanetary medium and in the atmosphere, where we provide an alternative explanation to the rare phenomenon of ball-lightning.

. Besides, they can be artificially realized as waveguides close to modal cutoff 30 , using surface phonon polaritons in GaAs quantum wells 31 , or by engineering subwavelength metallic nanowires, nano-spheres, or nano-circuits embedded in dielectric matrices. The latter strategy has enabled the development of epsilon-near-zero (ENZ) metamaterials, which have been investigated for applications such as enhanced transmission 32 , cloaking 33 , energy squeezing in narrow channels 34 , and subwavelength imaging 35,36 . The ENZ regime is inevitably associated with high dispersion and is therefore accompanied by absorption, which can be suppressed by embedding externally pumped active inclusions in the NZI medium 37 .
Here we theoretically investigate self-organization of light in NZI media with Kerr-like instantaneous nonlinearity. In particular, we reveal the existence of fully confined doughnut-shaped solitons with vanishing Poynting vector and angular momentum. In practice nonlinearity enables digging a three-dimensional cavity for light, which in turn remains frozen and self-trapped. We study the effect of loss on stationary light doughnuts by developing a fully numerical soliton perturbative theory, finding that they evaporate over time due to absorption: their amplitude decreases, their frequency blueshifts slightly, and their radius increases. Conversely, if externally pumped active inclusions with inversion of population are embedded within the NZI medium, the opposite scenario takes place and azimuthal doughnuts condensate over time. These findings demonstrate the possibility to freeze light beams in ENZ media, with potential applications in optical data processing and storage, quantum

Results and Discussion
Still light. We consider a generic NZI medium with Drude temporal response and instantaneous Kerr-like nonlinearity (see Methods). Both of these ingredients ensue from free-particle temporal dynamics, which is characteristic of plasmas, metals, transparent conductors, and ENZ metamaterials, all examples of NZI media. In particular, Kerr-like nonlinearity naturally arises from the ponderomotive force in plasmas and metals 41 .
In the linear limit, homogeneous transverse electromagnetic (TEM) waves are solutions of Maxwell's equations with a complex electric field given by  = is the frequency-dependent dielectric constant, which is given by the Fourier transform of the Drude temporal response function  τ ( ) (see Methods). The material dispersion basically depends on two constants: the plasma frequency ω p and the damping rate γ. The linear dispersion relation of TEM waves ω ( ) k is depicted in Fig. 1 in the lossless limit γ → 0, together with the phase and group velocities ω ω ω . Note the cutoff of TEM waves at the plasma frequency ω ω = p , where the medium enters the ENZ regime, the phase velocity diverges, and the group velocity vanishes 22-24 . Homogeneous nonlinear modes. Owing to the vanishing group velocity, nonlinear effects are dramatically enhanced in the ENZ regime 25,26 . For ω ω < p , homogeneous modes with vanishing wave-number, infinite phase velocity, and zero group velocity can be found by neglecting damping and setting , and χ 3 being the material's Kerr coefficient (see Methods). The resulting dispersion relation is plotted in Fig. 2(a). We find zero-index homogeneous modes to have a cutoff at the plasma frequency ω p , where the electric field amplitude drops to zero. In order to evaluate the stability of homogeneous modes, we perturb them with small-amplitude waves: are the perturbation amplitudes with wave-vector q and temporal growth eigenvalue α. Inserting this expression in the Maxwell's equations and retaining only the lowest-order terms in δE 1 and δE 2 , we find a homogeneous system of linear equations, whose non-trivial solutions are signaled by the vanishing of the secular determinant [see supplementary information (SI) for more details on the technical aspects of the theory]. This condition determines the complex temporal eigenvalues α. Instabilities are then associated with positive real parts of the eigenvalue α, indicating unbound amplification of the perturbation. We plot results of the stability analysis in Fig. 2(b-d), and in particular, we depict the maximum of the real part of the eigenvalue, α ′ M . In analogy to standard modulation instability in 1D paraxial systems 3 , the gain spectrum of the perturbations is non-vanishing within a finite wave-vector window and is peaked at a characteristic wave-vector modulus. However, in contrast to 1D paraxial systems, the gain spectrum is 3D and has a non-trivial dependence on polar and azimuthal angles (θ,φ) of the perturbation wave-vector q.
Still azimuthal doughnuts. The modulation instability scenario strongly suggests the presence of still 3D solitons in NZI media. In order to verify this hypothesis, we transform Maxwell's equations into spherical coordinates and search for azimuthally-polarized solutions: . As Maxwell's equations are invariant under a constant phase shift (see Methods), without any loss of generality we can assume that the electric field envelope is real , meaning that we are seeking non-propagating solutions which are not accompanied by a phase flow. Indeed, assuming that such solutions exist, we show that the Poynting vector vanishes thoroughly (see SI). Besides, we seek localized soliton-like solutions vanishing at → ∞ r and at = r 0, θ π = , 0 owing to the azimuthal polarization. Upon examination of the asymptotical expansion of Maxwell's equations for → ∞ r , we find that 3D soliton-like azimuthal solutions can actually exist only in the ENZ regime (see SI). Thus, we discretize derivatives with respect to the radius r n and the polar angle θ m and then transform the differential wave equation for the electric field into a nonlinear algebraic system for the electric field amplitudes φ, , E n m in the two-dimensional grid θ , r n m (see SI). We solve this nonlinear algebraic system by means of an iterative Newton-Raphson algorithm, and find a family of still azimuthal doughnuts [see Fig. 3(a)] for ω ω < p , which presents a cutoff at ω p , where the soliton loses localization and its amplitude vanishes. The frequency-dependent maximum amplitude and the corresponding radius of the still doughnut family are plotted in Fig. 3(b), while a r-θ contour-plot of the squared electric field profile ( )/  Fig. 3(c). The total dielectric permittivity profile Fig. 3(d). Importantly, in the soliton existence domain ω ω < p , the linear dielectric constant is negative  ω ( ) < 0, and thus, at long radius where the electric field amplitude is small, the NZI medium is metal-like. Conversely, in the volume around the radius r max for which the electric field is maximum, nonlinearity is non-negligible and the total dielectric permittivity is positive . From here we see that the existence of still azimuthal doughnuts originates in the extraordinary ability of nonlinearity to dig a dielectric-like 3D cavity within a metal-like environment. This scenario is unique of NZI media, which prevent propagation of the fields outside the induced-dielectric trapping cavity. We emphasize that modulation instability enables the excitation of non-propagating solitons starting directly from unstable homogeneous waves with frequency falling in the ENZ regime.
Doughnut evaporation/condensation. In standard transparent media, the main quantity accounting for optical propagation is the Poynting vector, representing the temporal rate of energy transfer per unit area. For our trapped solitons, the Poynting vector is thoroughly vanishing (see SI), so we describe doughnut self-trapping through the optical-cycle-averaged density of electromagnetic energy. Now, if absorption is taken into account, the energy density is expected to be damped and vanish exponentially over time. A numerical verification of this hypothesis could consist in temporally evolving Maxwell's equation with the doughnut initial condition. However, temporal evolution requires nonlinear 3D finite-difference-time-domain (FDTD) numerical simulations, which are computationally demanding. Besides, traditional approaches used in dielectric and plasmonic waveguides [42][43][44] relying on the slowly-varying-envelope approximation (SVEA) can not be used, as the SVEA does not hold in the ENZ regime 25 . Instead, we have developed a soliton perturbation theory (see SI) capable of accounting for both damping and amplification (e.g., in systems containing externally pumped active inclusions within the NZI medium) under the assumption that (i) damping γ ( > ) 0 or (ii) gain γ ( < ) 0 are much smaller than the soliton angular frequency ω. We further assume that the temporal evolution of the still doughnut adiabatically follows the soliton family, finding that the soliton amplitude (i) decays or (ii) increases over time following the exponential law 2 , where E 0 is the initial field amplitude and γ is a phenomenological absorption/pumping rate. Accordingly, the doughnut (i) expands and blueshifts or (ii) shrinks and redshifts in either case (see SI). The time-dependent field amplitude (blue left y-axis) and doughnut radius (red right y-axis) are plotted in Fig. 4(a) for a representative example, along with three snap-shots of the iso-surface | ( )/ | = × φ − E E r 1 10 S 2 3 at different times in Fig. 4(b-d), where we have assumed as initial condition the doughnut of Fig. 3(a,i) damping γ ( > ) 0 (For the full temporal evolution see movie in the SI). The doughnut evaporates over time, as its amplitude decreases and its radius increases. The gain scenario (ii) γ ( < ) 0 can be interpreted by inverting the temporal direction, so that the doughnut condensates over time, as its amplitude increases and its radius decreases. Ball-lightnings? BLs are rare lightning events with hitherto unknown theoretical explanation [38][39][40] . BLs emit broadband radiation and can either propagate or stand still. Initially considered as myth, BLs have puzzled scientists for centuries and their existence has been questioned until the first recent experiment able to measure their spectrum 40 . Understanding of the nature of BLs is still unsatisfactory as they can not be easily reproduced in laboratory. Among the several theories trying to explain their nature, the so-called maser-caviton theory 38 suggests that BLs are localized high-field solitons forming a cavity surrounded by plasma. Indeed, during thunderstorms, atmosphere can get ionized and become a NZI medium with a plasma frequency falling in the terahertz-microwave spectral region, where rotational levels of water can be excited. The ensuing emitted radiation is thought to remain self-trapped and heat up the air, thus emitting broadband blackbody radiation 38 . This theory explains several aspects of BLs, e.g., their typical size and their motion due to plasma density perturbations, but it does not provide any quantitative description of the self-induced soliton cavity. Following our rigorous calculations, as suggested by the maser-caviton theory, we speculate that BLs may actually ensue from a self-organization process in the ENZ regime, where we theoretically demonstrate the existence of still doughnut solitons, as discussed above. The actual spherical shape of BLs observed in experiments 40 may be due to mixed polarization, heating and higher order nonlinear effects, or the intrinsically incoherent nature of radiation emitted in the atmosphere. The ENZ condition would explain the infrequency of the phenomenon and provides an insightful signature for experimental investigations.

Conclusions
Our investigation of self-organization phenomena in NZI media with cubic nonlinearity has resulted in the demonstration that zero-index nonlinear waves are unstable in all spatial directions and that still azimuthally polarized self-trapped doughnuts can be excited. We have discussed the existence domain of this 3D soliton family with thoroughly vanishing Poynting vector and provided details on its characteristics. Besides, we have studied the effect of loss/amplification, finding that still light doughnuts evaporate/condensate over time, respectively. Our model applies to any NZI medium with cubic nonlinearity and our results are universal as they are rescaled to the relevant physical quantities (plasma frequency ω p , plasma wave-vector k p , Kerr coefficient χ ) 3 of any specific medium in this regime (e.g., metals, transparent conductors, plasmas, and metal-dielectric ENZ metamaterials). Our findings pave the way for the development of novel applications in optical data processing and storage, the realization of quantum optical memories, and the design of soliton-based lasers without cavities. Incidentally, NZI conditions can be found also in the interplanetary medium and in the atmosphere, and we have discussed possible relationships between our results and ball-lightning formation.

Methods
Model. In our investigations we have considered a generic NZI medium with Drude temporal response and instantaneous Kerr-like nonlinearity. Both of these ingredients ensue from free-particle temporal dynamics, which  of the time-evolving doughnut with initial condition as in Fig. 3(a) for where  0 is the vacuum permittivity, χ 3 is the nonlinear susceptibility of the medium, τ δ τ ω γ is the Drude temporal response function, δ τ ( ) is the Dirac delta-function, ω p is the plasma frequency, and γ is the temporal damping rate due to inelastic collisions. Optical propagation is governed by the wave equation where µ 0 is the vacuum permeability.
Still homogeneous waves. Homogeneous nonlinear modes with vanishing group velocity and diverging phase velocity have been calculated by inserting the Ansatz ( , ) = ( ) ω − t e r E r i t in Eq. (2), which in turn enables the calculation of the nonlinear dispersion.
In order to evaluate the stability of homogeneous modes, we have perturbed them with small-amplitude waves where δE 1 and δE 2 are the perturbation amplitudes with wave-vector q and temporal growth eigenvalue α. We have thus numerically calculated the complex eigenvalues α (which real parts represent the instability growth rates) of the ensuing homogeneous system of algebraic equations (see SI).

Still azimuthal doughnuts.
Given the isotropic nature of the system, the most natural coordinates to calculate 3D solitons are spherical θ φ ( , , ) r , where r is the modulus of the position vector, and θ ϕ , are its polar and azimuthal angles (see SI). In our calculations we have assumed that the electric field does not depend on the azimuthal angle, and thus it is polarized along the azimuthal direction φ = φÊ E . Still non-paraxial soliton-like solutions of Maxwell's equations have been calculated numerically by transforming the continuous variables θ , r into a discrete two-dimensional grid θ θ = ∆ , = ∆ r n r m n m with steps θ ∆ , ∆ r , and the azimuthal electric field θ ( , ) φ E r into an ordered vector θ ( , ) = ( , ) φ n m E r V n m . Approximating derivatives by finite differences, Eq. (2) becomes a nonlinear system of algebraic equations, which we have numerically solved through the Newton-Raphson method.