Numerical simulation of light structures in bulk ENZ media with Kerr nonlinearity

A simplified mathematical model is suggested to describe the dynamics of a quasi-monochromatic optical wave in the bulk of an effectively isotropic metamaterial with averaged dielectrical permittivity near zero (ENZ medium), in the presence of a weak spatial nonuniformity, Kerr nonlinearity as well as linear gain due to external pumping. The model is a vector Ginzburg-Landau equation of the general kind, with the dominating curl-curl term in the dispersive operator, and it resembles the equation for electromagnetic waves in plasma [E. A. Kuznetsov, 1974]. In the case of purely real Kerr coefficients, a split-step Fourier method is appropriate for numerical simulations. It makes possible to observe various variants of nontrivial evolution of both central-symmetric and toroidal vector wave structures trapped by a quadratic potential well, as well as nonlinear interaction between the longitudinal and transverse waves in the case of their combination.

💡 Research Summary

The paper presents a comprehensive theoretical and numerical study of light propagation in bulk epsilon‑near‑zero (ENZ) metamaterials that are effectively isotropic, weakly spatially non‑uniform, and exhibit Kerr‑type third‑order nonlinearity together with linear gain supplied by external pumping. Starting from Maxwell’s equations, the authors derive a vector Ginzburg‑Landau equation in which the dominant dispersive operator is the curl‑curl term, reflecting the fact that in ENZ media the phase velocity diverges and the rotational component of the electric field governs the dynamics. The equation includes a quadratic trapping potential U(r)=w r², real Kerr coefficients α (self‑phase modulation) and β (polarization coupling), linear loss γ compensated by a gain term G, and small dissipative coefficients ν⊥, ν∥ that model viscous‑like attenuation of transverse and longitudinal components, respectively.

After appropriate nondimensionalisation (time scaled by Q²τ, length by Qc/ω₀, field amplitude by 1/(Q√α′) with Q≈50) and neglect of the transverse dispersion D⊥, the final governing equation reads

i(E_t – G E) = (1 – iν⊥) rot rot E – (D∥ – iν∥)∇(∇·E) + U(r)E – (1 + iσ)|E|²E – η(E·E)E*

where η = β/α′ and σ = Im α/α′ are taken to be zero for the simulations, i.e., purely real Kerr response.

The authors adopt a split‑step Fourier method of second‑order accuracy in time. Each full time step Δt is split into: (i) a half‑step local nonlinear evolution given analytically by Eq. (7), (ii) M‑1 full‑step linear‑dispersive updates in Fourier space using the exact solution of the linear part (Eq. (9)), and (iii) a final half‑step nonlinear update. The algorithm is implemented on a periodic cubic domain of side 2π with a 192³ grid; modes with |k|>50 are filtered to avoid aliasing. Typical runs require several days of CPU time on a modern workstation and reach evolution times up to t≈300 in dimensionless units.

Three families of initial conditions are explored: (a) centrally symmetric “worm” solutions of the form E∥∝r exp