Structured light entities, chaos and nonlocal maps

TL;DR

This paper introduces nonlocal nonlinear maps derived from Green functions that emulate complex spatiotemporal phenomena like solitons and vortex lattices in optical resonators, linking iterative maps to PDEs of Ginzburg-Landau type.

Contribution

It generalizes classical maps to include convolution integrals, connecting iterative maps with PDEs, and demonstrates their application in modeling optical phenomena.

Findings

Nonlocal maps are equivalent to Ginzburg-Landau PDEs.

Maps can simulate spatial solitons and vortex lattices.

Noise helps stabilize entities and reduce artifacts.

Abstract

Spatial chaos as a phenomenon of ultimate complexity requires the efficient numerical algorithms. For this purpose iterative low-dimensional maps have demonstrated high efficiency. Natural generalization of Feigenbaum and Ikeda maps may include convolution integrals with kernel in a form of Green function of a relevant linear physical system. It is shown that such iterative $nonlocal$ $nonlinear$ $maps$ are equivalent to ubiquitous class of nonlinear partial differential equations of Ginzburg-Landau type. With a Green functions relevant to generic optical resonators these $nonlocal$ $maps$ emulate the basic spatiotemporal phenomena as spatial solitons, vortex eigenmodes breathing via relaxation oscillations mediated by noise, vortex-vortex and vortex-antivortex lattices with periodic location of vortex cores. The smooth multimode noise addition facilitates the selection of stable entities and elimination of numerical artifacts.

Equations

$$\frac{E_{n+1,m}-E_{n,m}}{\Delta t}=f(E_{n,m})E_{n,m}+$$

$$\chi \frac{E_{n,m+1}+E_{n,m-1}-2E_{n,m}}{{\Delta z}^2}.$$

$$E_{n+1}(\vec{r})={\int }_{-\infty }^{\infty }K(\vec{r}-{\vec{r}}^{,})f(E_n({\vec{r}}^{,}))d^3{\vec{r}}^{,}\approx$$

$$\sum _{mx,my,mz}K(\vec{r}-{\vec{r}}_m)f(E_n({\vec{r}}_m))S(mx,my,mz),$$

$$\frac{\partial E(z,\vec{r},t)}{\partial z}+\frac{1}{V}\frac{\partial E}{\partial t}+\chi \frac{{\partial }^{2}E}{\partial t^2}=f(E)E.$$

$$\frac{\partial E(z,\vec{r},t)}{\partial z}+\frac{n}{c}\frac{\partial E}{\partial t}+\frac{i}{2k}{\Delta }_{\perp }E=ik{n}_2|E{|}^{2}E.$$

$$i\hslash \frac{\partial \Psi (z,\vec{r},t)}{\partial t}=-\frac{{\hslash }^2}{2m}{\Delta }_{\perp }\Psi +U(\vec{r},t)\Psi +\frac{4\pi {\hslash }^2{a}_S}{m}|\Psi {|}^2\Psi .$$

$$\frac{\partial E(z,\vec{r},t)}{\partial z}+\frac{n}{c}\frac{\partial E}{\partial t}+\frac{i}{2k}{\Delta }_{\perp }E=ik{n}_2|E{|}^{2}E+$$

$$\sigma N_o(z,\vec{r})E\exp [-2\sigma {\int }_{-\infty }^t|E{|}^{2}d{t}^{{}^{\prime }}]-\gamma E,$$

$$E_{n+1}={\lambda }_F{E}_n(1-E_n),$$

$${\delta }_F=\underset{M\to \infty }{\lim }\frac{{\lambda }_M-{\lambda }_{M-1}}{{\lambda }_{M-1}-{\lambda }_{M-2}}\to 4.6692,$$

$$E_{n+1}=g(E_n)\{1-\tanh [g(E_n)]\},$$

$$E_{n+1}=\frac{g(E_n)}{1+{{\chi }_3}^2{L}^2|g(E_n){|}^2},$$

$$\displaystyle P(I)dI\cong\Bigl{[}1+\cos(\Delta\phi)\Bigr{]}d(\Delta\phi){\>},$$

$$P(I)dI=\frac{d(\Delta \varphi )}{\pi },P(I)=\frac{1}{\pi \sqrt{I(1-I)}}.$$

$$E_{n+1}=RE_{n}\cdot\exp\Bigl{[}ikn_{o}L_{c}+ikn_{2}|E_{n}|^{2}L_{nl}\Bigr{]}+TE_{0},$$

$$\displaystyle E_{n+1}=TE_{0}+RE_{n}\cdot\exp\Bigl{[}ikn_{o}L_{c}+\sigma N_{n}L_{nl}ikn_{2}|E_{n}|^{2}L_{nl}\Bigr{]}$$

$$\frac{N_{n+1}-N_n}{\Delta t}=+\frac{N_0-N_n}{T_1}-\sigma N_n|E_n{|}^2,$$

$$E_{n+1}(\vec{r})=\hat{D}\widehat{Fr}f(E_n(\vec{r})),$$

$$\hat{D}=1+\frac{\xi {\partial }^2}{2\partial t^2},\xi =\Delta z\frac{{\partial }^{2}k}{\partial {\omega }^2},$$

$$\widehat{Fr}=1+\frac{i\Delta z}{2k}{\Delta }_{\perp },$$

$$E_{n+m}(\vec{r})=\hat{D}\hat{Fr}{...}f\Bigl{[}\hat{D}\hat{Fr}f\Bigl{[}\hat{D}\hat{Fr}(f[E_{n}(\vec{r}))]\Bigr{]}\Bigr{]},$$

$$E_{n+m}(\vec{r},t)=\lim_{m\rightarrow\infty}\hat{D}\hat{Fr}{...}f\Bigl{[}\hat{D}\hat{Fr}f\Bigl{[}\hat{D}\hat{Fr}(f[E_{n}(\vec{r}))]\Bigr{]}\Bigr{]},$$

$$E(z+\Delta z,{\vec{r}}_{\perp },t)=\hat{D}\widehat{Fr}f(E(z,{\vec{r}}_{\perp },t)).$$

$$\displaystyle E(z+\Delta z,\vec{r}_{\bot},t)=E+\Delta z\frac{\partial E}{\partial z}=ik{\>}\Delta z{\>}E\Bigl{[}n_{0}+n_{2}|E|^{2}\Bigr{]}{\>}{\>}$$

$$+E\Delta z\frac{{\partial }^{2}k}{\partial {\omega }^2}\frac{{\partial }^{2}E}{\partial t^2}+\frac{i\Delta z{\Delta }_{\perp }E}{k},$$

$$\frac{\partial E(z,\vec{r},t)}{\partial z}+\frac{\partial^{2}k}{\partial\omega^{2}}\frac{\partial^{2}E}{\partial t^{2}}+\frac{i}{2k}\Delta_{\bot}E+ik\Bigl{[}n_{0}+n_{2}|E|^{2}\Bigr{]}E=0.$$

$$E(z+\Delta z,\vec{r}_{\bot})=\Bigl{\{}1+\frac{idz}{2k}\Delta_{\bot}\Bigr{\}}E(z,\vec{r}_{\bot}),$$

$$E(z=0,{\vec{r}}_{\perp })=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\exp [i\vec{\kappa }\cdot {\vec{r}}_{\perp }]\tilde{E}(z=0,\vec{\kappa })d^2\vec{\kappa },$$

$$E(\Delta z,\vec{\kappa})=\Bigl{\{}1-\frac{idz\vec{\kappa}^{2}}{2k}\Bigr{\}}E(0,\vec{\kappa}),$$

$$E(2\Delta z,\vec{\kappa})=\Bigl{\{}1-\frac{idz\vec{\kappa}^{2}}{2k}\Bigr{\}}^{2}E(0,\vec{\kappa}),$$

$$E(m\Delta z,\vec{\kappa})=\Bigl{\{}1-\frac{idz\vec{\kappa}^{2}}{2k}\Bigr{\}}^{m}E(0,\vec{\kappa}).$$

$$E(L,\vec{\kappa})=\lim_{m\rightarrow\infty}\Bigl{\{}1-\frac{iL\vec{\kappa}^{2}}{m2k}\Bigr{\}}^{m}E(0,\vec{\kappa}).$$

Full text

Keywords: Discrete maps, integral transforms, solitons, vortices, switching fronts, periodic solutions, vortex lattices, chaos, turbulence, probability density.

I Introduction

Complete spatial synchronization in nonlinear systems is replaced by turbulent states with a fast decay of correlations. The propagation dynamics ranges from switching waves to long-living localized excitations and spatiotemporal solitons Cross:1993 . Numerical modeling of these complex nonlinear distributed systems is based upon finite-difference schemes which emulate the basic solutions of partial differential equations such as localised solitonic entities Malomed:2006 , propagation fronts, stable phase-locked periodic configurations alike lattices and self-organized vortex clusters. The alternative numerical approach is outlined. It is shown that discrete-time and continuous space dynamical systems composed of sequence of local nonlinear point-to-point Feigenbaum and Ikeda-like maps Kaneko:1992 and nonlocal diffusion-dispersion integral transforms are equivalent in a quite general set of cases to conventional partial differential equations of Ginzburg-Landau type. Such $nonlocal$ $maps$ are itself the robust numerical schemes that exhibit the transitions from purely chaotic states to localization in momentum space inherent to corresponding coherent entities alike spatial solitons and vortex lattices in the presence of spatial noise. In some experimentally achievable cases the statistics generated by iterative maps is shown to have the quantitative similarity with statistics of multimode random fields Okulov:2009 .

This computational approach is relevant to spatially distributed nonlinear systems which demonstrate the noise-mediated origin of complex patterns. Superconductors of type II in external magnetic field demonstrate the hexagonal supercurrent lattices with nodal lines of penetrating magnetic induction $\vec{\bf B}$ Abrikosov:1957 . Thin slice nonlinear optical materials with feedback mirrors generate triangular intensity patterns Firth:1991 . Wide aperture solid state lasers emit coherent multivortex phase-locked beams of rectangular symmetry Chen:2001 ; Staliunas:1995 ; Okulov:2008 . Delayed feedback systems of coupled oscillators perform reservoir computations Gauthier:2015 with Boolean ring networks Gauthier:2017 .

The transitions from regular patterns to turbulent states and vice versa are controlled by a limited number of control parameters namely ambient temperature ${\bf T}$, magnetic induction, density of carriers $n$, gain $G$ and loss $\gamma$, quadratic-qubic-quintic nonlinearities $\chi_{2},\chi_{3},\chi_{5}$, carrier frequency $\omega=ck=2\pi/\lambda$, wavelength $\lambda$, dielectric permittivity $\epsilon$ , magnetic permeability $\mu$ and their spatial distributions. Geometrical control parameters are system dimensions at different scales, aspect ratios (degree of system asymmetry), say ratio of transverse dimensions $d,D$ to longitudinal extension $L$, focal lengths of optical components $F$, frame of reference parameters namely vectors of displacement velocity $\vec{\bf V}$ and reference frame rotation angular velocity $\vec{\bf\Omega_{\oplus}}$ Okulov:2013 . The variations of control parameters force system to switch between ordered structures of different symmetry or toggle system into spatially turbulent states with fast decay of correlations. The regions of stability are defined via calculation of instability increments (Lyapunov exponents) Okulov:2000 .

The above mentioned self-organized structures are produced in framework of master equations of Ginzburg- Landau type. These evolution equations contain first order derivative of the order parameter $E$ over time $t$ and diffusion terms $\Delta_{\bot}=\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+\frac{\partial^{2}}{\partial z^{2}}$ with second order derivatives over spatial coordinates, so the simplest explicit finite difference numerical scheme has the form of iterative mapping, which mix the field values in adjacent points of numerical mesh $z_{m-1},z_{m},z_{m+1}$ separated by interval $\Delta z$ at each time step $\Delta t$ as:

[TABLE]

Behavior of such discrete dynamical system (DDS) ranges from exactly solvable diffusion for purely real $\chi$ and diffraction for purely imaginary $\chi$ towards spatiotemporal instabilities Kaneko:2018 , localized structuresMalomed:2011 and turbulence. The role of nonlinear term $f(E_{n,m})$ in pattern formation is crucial Cross:1993 ; Anishchenko:2018 . In many cases the numerical schemes (I) are replaced for much more sophisticated DDS of implicit type Malomed:2007 or else convolution - like nonlocal integral transforms alike fast Fourier transform (FFT) under a proper choice of spatial filtering Siegman:1975 ; Okulov:1993 ; Okulov:1994 :

[TABLE]

where $S({mx,my,mz})$ is a finite volume element instead of infinitesimal one $d^{3}\vec{r}$, kernel $K(\vec{r}-\vec{r}^{,})$ is a Green function of linear version of (I, I), i.e. response for delta-function $\delta(\vec{r}-\vec{r}^{,})$ in right part of scalar diffusion equation, reads as

$K(\vec{r}-\vec{r}^{,})\approx(\chi/\Delta t)^{-1/2}exp(-\sqrt{\chi}|\vec{r}-\vec{r}^{,}|^{2}/\Delta t)$

or for scalar parabolic diffraction equation this reads as

$K(\vec{r}-\vec{r}^{,})\approx(\chi/\Delta t)^{-1/2}exp(-i\sqrt{\chi}|\vec{r}-\vec{r}^{,}|^{2}/\Delta t)$ Okulov:1988 .

Such discrete numerical approach proved to be extremely effective for modeling of nonlinear dynamics inherent to evolution partial differential equations of so-called parabolic type Konotop:1994 . Apart from nonlocal dispersive laplacian there exists a nonlinear term $f(E_{n})$ which depends on square modulus of field $E$ as in Kolmogorov-Petrovskii-Piskounov (KPP) equation Okulov:1988 :

[TABLE]

Here the diffusion term $\chi=F^{2}/k^{2}D^{2}$ may be responsible for spatial filtering of high transverse harmonics in laser cavity by intracavity diaphragm and iteratively repeated nonlinear transformation of the scalar field $E$ in ring laser with intracavity second harmonic generation or Raman scattering Okulov:1986 . More realistic models are based upon nonlinear Shrodinger equation, known also as Ginzburg-Landau equation (NLS - GLE) which captures the interplay of phase-amplitude modulation during propagation of complex field which is the source of modulational instability Konotop:1994 , solitons Malomed:2006 and collapse Cross:1993 . For propagation of light pulse in Kerr dielectric NLS - GLE reads as:

[TABLE]

The similar dynamics is embedded in Gross-Pitaevskii equation for macroscopic wavefunction $\Psi$ of Bose-Einstein condensates Pitaevskii:1999 :

[TABLE]

where $U(\vec{r},t)$ is confining potential of arbitrary complexity Lembessis:2018 , $a_{S}$ is scattering length, $m$ is mass of boson. In a presence of gain $G$ and losses $\gamma$ the NLS-GLE may have a form of NLS with Frantz-Nodvik resonant gain term Okulov_self:1988 , relevant to amplification with stimulated emission cross-section $\sigma$ of light pulse of duration $T_{2}<\tau<T_{1}$ in a rare earth doped dielectric with transverse $T_{2}$ and longitudinal relaxation times $T_{1}$:

[TABLE]

We will discuss discrete iterative maps (I) for numerical modeling of chaotic and regular spatiotemporal propagation inherent to equations (3, 4, 5, I). The article is organized as follows: section II is devoted to overview of regular and chaotic iterative dynamics of real $1\bf D$ Okulov:1984 and complex point maps Ikeda:1982 , section III outlines the direct link between evolution PDE s (3-4) and maps (1-2), section IV contains examples of localized solutions obtained with nonlocal maps being equivalent to equations (3-4) , the spatially periodic and chaotic lattices are obtained numerically in section V and nonlinear dynamics in the presence of multimode fluctuations is presented in section VI with conclusive discussion in section VII.

II Iterative maps with universal behavior.

Feigenbaum demonstrated the universality in iterations of real mapping with parabolic maximum of the unit interval into itself alike logistic map Feigenbaum:1978 :

[TABLE]

where the sole control parameter $\lambda_{F}$ completely determines the discrete time evolution of this simplest dynamical system. In particular he shown that two universal irrational numbers, namely $\delta_{F}=4.6692..$ scales the separation of the values of $\lambda_{F}=\lambda_{1},\lambda_{2},\lambda_{3}...\lambda_{M-1}.\lambda_{M}...$ where period-doubling bifurcations occur(fig.1):

[TABLE]

and $\alpha_{F}=2.5...$ scales the location of limit cycle points in phase-space.

There exists a variety of nonlinear optical systems whose dynamics might be approximated by iterates of maps with parabolic maxima. The intuitively attractive example is a toy model of ring laser with intracavity nonlinear losses Okulov:1984 (fig.2a). Radiation of electric field amplitude $E_{n}$ circulates between confining mirrors along closed trajectory and it passes repeatedly through gain element, diaphragms and nonlinear elements. The successive passages of field through fast amplifying medium with gain $E_{n+1}=g(E_{n})\rightarrow GE_{n}{\>}at{\>}small{\>}E_{n}$ and nonlinear medium with quadratic $\chi_{2}$ or cubic $\chi_{3}$ susceptibilities are described by following maps:

[TABLE]

[TABLE]

Both maps have parabolic maxima and their bifurcation points are condensed to different values $G=\lambda_{chaos}$ (fig.1) of laser system gain with the same universal speed $\delta_{F}=4.6692...$. At these bifurcation points the deterministic dynamical system generates chaotic time series. Thus in a model of single transversal mode solid-state laser with nonlinear losses Okulov:1984 this dynamical regime corresponds to generation of spatially coherent but temporally chaotic radiation.

Apart from universality the above $1D$ dynamical systems might be considered as deterministic source of random numbers. In contrast to logistic maps whose range of chaotic oscillation amplitudes is limited within interval $\lambda\in[0,1)$, the time series generated by chaotic optical cavities are produced by mapping of semi-infinite interval on itself Okulov:1988 ; Okulov:1986 ; Okulov:1984 .

Randomization of time series is so strong at chaotic accumulation points $\lambda_{chaos}$ that probability density functions ($\bf PDF$) for amplitudes $E_{n}$ at a certain values of bifurcation parameter $\lambda,G$ are very close to experimentally obtained histograms for interference of statistically independent Stockes pulses reflected from independent phase-conjugating Brillouin mirrors Basov_Lett:1980 . Noteworthy phases of these Stockes pulses are random because stimulated Brillouin scattering originates from thermal acoustic fluctuations. As a result the recorded interference pattern $I\cong 1+\cos(\Delta\phi)/2$ of the two beams with phase difference distributed uniformly at interval $\Delta\phi\in[0,\pi]$ is also random though light intensity $I$ has exact theoretical probability density $P(I)$ Okulov:1983 (fig.3c):

[TABLE]

Both histograms $\bf{(b,c)}$ at fig.3 are perfectly fitted with $P(I)=1/[\pi\sqrt{I(1-I)}]$ exact probability densities. For this particular case there exists a remarkable coincidence of dynamical chaos $\bf PDF$ and interference pattern $\bf PDF$.

In a more realistic models the nonlinear self-phase modulation inherent to the semiconductor lasers the Kerr cubic nonlinearity [16] may be directly introduced into point maps (9). For this purpose Ikeda constructed a map for complex envelopes of electric field of ultrashort pulses circulating in single transverse mode ring cavity Ikeda:1982 (fig.2b):

[TABLE]

where $E_{0}$ is external optical pump, $R$ is reflectivity of entrance mirror, $T$ is transmission of entrance mirror, $n_{0}$ is linear index of refraction, $n_{2}$ is Kerr component of nonlinear index of refraction, $L_{nl}$ is width of Kerr slice, $L_{c}$ is length of cavity. For more general model of laser cavity with nonstationary gain and population inversion lifetime $T_{1}$ the generalized Ikeda map had been introduced in Hollinger:1985 and subsequently generalized for wide area laser Okulov:2008 :

[TABLE]

where $N_{n}$ is population inversion at a given bounce of pulse $E_{n}$ in cavity, $\sigma$ is stimulated emission cross-section, $N_{0}$ is the pump rate of amplifying medium, $\Delta t={L_{c}n}/c$ is discrete time step of map.

III Iterations of discrete maps and parabolic

partial differential equations.

The spatiotemporal evolution of pulse envelops $E_{n}$ described by equations (9,10,12,13) may be represented by sufficiently long operator products as is shown in Okulov:1988 . Indeed each equation (9,10,12,13) in medium of length $L$ admits the decomposition for sequence of thin slices of linear dispersive elements and nonlinear nondispersive elements. One may expect that in the limit of infinite number $n\rightarrow\infty$ of an infinitely thin slices with $\Delta L=L/n$ such an artificial medium will be equivalent to perfect continuous medium. Within each slice the propagation of pulse $E_{n}$ is exactly integrable so that the following product of maps (2) is evident for passage through one slice:

[TABLE]

where

[TABLE]

is dispersion operator,

[TABLE]

is diffraction operator. In this picture the propagation of pulse $E_{n}(\vec{r},t)$ through $m$ slices is modeled as a product of operators:

[TABLE]

for continuous time variable $t=m\Delta t$ this product becomes:

[TABLE]

Consider infinitesimal slice $\Delta z=L/n$ and use the map (13) for calculation of pulse envelope after passage through it:

[TABLE]

Substitution of operators $\hat{D}$ and $\hat{Fr}$ in this product gives:

[TABLE]

where the second identity leads immediately to NLS-GLE equation:

[TABLE]

On the other hand the infinite chain of operators may be used for construction of exact solution of linear Shrodinger equation (NLS with $n_{2}=0$) at finite propagation distance $L=m\Delta z$. For this purpose consider pulse propagation in free space using the following map:

[TABLE]

Let us decompose envelope $E$ in Fourier integral:

[TABLE]

The substitution in (22) leads to Fourier components after first iterate, i.e. after propagation distance $\Delta z$:

[TABLE]

The next iterates are represented as follows:

[TABLE]

[TABLE]

In order to obtain solution at finite distance $L$ let us use limit $m\rightarrow\infty$ under apparent constraint $m\Delta z=L$ :

[TABLE]

After rearrangement of this formula we have:

[TABLE]

Next because ${\lim_{{\xi\rightarrow\infty}}}(1-\frac{1}{\xi})^{\xi}=e$ is known as Euler number, we have for Fourier components after propagation at finite distance $L$:

[TABLE]

next after return to coordinate space we have:

[TABLE]

after substitution of initial Fourier spectrum:

[TABLE]

we obtain exact solution as Fresnel-Kirchoff integral Okulov:1988 :

[TABLE]

IV Localized wavetrains as fixed points of nonlocal maps

Bifurcation diagrams of discrete maps show the location of fixed points (fig.1a). One may use constructive analogy between fixed points of finite-dimensional maps and stationary self-similar solutions of evolution PDE s Moloney:1988 . Original idea was to construct eigenfunctions of nonlinear resonators with Kerr medium from soliton solutions of NLS-GLE using boundary conditions Moloney:1983 . The complications of this technique arise from the fact, that exact soliton solutions of NLS-GLE are asymptotic objects on the whole propagation axis $z$. As a result such a generalization of conventional theory of solitons to finite space interval bounded by cavity mirrors is not trivial. The alternative approach is to use Fox-Lee method when diffraction is taken into account by calculation of Fresnel-Kirchoff integrals at each round-trip Okulov:1988 . For the simplest Fabry-Perot resonator the mapping of field at $n$-th passage into field at $n+1$-th passage is as follows:

[TABLE]

where $D(\vec{r})$ may be a smooth, say hypergaussian or even step-like Heaviside function $\theta(d-|\vec{r}|)$ , where $d$ is diaphragm width, $L$ is distance between mirrors. This product of convolution integral operator $\hat{Fr}$ and point map $f(E_{n}(\vec{r}))$ had been named as nonlocal nonlinear map Okulov:1988 .

The basic properties of nonlocal maps are visible clearly when different spatial scales are taken into account. Here the key parameter is Fresnel number $N_{f}=kD^{2}/L$ of resonator. For the simplest plane-parallel Fabry-Perot cavity $N_{f}$ is the number of Fresnel zones on a given mirror $M_{1}$ visible from opposite one $M_{2}$ Born_Wolf:1972 . There are two limits each have clear physical meaning. The first one is the limit corresponding to geometrical optics i.e. $\lambda\rightarrow\infty$ . In this case evolution of spatial structure $E_{n}(\vec{r})$ follows to point transformations Okulov:1986 . The other case is a single spatial mode limit $N_{f}\sim 1$ when spatial filtering during each passage through the cavity is strong enough so $E_{n}(\vec{r})$ has a predefined single mode shape, say $\bf{TEM}_{oo}$ whose amplitude evolve in time as (8-9) Okulov:1984 or as (11-12) Okulov:2008 in a presence of self-phase modulation. The most interesting case is the intermediate one $N_{f}>1$ when mode interactions are mediated by nonlinearities, diffraction and dispersion.

Consider first the formation of solitary waves as a result of phase-locking of transverse modes. Such a localized wavetrain is expected to be the eigenfunction of nonlocal map $E_{n+1}(\vec{r})=\hat{Fr}f[E_{n}(\vec{r})]=\Upsilon E_{n}(\vec{r})$ , where $\Upsilon$ is eigenvalue. To get exact solutions for localized transverse structures one may consider a specially configured nonlinear cavities. In particular the confocal cavities (fig.2c) Okulov:1988 ; Okulov:2000 which are known to have the set of degenerate eigenmodes with identical frequencies greatly facilitates the phase-locking of eigenmodes and formation of stable nonlinear localized wavetrains. Noteworthy the excitation of threshold solitons due to saturable absorption and gain Okulov:2000 ; Okulov:1999 had been realized experimentally in confocal resonator(fig.2c) Taranenko:1998 .

The other possibility is the excitation of thresholdless spatial solitons due to gain saturation or nonlinear parametric processes alike second harmonic generation Okulov:1988 . The cavity is again the confocal Fabry-Perot resonator with two spherical mirrors $\bf M_{1},\bf M_{2}$ both have the same focal length $F$.

For detailed numerical modeling it is worthwhile consider the saturable gain medium is described by complex maps similar to Ikeda equations Ikeda:1982 :

[TABLE]

The saturable absorbing medium described by analogous set of Ikeda-like maps:

[TABLE]

where $L_{am},L_{ab}<<L_{r}$ are the thicknesses of amplifying medium and absorber media correspondingly, relevant to experimental situation, $\sigma_{am},\sigma_{ab}$ are the stimulated emission cross-sections, ${T_{1}}_{am},{T_{1}}_{ab}$ are longitudinal relaxation times for amplifier and absorber correspondingly placed near opposite confocal mirrors or deposited upon their surfaces. The propagation of fields $E_{n}$, $\tilde{E}_{n}$ between mirrors is described by already defined Fresnel-Kirchoff integral nonlocal maps (IV) with parabolic phase-modulation $\exp[-i2k|\vec{r}|^{2}/F]$ induced by mirrors:

[TABLE]

The standard numerical evaluation with FFT to equilibrium stationary eigenmodes for dicrete time step $\delta t=2L_{r}n/c$ is achieved for a different levels of accuracy via 10 - 150 iterates Okulov:1994 . The $2D$ spatial solitons (fig.4) were obtained for the cavity (fig.2c). The key assumptions for generation of spatial solitons were threshold excitation guaranteed via smallness in linear regime of gain $\sigma_{amp}N_{n}L_{amp}$ compared to absorption $\sigma_{ab}\tilde{N}_{n}L_{ab}$, faster saturation of absorption compared to gain, Kerr self-focusing in gain slice and higher saturable self-defocusing in absorber slice. Under above restrictions the threshold excitations in the form of bright spatial soliton Malomed:2006 emerge after a light $\delta$-like click in a given section of computational mesh.

In order to get analytical solutions the following decomposition in Taylor series to the third order proved to be useful:

[TABLE]

where $G=\exp[\sigma_{am}N_{0}L_{am}]$. The exact analytical solution $E(x)=E_{n+1}(x)=E_{n}(x)$ in absorber plane (fig.2c) in $\bf 1D$ case has the following form:

[TABLE]

This solution is generalization of the conventional Rayleigh formula Born_Wolf:1972 for focal spot size known as $\Delta x\cong\lambda F/d$:

[TABLE]

The effective width is not standard Rayleigh one $W_{Ral}\cong\lambda F/d$. Indeed in our case the dissipative soliton width is modified in the presence of gain and losses as: $W_{sol1}\cong\lambda F\sqrt{G}/(d\sqrt{1-\alpha G})$. As a result the effective width of spatial soliton diverges when gain approaches to lasing threshold $G\rightarrow\alpha^{-1}$. The stationary solution $\tilde{E}(x)$ in the amplifier plane (fig.2c) is Fourier transform of $E(x)$ :

[TABLE]

Noteworthy both solutions in Fourier conjugated planes are hyperbolic secants Okulov:2000 as it follows from exact result for the $sech$ spectrum Ryzhik:2014 . The inherent solitonic relation Moloney:1988 between amplitude and soliton width known as ”area theorem” is embedded here in explicit form.

Linear stability analysis with respect to small amplitude spatial harmonics had been realized in Okulov:2000 following the standard perturbative technique Okulov_self:1988 with linearized master equation (IV):

[TABLE]

where $\zeta$ is instability increment, $\psi_{\zeta}$ are linear excitation modes in effective potential produced by solitons (38), (IV). The spectrum of these infinitesimal excitations $\psi_{\zeta}$ consists of the two sets Okulov:2000 . Noteworthy the above eigenvalue problem is isomorphic to quantum mechanical problem of scattering a particle of mass $m$ on $sech^{2}(\nu x)$ potential well Landau:1977 :

[TABLE]

The first set of excitations consists of unbounded running plane waves $\psi_{\zeta}\cong\exp(\zeta t+ipx)$ with continuum spectrum $\zeta$. All $\zeta$ in this set are negative, thus all unbounded plane wave excitations $\psi_{\zeta}$ quench exponentially in time. The other set of bounded excitations $\psi_{\zeta}$, contains both negative and positive energies $\zeta$ being equal to positive and negative instability increments correspondingly:

[TABLE]

One may get exact formulas for boundaries of zero $\zeta$ which separate areas of stable increments from unstable ones as is shown at fig.5. The vital consequence from this linear stability analysis is a necessity of specially arranged filtering of spatial harmonics $\psi_{\zeta}$, belonging to the bounded set of excitations with negative spectrum of energies. Filtering of these positive instability increments $\zeta$ will stop the growth of excitations and the stability of solitons (38), (IV) is guaranteed in this case.

The other confocal cavity configuration is ring Sagnac-like cavity with thin-slice nonlinear gain medium and spatial filter in Fourier conjugated planes Okulov:1988 where exact solutions for localized solitonic excitations do exist. The exact solution for spatial solitons had been obtained by searching eigenfunctions of similar nonlocal maps which also include Fox-Lee convolution integral:

[TABLE]

The stability analysis in this case is almost identical to the previous one summarized above in eqs. (38, IV, 43, 44) and fig.5. In both cavities with gain and losses considered above the $1D$ spatial solitons have generic link of width and amplitudes in accordance with area theorem. The interesting common feature of solitonic widths in both cases:

[TABLE]

is that their widths $W_{sol1}$ and $W_{sol2}$ are basically the generalization of Rayleigh formula for the width of the focal spot of a thin parabolic lens $W_{Ral}$ illuminated by a plane monochromatic wave of wavelength $\lambda$ and aperture $d$ Born_Wolf:1972 .

V Periodic transverse structures in

laser cavities obtained with nonlocal map numerical schemes

Nonlocal maps (IV, IV, IV) contain a rich self-organization dynamics. Apart from spatial solitons whose shape fit the exact solutions (38, IV) the eigenmodes of wide aperture optical resonators with intermediate values of Fresnel number $1<N_{f}<\infty$ demonstrate a various nonlinear phase-locking regimes from spatially periodic lattices to fully chaotic states which are speckle fields characterized by randomly spaced optical vortices collocated with zeros of complex field amplitudes $E_{n}(\vec{r})$ Okulov:2009 .

The simplest case of periodic structure formation is possible in a plane-parallel Fabry-Perot cavity with a thin slice gain medium having periodic gain distribution in transverse directions $G(\vec{r})=G(\vec{r}+\vec{p})$ , where $\vec{p}=\vec{x}p_{x}+\vec{y}p_{y}$ , and $\vec{x},\vec{y}$, are unit orts in Cartesian coordinates. The self-imaging or Talbot effect is inherent to periodic field distributions. Because Fox-Lee nonlocal maps (22-24) derived above are exact solutions of Maxwell- Bloch equations in paraxial approximation direct substitution of spatially periodic field $E_{n}(\vec{r})=E_{n}(\vec{r}+\vec{p})$ in (22- 24) immediately proves the Talbot identity of so-called self-imaging of spatially periodic fields at propagation distances $z_{T}=2mp^{2}/\lambda$, m is integer in a set of cases of commensurability of $p_{x},p_{y}$. The corrections due to finite asymmetric aperture having widths $d_{x},d_{y}$ are given by:

[TABLE]

The generalization of Talbot theorem which is exact result of conventional diffraction theory to nonlinear resonators is as follows Okulov:1990 . It had been shown that for the thin gain slice of arbitrary optical nonlinearity the spatially periodic field is self-imaged from one mirror upon another provided the distance between mirrors is a multiple of a Talbot one $L_{c}=mz_{T}=m2p^{2}/\lambda$ , $m$ is positive integer. The issue of stability of such periodic configurations requires a special attention. Nevertheless the numerical investigations with Fox-Lee nonlocal maps (IV, IV, IV) implemented with the aid of standard fast Fourier transform (FFT) routines upon [128x128] and [1024x1024] computational meshes have shown the stable eigenmodes composed of 5x5 and 8x8 phase-locked periodically spaced filaments Okulov:1993 ; Okulov:1994 with field distributions $E_{n+1}(\vec{r})\cong E_{n}(\vec{r})$ close to in-phase wavefronts and almost single-lobe far field Fourier spectra with suppressed side lobes.

The self-organized vortex lattices in laser output became known since 2001 Chen:2001 for microchip laser oscillators composed from diode-pumped thin Nd:YAG gain slice in stable Fabry-Perot cavity with long focus $F$ output mirror. The qualitative agreement with predicted periodic vortex lattices Staliunas:1995 had been found. The separation of longitudinal modes $\Delta\omega_{c}=\pi c/(nL)$ was large enough to facilitate the interaction of a dense set of transverse modes. Fresnel number was in the range $10^{2}<N_{f}<10^{3}$. The nonlocal map approach gave straitforward numerical technique for modeling of the phase-locked vortex lattices (fig.6). For accurate numerical simulation the computational mesh of a moderate resolution [512,512] proved to be sufficient to get the stationary vortex-antivortex state. The detailed pattern of phase dislocations, wavefronts and effective fields of velocities of such lattices Okulov:2004 demonstrated the vortex pairing (fig.6). The effective magnetic fields $\vec{\tilde{B}}$ which may be realized with such vortex configurations Lembessis:2016 ; Lembessis:2017 associated with the each phase-locked vortex are counter-directed for adjacent vortices. This happens in contrast to conventional Abrikosov vortex lattices Abrikosov:1957 where all vortices produce co-rotating currents around magnetic field lines.

Numerical results on melting of this vortex-antivortex lattices Okulov:2003 proved to be in a close similarity with Berezinskii- Kosterlitz- Thauless scenario V.L.Berezinskii:1971 ; Kosterlitz:1972 . The melting of vortex-antivortex lattices ignited by increase of optical gain $\sigma_{am}N_{0}L_{am}$ was envisioned as unbounding the vortices and subsequent loss of the long-range order. The observed dynamics and proliferation of vortices to random locations fits with conventional model of speckle fields generated by superposition of plane waves with random phases and directions of propagation Okulov:2009 .

VI Modeling of thermal noise mediated patterns

The inclusion of noise in dynamical equations (IV, IV, IV) for multimode laser dynamics is natural, because each computational building block of nonlocal maps has a clear physical meaning Anishchenko:2018 . In particular the spontaneous emission in cavity means the addition to $E_{n}(\vec{r})$ the multimode speckle field composed via superposition of randomly directed plane waves Okulov:1991 :

[TABLE]

where $A_{\delta}$ is normalization constant, $k_{n_{x}}^{2}+k_{n_{y}}^{2}+k_{z}^{2}=k^{2}$ . The noisy additions to population inversion $\delta N_{n}(\vec{r})$ and density of carriers $\delta\tilde{N}_{n}(\vec{r})$ are generated via the identical procedure. With these noise sources the master equations are transformed to somewhat more complicated one without decrease of computational speed:

[TABLE]

The equations for thin gain slice are also modified straightforwardly as follows:

[TABLE]

The same holds for the thin saturable absorber slice:

[TABLE]

For low Fresnel number $5<N_{f}<40$ solid state laser cavity the off-axis alignment of optical pump may lead to emission of stable topologically charged vortices (fig.7) Chen:2018 ; Chen:2018optexp ; Suslov:2013 . Such output patterns are highly desirable for metrological Dowling:2012 ; Dowling:2016 and secure free space applications Zeilinger:2015 . The temporal spectrum of such laser oscillator is broadened by a set of factors including relaxation oscillations Chen:2001 . The fig.8 demonstrates the typical behavior of laser output in the presence of noise. The stepwise switching on the population inversion leads to relaxation oscillation Letokhov_Suchkov:1967 with characteristic time scale $\tau_{rel}\cong\sqrt{\tau_{c}T_{1_{am}}}$ because of different time scales of $T_{1_{am}}$ and photon lifetime in a cavity $\tau_{c}$(fig.8a). As expected in a steady-state regime the power spectrum of microchip laser oscillator has a well visible relaxation peak due to perpetual disturbances of intracavity field $\delta E_{n}(\vec{r})$, population inversion $\delta N_{n}(\vec{r})$ and density of carriers in absorber $\delta\tilde{N}_{n}(\vec{r})$ (fig.8b).

VII Conclusions

The toy models of $\bf 1D$ map for unidirectional ring laser and $\bf 2D$ Ikeda map for standing wave lasers are shown to be easily modified into much more realistic models via simple Fresnel-Kirchoff convolution integral transform to mediate interaction of spatially distributed point maps. The exact and numerical solutions for spatial solitons were obtained and their stability analysis had been performed. For high Fresnel number the iterations of nonlocal maps have shown the fast convergence rate to stable square vortex lattices known from table-top experiments Chen:2001 . The inclusion of noise leads to realistic relaxation oscillations power spectra during sufficiently short iterations time intervals.

Diversity of dynamical regimes of self-organization in spatially distributed nonlinear systems and computationally fast generation of stable spatial structures via nonlocal maps had been demonstrated in this work. The close links with conventional evolution equations Malomed:2007 ; Konotop:1994 ; Letokhov_Suchkov:1967 had been established. The stability issues remain still a subject of a very complicated analytical studies, as this requires the accurate calculation of eigenvalues and eigenfunctions of perturbations of solitons, vortices and their lattices. The easiness of numerical implementation of nonlocal maps especially with the aid of FFT routines promise new results in this field. In particular the clear physical meaning of the Fox-Lee method Fox_Lee:1966 , the natural inclusion of boundary conditions in numerical and analytical schemes and well elaborated numerical procedures for each computational block of nonlocal maps gives the firm guaranties for avoiding numerical artifacts, better convergence rates and accurate comparison with exact results and alternative numerical approaches.