Optical spatial shock waves in nonlocal nonlinear media

ABSTRACT

Dispersive shock waves are fascinating phenomena occurring when nonlinearity overwhelms linear effects, such as dispersion and diffraction. Many features of shock waves are still under investigation, as the interplay with noninstantaneity in temporal pulses transmission and nonlocality in spatial beams propagation. Despite the rich and vast literature on nonlinear waves in optical Kerr media, spatial dispersive shock waves in nonlocal materials deserve further attention for their unconventional properties. Indeed, they have been investigated in colloidal matter, chemical physics and biophotonics, for sensing and control of extreme phenomena. Here we review the last developed theoretical models and recent optical experiments on spatial dispersive shock waves in nonlocal media. Moreover, we discuss observations in novel versatile materials relevant for soft matter and biology.

Graphical Abstract

KEYWORDS:

• Nonlocal nonlinear media
• thermal nonlinearity
• dispersive shock wave
• nonlocal Kerr effect
• nonlocal nonlinear Schroedinger equation

PACS:

• 42.65.-k
• 42.65.Hw
• 42.65.Sf
• 78.20.N-

1. Introduction

Dispersive shock waves (DSWs) are rapidly oscillating solutions of hyperbolic partial differential equations that contrast the generation of multivalued regions through the formation of undular bores [1–12]. This class of phenomena was investigated in several physical fields, initially in shallow water waves [13] and ion-acoustic waves [14], then in oceanography [15], pulses propagation in photonic fibers [16,17], Bose-Einstein condensates [18–23], quantum liquids [24], photorefractive media [25], plasma physics [26], viscous fluids [27], and diffracting optical beams [5,28–43].

In 1967 Gardner, Greene, Kruskal, and Miura developed a method to solve the Korteweg-de Vries (KdV) equation, called inverse scattering transform (IST) [44]. Among all the equations solvable by IST, which allowed to find the mathematical formulation of exact solutions of such nonlinear models, KdV and the nonlinear Schrödinger equation (NLSE) belong to the case with dispersive regularization of the aforementioned-multivalued singularity. NLSE is a universal model that describes many phenomena, in particular, a third-order nonlinear phenomenon in optics: the Kerr effect [45], a refractive index perturbation linearly scaling with the light intensity. Kerr effect can be generalized to the nonlocal case when the nonlinear response in a specific point depends on entire beam transverse profile. This occurs, e.g., in thermal media [4,11,29,34,46–60]. In these materials, light propagation is affected by a highly nonlocal Kerr nonlinearity, ruled by nonlocal NLSE.

Unfortunately, IST is still of little use for the nonlocal NLSE and other theoretical approaches must be conceived, despite some recent progress in 2D media [61,62]. For many years, Whitham modulation and hydrodynamic approximation have predominated in solving the nonlocal NLSE [2,29]. However, hydrodynamic approximation cannot describe light propagating beyond the shock point, and two new techniques have been developed, one coming from nuclear physics, the time asymmetric quantum mechanics (TAQM) [12,63–73], which models the nonlinear wave intrinsically irreversible propagation as a superposition of decaying resonances, and the wave turbulence theory [33,39,41,74–77], which uses a statistical interpretation of nonlinear optics.

This review aims to summarize all the current theoretical models to describe wave breaking of nonlocal NLSE solutions in diffracting optical beam propagation, and to highlight some of the most recent experimental observations of DSWs in spatial nonlinear photonics.

After an introductory section about the derivation of nonlocal NLSE in Section 2, we report the main theoretical approaches and results related to DSWs. Section 2.1 explains in detail the difference between the wave-breaking due to local Kerr effect, which causes shock both in phase and in intensity, and the one in nonlocal Kerr media, where the beam intensity follows the phase singularity adiabatically [29]. The most recent theoretical models of nonlinear wave propagation in highly nonlocal nonlinear media are treated in Sections 2.2, 2.3. Section 2.2 treats DSWs generated by laser beams and gives an analytical description of their intrinsic irreversibility, due to the complexity of the dynamics rather than losses [12]. Section 2.3 illustrates the importance of nonlocality in random dispersive waves nonlinear interaction to produce giant collective incoherent shock waves [39,41,76].

The second part of the manuscript is a collection of experiments on DSW generation in thermal media. Section 3 reports observations in Rhodamine solutions [29]. Output beam intensity profiles in Section 3.1 are modeled by TAQM both in two-dimensional experiments, where decaying states describe the longitudinal propagation [37,38], and in the one-dimensional approximation, having the proof that TAQM is an excellent approach also to analyze transverse intensity profiles beyond the shock point [40]. The interplay of nonlinearity and disorder is illustrated in Section 3.2. There, observations in Rhodamine with silica spheres [32] and in silica aerogel [35] exhibit the competition between randomness and nonlocal Kerr effect. DSW generation processes in chemical [36] and biological solutions [43] are illustrated in Section 3.3. Last but not least, Section 3.4 shows recent experiments on the transition from coherent shocklets to a giant incoherent DSW in a photon fluid, modeled by wave turbulence theory [39].

2. The nonlocal nonlinear Schrödinger equation

From Maxwell’s equations, considering a region with zero charge, current and magnetization, we obtain the following electric field wave equation:

(1) $-{\mathrm{\nabla }}^2\mathbf{E}+\frac{1}{c^2}{\mathrm{\partial }}_t^2\mathbf{E}=-\frac{1}{{ϵ}_0{c}^2}{\mathrm{\partial }}_t^2\mathbf{P},$(1)

with $\mathbf{E}$ the electric field and $\mathbf{P}$ the medium nonlinear polarization [78].

The relation between $\mathbf{P}$ and $\mathbf{E}$ depends on the material properties. Including all the nonlinear terms, we have:

(2) $\mathbf{P}={ϵ}_0({\chi }^{(1)}\mathbf{E}+{\chi }^{(2)}\mathbf{E}\mathbf{E}+{\chi }^{(3)}\mathbf{E}\mathbf{E}\mathbf{E}+\dots )={\mathbf{P}}^{(L)}+{\mathbf{P}}^{(NL)},$(2)

where $1+{\chi }^{(1)}=n_0^2$, $n_0$ is the medium refractive index, ${\chi }^{(2)}$ and ${\chi }^{(3)}$ are tensors denoted as second and third-order susceptibility, respectively.

One must take into account the temporal delay between the instant when the electric field reaches the medium and the medium response. For this reason, this radiation-matter interaction is more properly represented by the following non-instantaneous superposition of linear and nonlinear polarization [78]:

(3) $\mathbf{P}={ϵ}_0({\chi }^{(1)}\ast \mathbf{E}(t)+{\chi }^{(2)}\ast \mathbf{E}\mathbf{E}(t)+{\chi }^{(3)}\ast \mathbf{E}\mathbf{E}\mathbf{E}(t)+\dots )={\mathbf{P}}^{(L)}+{\mathbf{P}}^{(NL)},$(3)

where $\ast$ is the convolution product:

$\begin{array}{r}{\chi }^{(n)}\ast \mathbf{E}\dots \mathbf{E}(t)={\int }_{-\mathrm{\infty }}^t\mathrm{d}{t}_1{\int }_{-\mathrm{\infty }}^t\mathrm{d}{t}_2\dots {\int }_{-\mathrm{\infty }}^t\mathrm{d}{t}_n{\chi }^{(n)}(t-t_1,\dots )\mathbf{E}(\mathbf{R},t_1)\dots \mathbf{E}(\mathbf{R},t_n).\end{array}$

If we have a third-order isotropic and centrosymmetric material, the nonlinear polarization is:

(4) ${\mathbf{P}}^{(NL)}(\mathbf{R},t)={ϵ}_0{\int }_{-\mathrm{\infty }}^t\mathrm{d}{t}_1{\int }_{-\mathrm{\infty }}^t\mathrm{d}{t}_2{\int }_{-\mathrm{\infty }}^t\mathrm{d}{t}_3{\chi }^{(3)}(t-t_1,t-t_2,t-t_3)\mathbf{E}(\mathbf{R},t_1)\mathbf{E}(\mathbf{R},t_2)\mathbf{E}(\mathbf{R},t_3)$(4)

and the related dielectric tensor changes as:

(5) ${ϵ}_{new}=ϵ+{ϵ}_2⟨\mathbf{E}\cdot \mathbf{E}⟩,$(5)

where $⟨\mathbf{E}\cdot \mathbf{E}⟩=\frac{1}{2}|\mathbf{E}{|}^2$ is the square of the electric field time average. The final refractive index causes the Kerr effect [78], a phenomenon that consists of a perturbation of the medium refractive index, proportional to the field intensity:

(6) $n=\sqrt{{ϵ}_{new}}=\sqrt{ϵ+{ϵ}_2⟨\mathbf{E}\cdot \mathbf{E}⟩}\approx n_0+n_{2}I,$(6)

with $I=|\mathbf{E}{|}^2$ the field intensity and $n_2$ the Kerr coefficient.

The nonlocal Kerr effect is a third-order phenomenon, but the radiation-matter interaction depends on the whole intensity profile, as occurs in thermal media. In these materials, when an optical beam propagates, it locally heats the medium, and the resulting temperature gradient generates a variation of the density distribution and a refractive index perturbation [29,53,59]:

(7) $\mathrm{\Delta }n={\left(\frac{\mathrm{\partial }n}{\mathrm{\partial }T}\right)}_0\mathrm{\Delta }T,$(7)

with ${\left(\frac{\mathrm{\partial }n}{\mathrm{\partial }T}\right)}_0$ the thermo-optic coefficient of the sample at the steady-state. It turns out that the nonlinear response induced at a specific spatial point is carried away to the surrounding region, and the size of this extended region determines the range of nonlocality. The heat conduction in optical thermal materials was termed ‘response with an infinite range of nonlocality’ [51] until 2007, when A. Minovich et al. [53] proved experimentally and theoretically that the nonlocal response of thermal optical media can be accurately described by a localized well function dependent only on the sample geometry, and not on the nature of the material. This property allows to express the temperature variation, in a stationary limit, as governed by the following 3D heat equation [29, 46, 51,52,53, 59, 79, 80] with constant boundary conditions (at room temperature):

(8) $\left({\mathrm{\partial }}_X^2+{\mathrm{\partial }}_Y^2+{\mathrm{\partial }}_Z^2\right)\mathrm{\Delta }T(\mathbf{R})=-\gamma |\mathbf{E}(\mathbf{R}){|}^2,$(8)

where $\gamma =(L_{loss}{\rho }_0{c}_P{D}_T{)}^{-1}$, $L_{loss}$ is the loss characteristic length, ${\rho }_0$ is the material density, $c_P$ is the specific heat at constant pressure, $D_T$ is the thermal diffusivity and $\mathbf{R}=(X,Y,Z)=({\mathbf{R}}_{\mathrm{\perp }},Z)$. The solution can be written as:

(9) $\mathrm{\Delta }T(\mathbf{R})=\int \int \int \mathrm{d}{Z}^{\prime }\mathrm{d}{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }}G({\mathbf{R}}_{\mathrm{\perp }}-{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }},Z-Z^{\prime })|\mathbf{E}({{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }},Z^{\prime }){|}^2,$(9)

with $G({\mathbf{R}}_{\mathrm{\perp }})$ a Green function that depends only on the sample geometry and the boundary conditions, and expresses the nonlocality of the nonlinear effect. In principle, one can remove the $Z$-dependence by integrating along the longitudinal medium length $Z_0$ [53], but we are interested in the Green function itself, and the longitudinal behavior of $G$ becomes as complicated as $Z_0$ becomes comparable to $L_{loss}$, getting smaller and strongly asymmetric near the boundaries [59]. Physically, the reason why this happens is due to the choice of heat equation to describe the nonlinear radiation–matter interaction: it works only in a neighborhood of the sample midpoint $\hat{Z}=Z_0/2$, not in proximity of the borders. Mathematically, this is deciphered in a longitudinal parabolic approximation with characteristic width $L_{nloc}=\sqrt{\frac{|n_2|}{\gamma {\left|\frac{\mathrm{\partial }n}{\mathrm{\partial }T}\right|}_0}}\propto \sqrt{L_{loss}}$:

(10) $\mathrm{\Delta }T(\mathbf{R})=\left[1-\frac{{(Z-\hat{Z})}^2}{2L_{nloc}}\right]\mathrm{\Delta }{T}_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }}).$(10)

From Eqs. (8,10) we obtain the 2D heat equation:

(11) $\left({\mathrm{\partial }}_X^2+{\mathrm{\partial }}_Y^2\right)\mathrm{\Delta }{T}_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }})-L_{nloc}^{-2}\mathrm{\Delta }{T}_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }})=-\gamma I_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }}),$(11)

with $I_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }})=\frac{1}{Z_0}\int \mathrm{d}Z|\mathbf{E}({\mathbf{R}}_{\mathrm{\perp }},Z){|}^2$. Eq. (9) now reads

(12) $\mathrm{\Delta }{T}_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }})=\int \int \mathrm{d}{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }}{G}_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }}-{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }})I_{\mathrm{\perp }}({{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }}).$(12)

In low absorption regime ($Z_0<<L_{loss}$) $\mathrm{\Delta }T(\mathbf{R})\sim \mathrm{\Delta }{T}_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }})$ and ${\mathrm{\partial }}_{Z}I(\mathbf{R})\sim 0$ (intensity longitudinal changes are negligible as for solitary wave packets); therefore, we attain $n[I](\mathbf{R})=n_0+\mathrm{\Delta }n[I]({\mathbf{R}}_{\mathrm{\perp }})$, with the refractive index nonlocal perturbation:

(13) $\mathrm{\Delta }n[I]({\mathbf{R}}_{\mathrm{\perp }})=n_2\int \int \mathrm{d}{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }}K({\mathbf{R}}_{\mathrm{\perp }}-{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }})I({{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }}),$(13)

and $n_{2}K({\mathbf{R}}_{\mathrm{\perp }})=(\frac{\mathrm{\partial }n}{\mathrm{\partial }T}{)}_0{G}_{\mathrm{\perp }}({\mathbf{R}}_{\mathrm{\perp }})$.

By a comparison between the nonlocality length $L_{nloc}$ and the beam waist $W_0$, we can analyze two different limits: the standard Kerr effect in Eq. (6) when $L_{nloc}\ll W_0$ (local approximation), i.e., $K({\mathbf{R}}_{\mathrm{\perp }}-{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }})\sim \delta ({\mathbf{R}}_{\mathrm{\perp }}-{{\mathbf{R}}_{\mathrm{\perp }}}^{\mathrm{\prime }})$, and the opposite case $L_{nloc}\gg W_0$, that is, the highly nonlocal approximation (HNA), where $K\ast I(\mathbf{R})\sim K({\mathbf{R}}_{\mathrm{\perp }})P(Z)$, with $P(Z)=\int \mathrm{d}{\mathbf{R}}_{\mathrm{\perp }}I(\mathbf{R})$ the power.

For a monocromatic field $\mathbf{E}(\mathbf{R},\tau )={\hat{\mathbf{E}}}_{0}A(\mathbf{R})e^{-i\omega \tau }$ in a third-order thermal medium, in paraxial and slowly varying envelope approximations, introducing the delayed time $\tau =t-\frac{n_0}{c}Z$ and adding a linear loss of characteristic length $L_{loss}$, from Eq. (1) we find that the propagation along $Z$ is ruled by the nonlocal NLSE [29]:

(14) $2ik{\mathrm{\partial }}_{Z}A+\left({\mathrm{\partial }}_X^2+{\mathrm{\partial }}_Y^2\right)A+2k^2\frac{\mathrm{\Delta }n[|A{|}^2]}{n_0}A=-i\frac{k}{L_{loss}}A,$(14)

with $k=\frac{2\pi n_0}{\lambda }=\frac{\omega n_0}{c}$ the wavenumber.

2.1. Spatial dispersive shock waves in nonlocal Kerr nonlinearity

Spatial DSWs are rapidly oscillating waves which regularize an abrupt discontinuity in phase through diffraction, that is, through the formation of intensity undular bores on the beam borders. Scientific community paid close attention to the theoretical description [1,2,22] and experimental demonstration [14, 16, 18–21, 26, 27, 29, 30, 33] of optical DSWs. Here we summarize results on the defocusing DSWs in nonlocal media [29]. In such materials, the IST cannot describe the solutions, and we need other methods.

In next two sections, we detail two different methodologies for DSWs in nonlocal media: the TAQM [12,63–73] and the wave turbulence theory [74,76]. Both theories also prove that DSWs are intrinsically irreversible.

Starting from Eq. (14), through the scaling $x=\frac{X}{W_0}$, $y=\frac{Y}{W_0}$, $z=\frac{Z}{L}$, $\psi (x,y,z)=\frac{A(X,Y,Z)}{\sqrt{I_0}}$, with $I_0$ the intensity peak, $L=\sqrt{L_{nl}{L}_d}$, $L_{nl}=\frac{n_0}{k|n_2|I_0}$ the nonlinear length scale associated to a local Kerr effect, $L_d=k{W}_0^2$ the diffraction length, one obtains the normalized nonlocal NLSE:

(15) $iϵ{\mathrm{\partial }}_z\psi +\frac{{ϵ}^2}{2}\left({\mathrm{\partial }}_x^2+{\mathrm{\partial }}_y^2\right)\psi +\chi \theta \psi =-i\frac{\alpha }{2}ϵ\psi ,$(15)

with $ϵ=\frac{L_{nl}}{L}=\sqrt{\frac{L_{nl}}{L_d}}$ a small quantity in strongly nonlinear (or weakly diffracting) regime, as the one that we are considering, $\chi =\frac{n_2}{|n_2|}$, $\theta =\left|\frac{k\mathrm{\Delta }n{L}_{nl}}{n_0}\right|$, $\alpha =\frac{L}{L_{loss}}$. From Eq. (11)

(16) $-{\sigma }^2\left({\mathrm{\partial }}_x^2+{\mathrm{\partial }}_y^2\right)\theta +\theta =|\psi {|}^2,$(16)

where $\sigma =\frac{L_{nloc}}{W_0}$ is the nonlocality degree, which expresses the nature of the Kerr effect through the limits we have previously discussed: if $\sigma \ll 1$ we are considering the local limit $\theta \sim |\psi {|}^2$, instead, if $\sigma \gg 1$, by HNA $\theta \sim \kappa (x,y)p(z)$, with $p(z)=\int \mathrm{d}x\mathrm{d}y{\left|\psi (x,y,z)\right|}^2=\frac{P(Z)}{W_0^2{I}_0}$.

The fundamental laser mode (Gaussian TEM${ }_{00}$) is described by an axisymmetric Gaussian input ${\psi }_0(r)=exp(-r^2)$, with $r=\sqrt{x^2+y^2}$, which evolves in the WKB approximation [81] as $\psi (r,z)=\sqrt{\rho (r,z)}exp\left[i\frac{\varphi (r,z)}{ϵ}\right]$. For $D=2$ the transverse dimensionality and $u={\mathrm{\partial }}_r\varphi$ the phase chirp, from Eqs. (15, 16) one obtains

(17) $\begin{array}{ccc}{\mathrm{\partial }}_z\rho +\left[\frac{D-1}{r}\rho u+{\mathrm{\partial }}_r(\rho u)\right]& =& -\alpha \rho ,\\ {\mathrm{\partial }}_{z}u+u{\mathrm{\partial }}_{r}u-\chi {\mathrm{\partial }}_r\theta & =& 0,\\ -{\sigma }^2\left[{\mathrm{\partial }}_r^2\theta +\frac{D-1}{r}{\mathrm{\partial }}_r\theta )\right]+\theta & =& \rho .\end{array}$(17)

Figure 1 reports phase chirp and field amplitude for $D=1$, so for ${\mathrm{\partial }}_y\sim 0$ and $r\to x$, in a defocusing medium $(\chi =-1)$ without losses $(\alpha =0)$. The local case $(\sigma =0)$ is illustrated in Figure 1(a,c) and follows from system (17):

(18) ${\mathrm{\partial }}_z\rho +{\mathrm{\partial }}_x(\rho u)=0,{\mathrm{\partial }}_{z}u+u{\mathrm{\partial }}_{x}u=-{\mathrm{\partial }}_x\rho .$(18)

Figure 1. Phase chirp $u(x)$ (a,b), and amplitude $\sqrt{\rho (x,z)}$ (c,d) for transverse dimensionality $D=1$ and different values of $z$, as indicated. (a,c) are obtained by Eqs. (18) with $ϵ={10}^{-3}$. (b,d) are simulations of the result of the system (17) with $D=1,\alpha =0,\chi =-1,{\sigma }^2=5$.

Figure reprinted with permission from [29]. Copyright 2007 by the American Physical Society.

Eqs. (18) are equivalent to Euler and continuity equations, respectively, for a fluid of speed $u$, mass density $\rho$ and pressure proportional to ${\rho }^2$. In the reported dynamics, the diffraction, initially of order ${ϵ}^2$, starts to play a relevant role in proximity of the wave breaking. In fact, it regularizes such a discontinuity by rapid oscillations of wavelength $\sim ϵ$, which appear simultaneously in phase chirp $u$ and intensity $\rho$. For large values of $\sigma$, the normalized refractive index variation, here expressed by $\theta (x)$, is wider than the Gaussian input. As shown in Figure 1(b,d) the shock oscillations are essentially driven by the phase chirp $u$, while the intensity $\rho$ adiabatically follows. Major details are given in [29].

An in-depth description of the difference between DSWs in local and nonlocal Kerr media is also provided by turbulence theory, in particular by the Vlasov formalism, briefly summarized in Section 2.3. The analysis made for random optical waves in Section 2.3 is also relevant to the coherent problem considered here, since the reduced hydrodynamic equations derived from the Vlasov model [39], Eqs. (3,4)] coincide with Eq. (17). Following this approach, DSWs in thermal nonlinearity were interpreted for the first time as ‘annular collapse singularities’ in [39]. By looking at the M-shaped field amplitude in Figure 1d (and, in the following, at the intensity profiles in Figures 4a–6(d,e), Figures 8(a,b) and 16(b,c)), and comparing this to the fast oscillations in Figure 1c [or Figure 2 in [30]], the feature of the collapse singularity in nonlocality appears evident. Indeed, the corresponding hydrodynamic model in the limit of a local nonlinearity [Eq. (18)] recall the shallow water equations, which exhibit a pure shock without collapse.

Figure 2. Pictorial representation of an energy landscape. When the system is in the proximity of a local maximum it obeys the RHO Hamiltonian, in figure ${\hat{H}}_{RO}$. In the proximity of the minimum the system obeys the Hamiltonian of a harmonic oscillator, in figure ${\hat{H}}_{HO}$. The two Hamiltonians are explicitly written in the two the corresponding text boxes, with the related dynamical systems and the discrete eingenvalues. Insets show the transverse profiles of the respective eigenfunctions, bounded on right-hand side for the harmonic oscillator, unbounded on the left-hand side for the RHO.

Reprinted by permission from Macmillan Publishers Ltd. from [38]. Copyright 2015.

2.2. High nonlocality and time asymmetric quantum mechanics

Let us consider the nonlocal NLSE in Eq. (14) with a medium response function $K(X,Y)=exp\left[-(|X|+|Y|)/L_{nloc}\right]/(2L_{nloc}{)}^2$. Being $K$ separable, i.e., $K(X,Y)=\tilde{K}(X)\tilde{K}(Y)$, through the approximation ${\mathrm{\partial }}_Y\sim 0$ (as in the previous section) we can consider only one transverse dimension, since analyzing propagation along $Y$ is no more interesting for our purposes. We rewrite Eq. (14) in terms of $1+1$ dimensionless variables by using the same scaling of Eq. (15) and choosing $I_0$ such that $L_{nl}=L_d$:

(19) $i{\mathrm{\partial }}_z\psi +\frac{1}{2}{\mathrm{\partial }}_x^2\psi -\kappa \ast |\psi {|}^2\psi =-i\frac{\alpha }{2}\psi .$(19)

with $\kappa (x)=W_0\tilde{K}(x{W}_0)=exp\left(-|x|/\sigma \right)/(2\sigma )$.

We take into account a medium where the nonlocality length is much larger than the beam waist. By HNA we have [55,82]

(20) $\kappa \ast |\psi {|}^2\sim \kappa (x)p(z),$(20)

where $\kappa$ is a function no more depending on $|\psi {|}^2$. In a system without loss, that is, $\alpha =0$, the normalized power $p$ is conserved and the NLSE is mapped into a linear Schrödinger equation $i{\mathrm{\partial }}_z\psi =\hat{H}\psi$, with the Hamiltonian $\hat{H}=\frac{1}{2}{\hat{p}}^2+p\kappa (x)$ $\left(\hat{p}=-i{\mathrm{\partial }}_x\right)$. When we express the even function $\kappa$ as its second-order expansion, that is, $\kappa (x)={\kappa }_0^2-\frac{{\kappa }_2^2}{2}{x}^2$, where ${\kappa }_0^2=\frac{1}{2\sigma }$ and ${\kappa }_2^2=\frac{1}{\sqrt{\pi }{\sigma }^2}$, we obtain the reversed harmonic oscillator (RHO) Hamiltonian [68,70,73]:

(21) $\hat{H}=p{\kappa }_0^2+{\hat{H}}_{RHO},{\hat{H}}_{RHO}=\frac{{\hat{p}}^2}{2}-\frac{{\gamma }^2{\hat{x}}^2}{2},{\gamma }^2=p{\kappa }_2^2.$(21)

If $\psi =exp\left(-i{\kappa }_0^{2}pz\right)\varphi$, then $i{\mathrm{\partial }}_z\varphi ={\hat{H}}_{RHO}\varphi$.

Figure 2 sketches the relation between the harmonic and the reversed oscillators. For a harmonic oscillator (HO), the spectrum is discrete and the corresponding eingenstates form an orthonormal basis (both an orthogonality and a completeness relations hold):

(22) $\begin{array}{c}{\hat{H}}_{HO}=\frac{{\hat{p}}^2}{2}+\frac{{\omega }^2}{2}{\hat{x}}^2\\ {\hat{H}}_{HO}\psi (x)=E\psi (x),E_n=\omega \left(n+\frac{1}{2}\right),\\ {\psi }_n(x)=\sqrt[4]{\frac{\omega }{\pi }}\frac{1}{\sqrt{2^{n}n!}}{H}_n\left(\sqrt{\omega }x\right),\end{array}$(22)

with $H_n(x)=(-1{)}^n{x}^2\frac{d^n}{d{x}^n}{e}^{-x^2}$ the Hermite polynomials. On the other hand, RHO has complete continuous spectrum, but one derives a generalized discrete spectrum from HO spectrum by a complex analytic prolongation in the rigged Hilbert space [67,73,83] through the transformation $\omega \to i\gamma$, $\hat{x}\to e^{-i\frac{\pi }{4}}\hat{x}$, $\hat{p}\to e^{i\frac{\pi }{4}}\hat{p}$ [70,73]. The new stationary Schrödinger equation is ${\hat{H}}_{RHO}{\mathfrak{f}}^{\pm }(x)=i\frac{\mathrm{\Gamma }}{2}{\mathfrak{f}}^{\pm }(x)$, solved by the spectrum $\frac{{\mathrm{\Gamma }}_n}{2}=\gamma \left(n+\frac{1}{2}\right)$ and the non-normalizable eigenfunctions

(23) ${\mathfrak{f}}_n^{\pm }(x)=\frac{\sqrt[4]{\pm i\gamma }}{\sqrt{2^{n}n!\sqrt{\pi }}}{H}_n(\sqrt{\pm i\gamma }x)exp(\mp i\frac{\gamma }{2}{x}^2),$(23)

namely, the RHO Gamow vectors (GVs) [84,85].

We can express every wavefunction as a truncated superposition of GVs added to a background function, which dispersively oscillates at infinite as a polynomials [73]:

(24) $\varphi (x)={\varphi }_N^G(x)+{\varphi }_N^{BG}(x)$(24)

with

(25) ${\varphi }_N^G(x)=\sum _{n=0}^N{\mathfrak{f}}_n^{-}(x)⟨{\mathfrak{f}}_n^{+}|\varphi (x,0)⟩.$(25)

Figure 3 shows the GV square norms (Figure 3a) and phase chirps (Figure 3b). The evolution of the normalized field $\psi$ presents a Gamow part resulting as a superposition of exponential decays with quantized decay rates [73]:

(26) ${\psi }_n^G(x,z)=\sum _{n=0}^N⟨{\mathfrak{f}}_n^{+}|\psi (x,0)⟩{\mathfrak{f}}_n^{-}(x)e^{-i{\kappa }_0^{2}pz}{e}^{-\frac{{\mathrm{\Gamma }}_n}{2}z}.$(26)

Figure 3. (a) $||{\mathfrak{f}}_n^{-}(x)|{|}^2$ in Eq. (23) for increasing even order $n$; (b) corresponding phase chirps ${\mathrm{\partial }}_x\mathrm{A}\mathrm{r}\mathrm{g}\left[{\mathfrak{f}}_n^{-}(x)\right]$; (c) weights $p_n(0)$ [Eq. (27)] of the GV expansion of a Gaussian wave packet.

Figure reprinted with permission from [37]. Copyright 2015 by the American Physical Society.

Eq. (26) proves an intrinsical irreversibility of DSWs, where a backward propagation beyond the shock point is no physically possible because of the exponentially decaying evolution. This explains why the quantum representation of wave propagation theory in a rigged Hilbert space is called TAQM (here time is replaced by $z$).

In the probabilistic interpretation of TAQM [37], the projection of Eq. (26) over $\sqrt{{\mathrm{\Gamma }}_n}{\mathfrak{f}}_n^{+}$ gives the probability $p_n(z)$ of finding the system in a decaying GV:

(27) $p_n(z)={\mathrm{\Gamma }}_n\left|⟨{\mathfrak{f}}_n^{+}\right|\psi (x,0)⟩{|}^2{e}^{-{\mathrm{\Gamma }}_{n}z},$(27)

which gives the $z-$dependent weight of the $n$-order GV. Initial weights $p_n(0)$ are reported in Figure 3c as functions of ${\gamma }^2$. Since a Gaussian beam $\psi (x,0)=\phi (x)=exp(-x^2/2)/\sqrt[4]{\pi }$ is an even input, all the odd terms in Eq. (26) vanish due to the $x-$parity. Figure 4a shows the numerical solution of Eq. (19). Yellow lines give the transverse intensity profile. We see that these are modeled by a superposition of exponential decays, where the plateau is given by the groundstate GV, and the peaks are given by higher-order GVs. Simulations of weights $p_n(z)$ are in Figure 4(b,c). While dotted profiles are numerical results from Eq. (27), continuous lines result from the general projection definition $p_n(z)={\mathrm{\Gamma }}_n\left|⟨{\mathfrak{f}}_n^{+}\right|\psi (x,z)⟩{|}^2$, with $\psi (x,z)$ numerical solution of Eq. (19).

Figure 4. (a) Numerical solution of Eq. (19) with $p={10}^4$ and ${\sigma }^2=10$; (b) projection on GVs for increasing order $n$ for $\alpha =0.3$ and $\gamma =8$; continuous lines are from Eq. (19), dots are from Eq. (27); (c) as in panel (b) for $\gamma =24$.

Figure reprinted with permission from [37]. Copyright 2015 by the American Physical Society.

2.3. Random optical waves and turbulence theory

We have summarized how TAQM describes the intrinsic irreversibility of DSWs generated by a laser input. We highlight that nonlinear optics offers other examples of intrinsical irreversibility, as photon fluids, i.e., turbulent flows of a conservative system of random optical waves, in which the statistical interpretation has produced significant results [76].

Statistical nonlinear optics is related to wave turbulence theory, whereby the kinetic wave description provides a thermodynamic treatment of turbulence [33,39,41,74–76,86–89]. If one considers the nonlinear propagation of incoherent optical waves characterized by fluctuations that are statistically inhomogeneous in space, DSWs arise only in the strong turbulence (strongly nonlinear) regime [39,41]. The wave breaking emerges from a turbulent field, whose local-averaged spectrum follows a specific Vlasov-like kinetic equation, which can be a traditional Vlasov equation (local nonlinearity), a short-range Vlasov equation (SRVE) (quasi-local nonlinearity) or a long-range Vlasov equation (LRVE) (highly nonlocal nonlinearity), derived by the zero-loss NLSE through a multiscale expansion [76]. These Vlasov equations, being reversible kinetic equations, do not exhibit a H-theorem of entropy growth, and so do not describe the process of irreversible thermalization to equilibrium. They can be interpreted as a mean-field theory that does not explain, by itself, irreversible processes. For the sake of completeness, we stress that the thermalization process arises at the next order of a weakly nonlinear expansion procedure. This analysis leads to a different kind of time asymmetric behavior with respect to DSWs reported above, because it is not due to strong nonlinearity. In this framework, by going to next order, the theory reveals that the Vlasov equation is corrected by a collision term involving the nonlocal interaction. This collision term has a form analogous to the wave turbulence kinetic equation [86] and is written in explicit form in Methods of [39]. The mathematical statement of such irreversibility relies on the H-theorem of entropy growth. Wave thermalization can be characterized by a self-organization process, that is, the system spontaneously generates a large-scale coherent structure. A remarkable example of this self-organization process is the wave condensation, whose thermodynamic equilibrium properties are similar to those of quantum Bose–Einstein condensates. We stress that nonlocality significantly decelerates the rate of the thermalization process, in complete analogy with gravitational systems, where the dynamics of stars are described by a collision-less (Vlasov-like) equation toward quasi-stationary nonequilibrium states (e.g., galaxies), which are of fundamental different nature than usual thermodynamic equilibrium states [90]. A detailed treatise of these phenomena is reported in [39] and in the references therein.

In what follows, we analyze long-range interactions in strongly nonlinear wave systems operating far from thermodynamic equilibrium. Starting from the NLSE [Eq. (14)], with nonlinearity expressed by Eq. (13) and no loss, we obtain Eq. (15), here explicitly written as:

(28) $i{\mathrm{\partial }}_z\psi +\frac{ϵ}{2}\left({\mathrm{\partial }}_x^2+{\mathrm{\partial }}_y^2\right)\psi +\chi \upsilon \psi U\ast |\psi {|}^2=0,$(28)

with $\upsilon =\frac{|n_2|I_0}{n_0}k\sqrt{L_{nl}{L}_d}$, $U(\mathbf{r})=W_0^{2}K(W_0\mathbf{r})$, $\mathbf{r}=(x,y)$. From Eq. (28), we want to attain a kinetic equation, that is, an equation describing the evolution of the spectrum during the related field propagation in the nonlinear medium. The structure of a kinetic equation depends on the statistics of the random wave. For this reason, we consider the field autocorrelation function:

(29) $\begin{array}{ccc}B(\mathbf{r},\xi ,z)=⟨\psi (\mathbf{r}+\frac{1}{2}\xi ,z){\psi }^{\ast }(\mathbf{r}-\frac{1}{2}\xi ,z)⟩,& \mathbf{r}=({\mathbf{r}}_{\mathbf{1}}+{\mathbf{r}}_{\mathbf{2}})/2,& \xi ={\mathbf{r}}_{\mathbf{1}}-{\mathbf{r}}_{\mathbf{2}}.\end{array}$(29)

The statistics is said to be homogeneous if $B$ depends only on the distance $|{\mathbf{r}}_{\mathbf{1}}-{\mathbf{r}}_{\mathbf{2}}|$. Following Eq. (28):

(30) $i{\mathrm{\partial }}_{z}B(\mathbf{r},\xi ,z)+ϵ{\mathrm{\nabla }}_{\mathbf{r}}\cdot {\mathrm{\nabla }}_{\xi }B(\mathbf{r},\xi ,z)+\chi \upsilon \left[P(\mathbf{r},\xi ,z)-Q(\mathbf{r},\xi ,z)\right]=0,$(30)

with

(31) $\begin{array}{rl}& P(\mathbf{r},\xi ,z)=B(\mathbf{r},\xi ,z)\int \mathrm{d}{\mathbf{r}}^{\prime }U({\mathbf{r}}^{\prime })\left[N(\mathbf{r}-{\mathbf{r}}^{\prime }+\frac{1}{2}\xi ,z)-N(\mathbf{r}-{\mathbf{r}}^{\prime }-\frac{1}{2}\xi ,z)\right],\\ & Q(\mathbf{r},\xi ,z)=\int \mathrm{d}{\mathbf{r}}^{\prime }U({\mathbf{r}}^{\prime })\left[B(\mathbf{r}-\frac{1}{2}{\mathbf{r}}^{\prime }+\frac{1}{2}\xi ,{\mathbf{r}}^{\prime },z)B(\mathbf{r}-\frac{1}{2}{\mathbf{r}}^{\prime },\xi -{\mathbf{r}}^{\prime },z)\right.\\ & \left.+-B(\mathbf{r}-\frac{1}{2}{\mathbf{r}}^{\prime },\xi +{\mathbf{r}}^{\prime },z)B(\mathbf{r}-\frac{1}{2}{\mathbf{r}}^{\prime }-\frac{1}{2}\xi ,-{\mathbf{r}}^{\prime },z)\right],\\ & N(\mathbf{r},z)=B(\mathbf{r},\xi =\mathbf{0},z)=⟨|\psi {|}^2⟩(\mathbf{r},z).\end{array}$(31)

In the last equation, $N$ is the averaged field power, also depending on $\mathbf{r}$ because of the inhomogeneity of the statistics. By defining the length scale of random fluctuations as the coherence length ${\lambda }_c$ and $W_0$ the incoherent beam waist, we can assume that the statistics is quasi-homogeneous if ${ϵ}_c={\lambda }_c/W_0<<1$. For a Gaussian response function $U(\mathbf{r})=\frac{1}{2\pi {\sigma }^2}exp\left(-\frac{|\mathbf{r}{|}^2}{2{\sigma }^2}\right)$, we get the SRVE if $\sigma <<1$ and the LRVE if $\sigma >>1$ through two different multiscale expansions with respect to ${ϵ}_c$:

(32) ${\mathrm{\partial }}_z{n}_{\mathbf{k}}(\mathbf{r},z)+{\mathrm{\nabla }}_{\mathbf{k}}{\tilde{\omega }}_{\mathbf{k}}(\mathbf{r},z)\cdot {\mathrm{\nabla }}_{\mathbf{r}}{n}_{\mathbf{k}}(\mathbf{r},z)-{\mathrm{\nabla }}_{\mathbf{r}}{\tilde{\omega }}_{\mathbf{k}}(\mathbf{r},z)\cdot {\mathrm{\nabla }}_{\mathbf{k}}{n}_{\mathbf{k}}(\mathbf{r},z)=0,$(32)

for the local spectrum $n_{\mathbf{k}}(\mathbf{r},z)=\int \mathrm{d}\xi B(\mathbf{r},\xi ,z)exp(-i\mathbf{k}\cdot \xi )$, with averaged power $N(\mathbf{r},z)=\frac{1}{{(2\pi )}^2}\int \mathrm{d}\mathbf{k}{n}_{\mathbf{k}}(\mathbf{r},z)$. The nonlocal features of Eq. (32) are traced by the generalized dispersion relation. Once defined the Fourier transform of the response function $\tilde{U}(\mathbf{k})=\int \mathrm{d}\mathbf{r}U(\mathbf{r},z)exp(-i\mathbf{k}\cdot \mathbf{r})$, for the LRTE:

(33) $\begin{array}{rl}& {\tilde{\omega }}_{\mathbf{k}}(\mathbf{r},z)=\omega (\mathbf{k})+V(\mathbf{r},z),\\ & \omega (\mathbf{k})=\frac{ϵ}{2}|\mathbf{k}{|}^2,\\ & V(\mathbf{r},z)=-\chi \upsilon \int \mathrm{d}{\mathbf{r}}^{\prime }U(\mathbf{r}-{\mathbf{r}}^{\prime })N({\mathbf{r}}^{\prime },z).\end{array}$(33)

One can also compute the momentum both from the NLSE and the LRVE:

(34) $\begin{array}{ccc}{\mathbf{p}}_{NLSE}(\mathbf{r},z)& =& {(2\pi )}^2\mathrm{\Im }\left({\psi }^{\ast }\mathrm{\nabla }\psi \right)/N(\mathbf{r},z),\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ {\mathbf{p}}_{LRVE}(\mathbf{r},z)& =& \int \mathrm{d}\mathbf{k}{n}_{\mathbf{k}}(\mathbf{r},z)\mathbf{k}/N(\mathbf{r},z).\end{array}$(34)

Major details are given in [76].

By increasing the range of the nonlocality, we pass from a stage where the field evolution is ruled by stochastic generation of small-scale DSW structures, naturally denoted as dispersive shocklets, to the emergence of an unexpected global collective behavior. The latter phenomenon is characterized by a strong non-homogeneous redistribution of the spatial fluctuations, whose description is provided by the NLSE and the LRVE, and unveils the formation of a giant shock singularity. From theoretical analysis [41], it turns out that – in the short-range regime – wave breaking occurs at random positions in the turbulent field, predominantly around high-amplitude fluctuations, leading to a gas of coherent dispersive shocklets in the midst of turbulent fluctuations. On the other hand, for highly nonlocal Kerr samples, the regularization of the global incoherent shock does not require the formation of a regular oscillating DSW structure because the self-organization of the turbulent waves ensemble makes the field oscillate as a whole. The momentum of the speckled beam [Eq. (33)] is radially outgoing and exhibits a shock-like singularity, while the envelope of the intensity of the beam experiences a ring-shaped collapse-like behavior. The fluctuations of the incoherent wave then result to be pushed towards the annular shock front, which leaves behind itself an internal region of the beam with a high degree of coherence. In other terms, the dynamics is featured by a dramatic degradation of the coherence properties on the annular boundary of the beam, while its internal region exhibits a significant coherence enhancement. Experimental observations of incoherent DSWs are reported in Section 3.4, together with related numerical simulations.

This alternation between coherence degradation on the boundaries and relative enhancement on the internal region also emerges in the spectrograms achievable by LRVE simulations. It turns out that, as the pump power increases, the spectrogram is affected by a Z-shaped distortion: the coherence length ${\lambda }_c\sim \mathrm{\Delta }{k}^{-1}$ decreases at shock front ($\mathrm{\Delta }{k}_{shock}$ determines the boundaries of the admitted $k$ values), and it increases in the internal region of the beam. This dynamics is conservative but irreversible, because the non-equilibrium thermodynamic condition is stabilized by an irreversible self-organization of the random waves. Indeed, in [77] authors showed that a spontaneous long-range phase coherence emerges among incoherent waves in nonlocal nonlinearity, when the speckles self-organize into giant collective waves. Their theory reveals that this phenomenon constitutes a generic property of a conservative (Hamiltonian) system of highly nonlocal random waves that evolve in a strongly nonlinear regime. Moreover, the field exhibits intensity fluctuations whose coherence length increases dramatically during the evolution of the system.

3. Experimental observations in thermal media

Shock waves described by Eq. (14) have been originally shown in an experiment from [29]. The sample is a cell of length 1 mm filled with an aqueous solution of rhodamine B (RhB), with a concentration of 0.6mM. Measurements of the shock profiles are in Figure 5. A Gaussian CW laser beam of intensity waist $W_0=20\mu$m, at wavelength $\lambda =532$nm, propagates in a material with linear refractive index $n_0=1.3$, defocusing Kerr coefficient $n_2=-7\times {10}^{-7}$cm${ }^2$W${ }^{-1}$, loss length $L_{loss}^{-1}=62$cm${ }^{-1}$. For water $D_T=1.5\times {10}^{-7}$m${ }^2$s${ }^{-1}$, ${\rho }_0={10}^3$kg m${ }^{-3}$, $c_p=4\times {10}^3$J kg${ }^{-1}$K${ }^{-1}$, ${\left|\frac{\mathrm{\partial }n}{\mathrm{\partial }T}\right|}_0={10}^{-4}$K${ }^{-1}$. The degree of nonlocality is estimated as $\sigma =0.3$.

Figure 5. Experimental transverse intensity profiles of an initial Gaussian beam propagating in a thermal medium. Measurements are performed for varying input power $P=\pi W_0^2{I}_0$. Insets show the 2D output patterns.

Figure reprinted with permission from [29]. Copyright 2007 by the American Physical Society.

The beam exhibits, beyond the shock point, the formation of undular bores moving outward with increasing power. The next subsections report different experiments exhibiting DSWs in nonlocal samples.

3.1. Rhodamine and time asymmetric quantum mechanics interpretation

In this section we report two experiments in order to validate the presence of GVs in DSWs: a 2D propagation pattern to observe GVs decay rates ${\mathrm{\Gamma }}_n$ [38], and a 1D experiment to show that GVs describe also the M-shaped profile in the far-field of a DSW in HNA [40], in [39] identified for the first time as collapse singularity. These are validations of TAQM in describing DSW propagation.

The experimental setup is illustrated in Figure 6a. Samples are prepared by dispersing 0.1mM of RhB in water. The solution is placed in a cuvette 1 mm thick in the propagation direction. The measured defocusing Kerr coefficient is $|n_2|=2\times {10}^{-12}$m${ }^2$W${ }^{-1}$ and the absorption length is $L_{loss}\simeq 1.6$mm at the laser wavelength $532$nm [32]. The CW laser beam is focused through a lens into a sample. Light is collected by a spherical lens and a Charged Coupled Device (CCD) camera. A microscope is placed above the sample in order to capture top-view images of the laser beam along the propagation direction $Z$. The difference between the two experimental apparatus is the choice of the first lens (L1). In the 2D experiment [38], L1 is spherical with focal length 10 cm, and a focus spot size of $10\mu$m. The setup was placed having the beam propagating vertically through the sample, reducing thermal convection in the water. In the 1D experiment [40], authors used a cylindrical lens as L1, with focal length $f=20$cm in order to mimic a nearly one-dimensional propagation. Being $Z$ the propagation direction, the lens focuses the beam in the $X$ direction. The input spot dimension is $1.0$mm in the $Y$ direction and $35\mu$m in the $X$ direction. These geometrical features make the one-dimensional approximation valid and allow us to compare experimental results with the theoretical one-dimensional model. The diffraction length in the $X$ direction is $L_d=3.0$mm. This time, the setup was placed horizontally.

Figure 6. (a) Experimental setup. Authors of [38,40] collected the transmitted and fluorescence images of the laser beam propagating in RhB samples. Two types of launching lenses L1 were used: a cylindrical and a spherical, for the 1D and 2D experiments, respectively. The top fluorescence image of the propagating beam was collected by a microscope placed above the RhB samples. The second lens is spherical and was used to collect the transverse output profile. (b, c) Top-view intensity distribution as obtained from 2D experiment (b) and numerical simulations (c). Respectively experimental (d) and numerical (e) sections of the images (b) and (c) taken at $z=0.2$ (red), 0.6 (green) and 0.9 mm (blue).

Reprinted by permission from Macmillan Publishers Ltd. from [38]. Copyright 2015.

Figure 6(b,c) reports the observed laser beam propagation top-view, detected by a microscope through RhB fluorescence, and the numerical calculation from the NLSE, respectively. The beam displays the characteristic strongly defocusing and the M-shaped behavior, also evident in the transverse sections of the intensity in Figure 6(d,e). These are signatures of DSWs in nonlinear media at high power.

Decay rates in Figure 7 are detected by slicing the intensity profile $I(X,Z)$ at $X\simeq 0.1$mm (yellow line in Figure 6b) and fitting the intensity versus $Z$ with two exponential functions. Different power levels exhibit very different dynamics. The presence of double-exponential decays, that is, the superposition of the first two GVs, is more evident at high power. It was observed and calculated that double-exponential decay dynamics obeys the quantized spectrum scaling ${\mathrm{\Gamma }}_2/{\mathrm{\Gamma }}_0=5$ at all investigated power levels, as shown in Figure 7d. This demonstrates that authors of [38] excited the fundamental state ${\mathfrak{f}}_0^{-}$ and the first excited state ${\mathfrak{f}}_2^{-}$. Odd states are not excited, as expected from Gaussian TEM${ }_{00}$ $x-$parity. Each of the two rates has a square root dependence on $P$, signature of the underlying nonlinearity. This power dependence distinguishes RHO dynamics from linear loss, due to absorption and scattering.

Figure 7. (a) Observed intensity decay at different laser powers, obtained by slicing along $X\simeq 0.1$mm the top-view intensity distribution the propagation direction (see the yellow line in Figure 6b). (b) Numerically calculated decays in the conditions of panel A. (c) Peak region of the experimental curve at $P=450$mW. The superposition of the first two exponential decays unveils the presence of two GVs, the fundamental state, $n=0$ (slowly decaying) and the first excited state, $n=2$ (fastly decaying). (d) Decay rates vs $P$ for the fundamental state, ${\mathrm{\Gamma }}_0$ (filled circles) and the excited state, ${\mathrm{\Gamma }}_2$, (triangles).

Reprinted by permission from Macmillan Publishers Ltd. from [38]. Copyright 2015.

The RHO eigenstates are quasi-eigenstates of the Fourier transform operator, which in optics represents the far-field. Let us consider the RHO Hamiltonian in the momentum basis ($\hat{p}\to \mathfrak{p}$ and $\hat{x}\to i{\mathrm{\partial }}_{\mathfrak{p}}$)

(35) ${\hat{H}}_{RHO}(\mathfrak{p},i{\mathrm{\partial }}_{\mathfrak{p}})=\frac{{\mathfrak{p}}^2}{2}+\frac{1}{2}{\gamma }^2{\mathrm{\partial }}_{\mathfrak{p}}^2=-{\hat{H}}_{RHO}(-i{\mathrm{\partial }}_x,x).$(35)

Pure GVs are unfeasible to describe a physical experiment, because one cannot neglect that GVs have an infinite support, i.e., the $x$-region where the eigenfunction is not null, is not finite. In order to account for the spatial confinement of the experiment, authors of [40] introduced the windowed GVs:

(36) ${\varphi }_G^W(x)=\sum _{n=0}^N\sqrt{{\mathrm{\Gamma }}_n}{\mathfrak{f}}_n^{-}⟨{\mathfrak{f}}_n^{+}|\psi (x,0)⟩\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}{}_W(x),$(36)

where

$\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}{}_W(x)=\left\{\begin{array}{cc}0& \mathrm{f}\mathrm{o}\mathrm{r}|x|\ge W\\ 1& \mathrm{f}\mathrm{o}\mathrm{r}|x|<W\end{array}\right.,$

with $W$ is the finite size of the physical system. During the evolution, the Gamow ground state has the lowest decay rate, i.e., $\gamma /2$. This allows to consider, in the long-term evolution, only the fundamental GV, and so only the Fourier transform $\mathcal{F}$ of the fundamental state of Eq. (36):

(37) $\tilde{\psi }(k_x)=F\left[{\mathfrak{f}}_0^{-,W}(x)\right]=\left(\frac{1}{4}+\frac{i}{4}\right)e^{-\frac{ik{x}^2}{2\gamma }}\frac{{(-i\gamma \pi )}^{1/4}}{W}\times \left\{-\text{Erf}\left[\frac{\left(\frac{1}{2}-\frac{i}{2}\right)(k_x-W\gamma )}{\sqrt{\gamma }}\right]+\text{Erf}\left[\frac{(\frac{1}{2}-\frac{i}{2})(k_x+W\gamma )}{\sqrt{\gamma }}\right]\right\}.$(37)

Eq. (37) provides an analytical expression of the far-field, which is compared below with the experiments. Indeed, Eq. (37) allows to predict in closed form the typical M-shaped shock profile: it describes the internal undular bores and the correct scaling of the undulation period with respect to the power, i.e., the period $T$ is predicted to scale with the square root of $\gamma$, and hence with the fourth square root of the beam input power.

Figure 8 reports experimental results in RhB, through the previously described setup, and the comparison with the numerical results. Images of the beam in the far-field (corresponding to the square modulus of the spatial intensity Fourier transform) for different input powers were collected and shown in Figure 8(a,b). For low power (not reported) the elliptical beam profile remains Gaussian along propagation. A different phenomenon occurs while increasing the power: the beam transverse section along $X$ broadens and develops intensity peaks on its lateral edges. Essentially, it becomes M-shaped. These results are in remarkable agreement with Eq. (37), as shown in Figure 8(c,d).

Figure 8. (a, b) CCD image of the light beam at laser powers $P=2$W and $4$W, respectively; the bottom panels show the normalized intensity profile at the maximum waist along $Y=0$. (c) Analytical solution obtained by Eq. (37) changing Gaussianly the power $P$ in the $y$ direction. (d) As in (c) but for higher powers; the bottom panels show the slice of panel (c) and (d) at $y=0$, i.e., Eq. (37) square modulus for $W=1.5$ and $\gamma \simeq 12$ and $\gamma \simeq 40$, respectively. (e) Log-scale normalized intensity as a function of power, as obtained by slicing along $Y$ a region in panel (b). The slopes of the straight lines give the GV decay rates (${\gamma }_1=-8\pm 0.4$ and ${\gamma }_2=-1.6\pm 0.1$). Their quantized ratio is $5.0\pm 0.4$ as expected from theory [37]. (f) Intensity oscillations for different power values. (g) Measured oscillations period $T$ as a function of power; continuous line is the fit function $T\propto \sqrt[4]{P}$, as expected by the theory; the inset shows the same curve of (g) with $P^{1/4}$ as abscissa axis.

Reprinted with permission from [40]. Copyright 2016 Optical Society of America.

Different positions in the $Y$ direction correspond to different power levels. Any power level furnishes a different value of $\gamma$, being $\gamma =\sqrt{\frac{p}{\sqrt{\pi }{\sigma }^2}}$. The Gaussian beam profile in the $Y$ direction, that is, $p\propto exp(-y^2)$, provides the link between $Y$, $P$ and $\gamma$. This implies that, observing a CCD image, intensity profiles at different $Y$ correspond to different powers. Therefore, the expected exponential trend with respect to the power can be extracted from a single picture by looking at different $Y$ positions. Figure 8e exhibits a fitting with two exponential decays in an intensity profile versus power. The extracted ratio of the related two decay rates is 5 and hence in agreement with the expected quantized theoretical values described in Section 2.2.

Undular bores of DSWs were analyzed and exhibited in Figure 8f, while the field intensity undulation period $T$ versus $P$ is shown in Figure 8g. In order to demonstrate univocally that $T\propto \sqrt[4]{P}$, inset in Figure 8g reports the period $T$ as function of $\sqrt[4]{P}$. The resulting linear behavior confirms the theoretical results.

3.2. Nonlinearity and disorder in thermal media

Thermal media have been investigated also in their interplay with disorder. Theoretical studies demonstrated that, even if solitons are stable under a certain amount of randomness, the latter competes with nonlinearity, while nonlocality filters disorder-induced scattering effects and soliton random walk can be efficiently suppressed in highly nonlocal media [54,56,91]. DSWs are nonlinear coherent oscillations, and the phenomenon of light scattering affects their formation in significant way [32].

In this section, we report experiments in two different optical systems that combine third-order nonlinearity (high-power laser beams) with nonlocality (thermal material response) and disorder (scattering particles). The first thermal medium is a dispersion of silica spheres of $1\mu$m diameter in 0.1mM aqueous solution of RhB. The second one is a $1$mm$\times 1$mm$\times 8.5$mm parallelepiped of silica aerogel. Despite observations of DSWs in disordered thermal media, a theoretical model that comprehends both nonlinearity, nonlocality, and disorder has been developed only for solitons [54]. The existing theoretical model for DSWs is summarized below and neglects the nonlocality contribution. It approximates thermal nonlinearity to a local Kerr effect, and adds a random potential [32].

We start from Eq. (14) with $\mathrm{\Delta }n[|A{|}^2]=n_2|A{|}^2+\mathrm{\Delta }{n}_R(X,Y,Z)$ and $L_{loss}\sim \mathrm{\infty }$ (no loss). Through the same scaling of Section 2.1 and approximation to cylindrical symmetry ${\mathrm{\partial }}_y\sim 0$, we obtain

(38) $iϵ{\mathrm{\partial }}_z\psi +\frac{{ϵ}^2}{2}{\mathrm{\partial }}_x^2\psi -|\psi {|}^2\psi +U_R\psi =0,$(38)

with $U_R(x,y,z)=\frac{\mathrm{\Delta }{n}_R(X,Y,Z)}{n_2{I}_0}$ taken as a random dielectric noise mainly acting on the phase [32]. In the hydrodynamic limit $ϵ\sim 0$, the phase chirp behaves like a moving unitary mass particle [32]:

(39) $\frac{\mathrm{d}{}^{2}x}{\mathrm{d}{z}^2}=-\frac{\mathrm{d}U}{\mathrm{d}x}+{\eta }_R,$(39)

with $U=exp\left(-x^2/2\right)$ the deterministic potential for a Gaussian TEM${ }_{00}$ given by the nonlinearity, and ${\eta }_R=-\frac{\mathrm{d}{U}_R}{\mathrm{d}x}$ a Langevin force with Gaussian distribution, such that $⟨{\eta }_R(z){\eta }_R(z^{\prime })⟩={\eta }^2\delta (z-z^{\prime })$ and $\eta =\sqrt{⟨{\left(\frac{\mathrm{d}{U}_R}{\mathrm{d}x}\right)}^2⟩}\simeq \sqrt{⟨{\left(\mathrm{\Delta }{n}_R\right)}^2⟩}(|n_2|I_0{)}^{-1}$ the disorder strength. Brackets $⟨,⟩$ denote the statistical average, and the dependence of ${\eta }_R$ on $x,y$ is neglected for stochastic independence and cylindrical symmetry, respectively, thus ${\eta }_R\simeq {\eta }_R(z)$.

Figure 9 shows trajectories $x(z)$ (Figure 9(a,b)) and phase space $(x,v)$ (Figure 9(c,d)), where $v=\frac{\mathrm{d}x}{\mathrm{d}z}$, respectively, without ($\eta =0$) and with ($\eta =0.1$) disorder, the latter obtained by a stochastic Runge-Kutta algorithm [92,93]. In absence of disorder (Figure 9(a,c)) the shock is signaled by the intersection of multiple trajectories $x(z)$ and, in the phase space, this corresponds to the induced wave-breaking phenomenon, that is, the folding of the velocity profile into a multivalued function for increasing $z$. In the presence of disorder, Figure 9(b,d), the particle-like dynamics tends to diffuse, as is evident from the related trajectories and phase space. Correspondingly, the propagation distance before the intersections is greater for the disordered case and the shock is delayed in the $z$ direction.

Figure 9. (a, b) Trajectories $x(z)$ and (b, d) phase space $(x,v)$, respectively, with disorder strength $\eta =0$ and $\eta =0.1$. $z$ varies from $z=0$ to $z=3$.

Figure reprinted with permission from [32]. Copyright 2012 by the American Physical Society.

We report here the experiments in RhB with silica sphere dispersions [32]. In order to vary the degree of disorder, several silica concentrations were prepared, ranging from $0.005$w$/$w to $0.03$w$/$w, in units of weight of silica particles over suspension weight. The experimental setup is similar to that illustrated in Figure 6a. The first lens focuses the beam on the input facet of the sample, reaching a beam waist $W_0\simeq 10\mu$m. The aqueous solutions are put in 1mm$\times 1$cm$\times 3$cm glass cells with propagation along the 1 mm vertical direction (parallel to gravity) to moderate the effect of heat convection. All measurements are performed after the temperature gradient has reached the stationary state and the particle suspensions are completely homogeneous. In [32], the main loss mechanisms are absorption and scattering. The measured loss length (absorption plus scattering) varies in a range from $1.2$mm to $1.6$mm (highest value is for pure dye solution). These values are obtained by fitting with exponential decay the beam intensity vs propagation distance $Z$. The fact that the loss length is always greater than the position of the shock point [32] allowed authors to neglect losses at a first approximation in their theory. In addition, they found that the scattering mean free path is of the order of millimeters for all the considered samples. In Figure 10 images of the transmitted beam on the transverse plane for different input laser powers $P$ and various concentrations $c$ are shown. The number and the visibility of the DSW oscillations increase with $P$ and decrease with $c$, evidence of DSWs enhancement by nonlinearity and inhibition by disorder.

Figure 10. Transverse intensity patterns for different input power $P$ and silica spheres concentration $c$: (a) $P=5mW$, $c=0$w$/$w, (b)$P=400mW$, $c=0$w$/$w, (c) $P=5mW$, $c=0.017$w$/$w, (d) $P=400mW$, $c=0.017$w$/$w, (e) $P=5mW$, $c=0.030$w$/$w, (f) $P=400mW$, $c=0.030$w$/$w. White 1D curves show the measured section of the intensity profiles vs $X$.

Figure reprinted with permission from [32]. Copyright 2012 by the American Physical Society.

Experimental observations have been also performed in silica aerogel [35]. The silica aerogel samples are prepared following a base-catalyzed sol-gel procedure [94], and in-depth details are given in [35]. It turns out that the sample used in the experiment has mass density $\rho =0.215$g$/$cm${ }^3$ and refractive index $n_0=1.074$. Experimental setup is very similar to the previous ones (Figure 1a), except for the sample. In [35], authors vary the input beam waist $W_0$, the input laser power $P_{in}$, and record the transmitted intensity distribution $I(X,Y,Z=8.5\mathrm{m}\mathrm{m})$ by the CCD camera. Observations are shown in Figure 11. Images in the second and third rows of Figure 11 correspond to the same experimental conditions in term of incident laser power and beam size, but the incident laser beam impinges on different points. In correspondence of regions of the silica aerogel sample displaying low enough disorder (second row), a transition from scattering dominated regimes to nonlinear regimes is present: at moderate powers DSWs are not observed because of scattering losses, at high powers DSWs can be generated.

Figure 11. Far-field intensity profiles at the output of the silica aerogel for $P_{in}$ ranging from $1$mW to $1$W, and input beam waist $w_0$ ranging from $43\mu$m to 1.4 mm. Images in the second and third rows correspond to the same incident laser power and beam size, but different positions of the incident laser beam.

Reprinted with permission from [35]. Copyright 2014 Optical Society of America.

3.3. Dispersive shock waves in biological suspensions and chemical compounds

The study of optical effects in light propagation through chemical and biological solutions is a field of growing interest [36,43,58,95–98], both from a linear and a nonlinear perspective. However, although observations of nonlinear optical phenomena in chemical and soft-matter systems can be found in a large literature [29,32,34,35,38,40,49,58,99–102], and new experiments in chemical media are useful only if the material owns very specific properties, little is known about nonlinearity in biological fluids and the related literature is very recent [43,97]. Bio-materials can be very interesting, because both chemical and biological compounds can be excellent tunable thermal media, and DSWs were already observed [29,36,43].

For the sake of completeness, in this section, we report two experiments. The first one is in M-Cresol/Nylon, a chemical solution that exhibits an isotropic giant self-defocusing nonlocal nonlinearity, tunable by varying the nylon concentration [36]. The second one is in human red blood cell suspensions, where the concentration of hemoglobin (Hb) and the input laser beam power make the nonlinearity change from self-focusing to nonlocal defocusing [43].

Figure 12 shows transverse profiles of output beam intensity after a propagation of $2$mm in M-Cresol/Nylon. M-Cresol/Nylon is made up of an organic solvent (m-cresol) and a synthetic polymeric solute (nylon). When it is enlightened by a CW laser beam, light absorption induces local temperature variations, which reduces the refractive index, that is, the material experiences a nonlinear thermo-optical effect. In particular, authors' [58] measured the M-Cresol/Nylon nonlinear Kerr coefficient $n_2$ and found that, if for pure m-cresol it is $-9\times {10}^{-8}$cm${ }^2/$W, for a nylon mass concentration of $3.5\mathrm{\%}$ it is $-1.6\times {10}^{-5}$cm${ }^2/$W, an order of magnitude higher than most of the other thermal nonlinear materials reported in the literature. Authors generated the DSWs in Figure 12 by focusing the input beam (a CW laser of wavelength $532$nm) to $20\pm 1\mu$m onto the surface of M-Cresol/Nylon solution of $3.5\mathrm{\%}$ nylon concentration. The input laser power was varied ranging from $2\mu$W to $20$mW and, when it reached $5$mW, the wave-breaking occurred.

Figure 12. Output beam intensity transverse profiles, coming out from a $2$mm long M-Cresol/Nylon solution. Input power varies: (a) $P_{in}=2\mu$W, (b) $P_{in}=5$mW, (c) $P_{in}=10$mW, (d) $P_{in}=20$mW.

Reprinted with permission from [36]. Copyright 2014 Optical Society of America.

Figure 13 reports a part of the results obtained in lysed human red blood cells aged samples, where free Hb determines sign and nonlocality of the optical nonlinearity from self-focusing (and self-trapping) to strong thermal defocusing effects, regime in which DSWs occur [43]. Beyond the biological issues related to human red blood cells, holding uncountable applications to life sciences and medicine, red blood cells refractive index tunability makes this medium be incredibly interesting also from a physical point of view [103–106]. In normal conditions, red blood cells are disc-shaped malleable cells, averagely with $8\mu$m of diameter and $2\mu$m of thickness, which have a spatially uniform refractive index because of the lack of nuclei and most organelles [104,106]. To enable the passage through veins and narrow microcapillaries, red blood cells exhibit distinctive deformability. Since their optical properties depend on the shape and refractive index of cells, they can be used as tunable optofluidic microlenses [105].

Figure 13. Output beam waist for varying hemoglobin concentration and input power. (a) Detected beam diameter as function of input power through the hemoglobin solutions for four different concentrations (Hb1-Hb4): $2.4$, $5.1$, $8.6$, and $15.0$ million cells per mL. Nonlinear self-focusing of the beam occurs around $100$mW for high concentrations of hemoglobin, but it subsequently expands into thermal defocusing rings at high powers. (b-e) Output beam transverse intensity profiles for (b) self-trapped beam at high concentration and low power, (c) DSW at high concentration and high power, (d) self-trapped beam at low concentration and low power, (e) DSW at low concentration and high power.

Reprinted by permission from Macmillan Publishers Ltd. from [43]. Copyright 2019.

The red blood cell refractive index is mainly determined by Hb, which is the largest part of the erythrocyte dry content by weight [103]. Figure 13a shows the output beam waist as a function of input power through the Hb solutions for four different concentrations, from $2.4$ to $15.0$ million cells per mL. Experiments in [43] are performed by using a linearly polarized CW laser beam with a wavelength of $532$nm focused through a lens of $125$mm focal length into a $3$cm long glass cuvette filled with the red blood cell suspensions. In particular, the focused beam has initial waist $W_0=28\mu$m at the focal point, which was located at $1$cm away from the input facet of the cuvette to avoid heating and surface effects [97]. Outputs from the sample were monitored with a CCD camera and a power detector, and are reported in Figure 13(b–e), at variance of Hb concentration and input power. DSWs occur at high power (Figure 13(c,e)), more visible in high Hb concentration regime (Figure 13c).

3.4. Incoherent dispersive shock waves

In this section, we report experimental evidence of incoherent DSWs theoretically introduced in Section 2.3. In particular, we show observations of the transition from shocklets to collective incoherent DSWs [39].

By varying the effective range of nonlocality, authors of [39] performed experiments both in the quasi-local and in the highly nonlocal regime. The experimental setup is sketched in Figure 14. A CW laser ($\lambda =532$nm) is made incoherent by passing through a 4-f telescope with a ground-glass plate in the middle. The initial coherence length ${\lambda }_c^0\sim 200\mu$m, that is, the size of the speckles induced in the beam at the sample input facet, is controlled by changing the waist impinging on the ground-glass plate. The sample is a dilute solution of methanol and graphene nanoscale flakes. The latter provide optimal conversion of absorbed laser energy into heat, thanks to the absence of fluorescence mechanisms and a negligible absorption. A CCD camera detects intensity images at the output.

Figure 14. Experimental setup. A CW laser beam, $\lambda =532$nm, is sent through a 4-f telescope. A ground-glass plate G, placed in the midst of the telescope (on the focus of the first lens), generates a speckle pattern. The incoherent beam impinges the samples, a cylindrical tube filled with a solution of methanol and graphene nanoscale flakes, with waist $W_0=2.3$mm, while the initial coherence length ${\lambda }_c^0$ is controlled by changing the beam size on G. A CDD camera detects the output.

Reprinted by permission from Macmillan Publishers Ltd. from [39]. Copyright 2015.

Figure 15 shows the shocklets formation in a quasi-local regime. To reach this stage, concentration of graphene nano-flakes was increased, to increase the absorption and reduce the nonlocality, and a higher coherence length was chosen (${\lambda }_c^0\sim 250\mu$m), to inhibit collective behaviors. At low power (Figure 15a for experiments, Figure 15e for NLSE simulations) propagation is linear, so dominated by diffraction. As the input power increases (Figure 15b for experiments, Figure 15f for NLSE simulations), each speckle develops its own DSW, with a regular undular pattern. These are more evident in experimental (Figure 15(c,d)) and numerical (Figure 15(g,h)) zooms of shocklets undular bores.

Figure 15. Experimental observation at short-range regime. (a, b) Experimental beam profiles of the output intensity recorded at power a $P=0.05$W, (b) $P=2.50$W. (c, d) Zooms on details of (b) that evidence the development of shocklets. (e, f) Numerical simulations of NLSE equation, and (g, h) corresponding zooms.

Reprinted by permission from Macmillan Publishers Ltd. from [39]. Copyright 2015.

Figure 16 reports the emergence of a giant collective incoherent DSW in highly nonlocal regime. Such a transition was made by decreasing the concentration of graphene nano-flakes and the coherence length (${\lambda }_c^0\sim 200\mu$m). At low power, no collective behavior was observed, neither in the intensity profile (Figure 16(a,d) for experiments, Figure 16g for NLSE simulations), nor in the spectrogram (Figure 16j for experiments, Figure 16m for LRVE simulations). At nonlinear regime, the occurrence of the annular reshaping of the beam is visible, with high frequencies piled on the boundaries, and low frequencies dominating the central part (Figure 16(b,c,e,f) for experiments, Figure 16(h,i) for NLSE simulations). The corresponding spectrograms show the linear behavior (Figure 16(j,m)) at low power and Z-shape distortion at higher powers. The latter phenomenon is a signature of incoherent DSWs (see Section 2.3), and is observable in Figure 16(k–o), respectively in experiments and simulations.

Figure 16. Experimental observation at long-range regime. (a-c) One-dimensional intensity transverse profiles along $x$ at $y=0$. (d-f) Two-dimensional intensity profiles. The asymmetry in the lower part of the beam is due to convection within the sample. (g-i) Numerical simulations of NLSE. (j-l) Experimental and (m-o) numerical spectrograms: the Z-shaped distortion reveals a dramatic coherence degradation on the annular boundaries of the beam (the coherence length decreases at the shock front), while a significant coherence enhancement occurs in the internal region of the beam. Input beam power: (a, d, g, j, m) $P=0.05$W, (b, e, h, k, n) $P=1.25$W, (c, f, i, l, o) $P=2.50$W.

Reprinted by permission from Macmillan Publishers Ltd. from [39]. Copyright 2015.

All the shock phenomena reported here, as DSWs in general, do not arise in the weakly nonlinear (weak turbulence) regime, but solely in the strongly nonlinear (strong turbulence) regime. The Vlasov equation is valid beyond the weakly nonlinear regime and thus describes the collective incoherent shocks in the highly nonlocal case [39]. However, in the weakly nonlocal regime, there is no theory that describes the development of the shocklets reported in Figure 15: despite intense efforts from several decades, strong turbulence still constitutes a challenging-unsolved problem of classical physics.

4. Conclusions

We reviewed the most widespread current theoretical models that describe nonlocal NLSE DSWs in spatial optical beam propagation. Moreover, we discussed their experimental observations.

In Section 2 the derivation of nonlocal NLSE was detailed, and main features of wave breaking in thermal Kerr media were reported [29]. In order to exhibit the theoretical interpretations of these phenomena as intrinsically irreversible, TAQM and turbulence wave theory approaches were summarized [12,39,41].

Section 3 is a collection of experiments on DSW generation in thermal media, first about a quite rich literature on observations in Rhodamine [29], and their TAQM explanation [37,38,40]. As second instance, we analyzed the interaction between disorder and nonlinearity in Rhodamine with silica spheres [32] and in silica aerogel [35], where the randomness inhibits the DSWs occurrence. Moreover, we reviewed very recent works on generation of photonic wave breaking in chemical [36] and biological solutions [43], fields where DSWs are emerging as surprising tools, useful for sensing and control of extreme phenomena. Finally, we showed emergence of giant collective incoherent shock waves from random wave propagation in highly nonlocal media [39].

The study of nonlinear optics in new materials, like soft matter and biological suspensions, is opening the way to a new branch of photonics. Many applications include, but are not limited to, spectroscopy, medicine, life sciences, non-invasive diagnosis and time-resolved low-power probes. May the physics of nonlinear waves support the development of these new directions? Can novel mathematical tools deepen our understanding of nonlinear radiation–matter interaction? This manuscript is intended to sustain the improvement of theory and experiments concerning nonlinear optical propagation in highly nonlocal and complex matter.

We thank A. Picozzi, S. Trillo and G. Xu for the thoughtful discussions and suggestions on the manuscript, and MD Deen Islam for the assistance in the laboratory. We acknowledge support from the QuantERA ERA-NET Co-fund grant number 731473 (Project QUOMPLEX), H2020 project grant number 820392 (Project PHOQUS), Sapienza Ateneo 2016 and 2017, PRIN 2015 NEMO (2015KEZNYM), PRIN 2017 PELM (20177PSCKT), Joint Bilateral Scientific Cooperation CNR-Italy/RFBR-Russia 2018–2020, and the National Key R&D Program of China (2017YFA0303800).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the QuantERA ERA-NET Co-fund, project QUOMPLEX [731473]; H2020, project PHOQUS [820392]; National Key R&D Program of China [2017YFA0303800]; Sapienza Ateneo [2016, 2017]; Joint Bilateral Scientific Cooperation CNR-Italy/RFBR-Russia 2018-2020 [NA]; PRIN 2015 NEMO [2015KEZNYM]; PRIN 2017 PELM [20177PSCKT].