Optical Spatial Shock Waves in Nonlocal Nonlinear Media

TL;DR

This paper reviews recent theoretical models and experiments on optical spatial dispersive shock waves in nonlocal media, highlighting their unique properties and applications in soft matter and biological systems.

Contribution

It provides a comprehensive overview of the latest developments in modeling and experimental observation of spatial dispersive shock waves in nonlocal optical media.

Findings

Identification of unconventional properties of shock waves in nonlocal media

Recent experimental demonstrations in colloidal and biological materials

Insights into potential sensing and control applications

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.

Equations

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

$$\mathbf{P}=\epsilon_{0}\Bigr{(}\chi^{(1)}\mathbf{E}+\chi^{(2)}\mathbf{E}\mathbf{E}+\chi^{(3)}\mathbf{E}\mathbf{E}\mathbf{E}+\dots\Bigl{)}=\mathbf{P}^{(L)}+\mathbf{P}^{(NL)},$$

$$\mathbf{P}=\epsilon_{0}\Bigr{(}\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\Bigl{)}=\mathbf{P}^{(L)}+\mathbf{P}^{(NL)},$$

$${\chi }^{(n)}\ast \mathbf{E}\dots \mathbf{E}(t)={\int }_{-\infty }^t\mathrm{d}{t}_1{\int }_{-\infty }^t\mathrm{d}{t}_2\cdots {\int }_{-\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).$$

$${\mathbf{P}}^{(NL)}(\mathbf{R},t)={ϵ}_0{\int }_{-\infty }^t\mathrm{d}{t}_1{\int }_{-\infty }^t\mathrm{d}{t}_2{\int }_{-\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)$$

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

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

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

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

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

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

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

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

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

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

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

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

$$\begin{array}[]{rcl}\partial_{z}\rho+\left[\frac{D-1}{r}\rho u+\partial_{r}(\rho u)\right]&=&-\alpha\rho,\\ \partial_{z}u+u\partial_{r}u-\chi\partial_{r}\theta&=&0,\\ -\sigma^{2}\left[\partial_{r}^{2}\theta+\frac{D-1}{r}\partial_{r}\theta)\right]+\theta&=&\rho.\\ \end{array}$$

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

$$ı{\partial }_z\psi +\frac{1}{2}{\partial }_x^2\psi -\kappa \ast |\psi {|}^2\psi =-ı\frac{\alpha }{2}\psi .$$

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

$$\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.$$

$$\begin{array}[]{l}\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}$$

$${\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),$$

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

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

$${\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{{\Gamma }_n}{2}z}.$$

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

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

$$\begin{array}[]{ccc}B(\mathbf{r},\mathbf{\xi},z)=\langle\psi(\mathbf{r}+\frac{1}{2}\mathbf{\xi},z)\psi^{*}(\mathbf{r}-\frac{1}{2}\mathbf{\xi},z)\rangle,&\mathbf{r}=(\mathbf{r_{1}}+\mathbf{r_{2}})/2,&\mathbf{\xi}=\mathbf{r_{1}}-\mathbf{r_{2}}.\end{array}$$

Full text

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.

††articletype: REVIEW

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, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 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, 19, 20, 21, 22, 23], quantum liquids [24], photorefractive media [25], plasma physics [26], viscous fluids [27], and diffracting optical beams [28, 5, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 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 [46, 47, 48, 49, 50, 51, 52, 4, 29, 53, 54, 55, 56, 57, 34, 58, 59, 60, 11]. 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) [63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 12], which models the nonlinear wave intrinsically irreversible propagation as a superposition of decaying resonances, and the wave turbulence theory [74, 75, 33, 76, 39, 41, 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 Sec. 2, we report the main theoretical approaches and results related to DSWs. Section 2.1 explains in details 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 Secs. 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 [76, 39, 41].

The second part of the manuscript is a collection of experiments on DSW generation in thermal media. Sec. 3 reports observations in Rhodamine solutions [29]. Output beam intensity profiles in Sec. 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 Sec. 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 Sec. 3.3. Last but not least, Sec. 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

[TABLE]

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

[TABLE]

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]:

[TABLE]

where $\ast$ is the convolution product

[TABLE]

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

[TABLE]

and the related dielectric tensor changes as

[TABLE]

where $\langle\mathbf{E}\cdot\mathbf{E}\rangle=\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 in a perturbation of the medium refractive index, proportional to the field intensity:

[TABLE]

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]:

[TABLE]

with $\left(\frac{\partial n}{\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 [46, 79, 80, 52, 51, 29, 53, 59] with constant boundary conditions (at room temperature):

[TABLE]

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}_{\perp},Z)$. The solution can be written as

[TABLE]

with $G(\mathbf{R_{\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{\partial n}{\partial T}\right|_{0}}}\propto\sqrt{L_{loss}}$:

[TABLE]

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

[TABLE]

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

[TABLE]

In low absorption regime ($Z_{0}<<L_{loss}$) $\Delta T(\mathbf{R})\sim\Delta T_{\perp}(\mathbf{R_{\perp}})$ and $\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}+\Delta n[I](\mathbf{R_{\perp}})$, with the refractive index nonlocal perturbation

[TABLE]

and n_{2}K(\mathbf{R_{\perp}})=\Bigl{(}\frac{\partial n}{\partial T}\Bigr{)}_{0}G_{\perp}(\mathbf{R_{\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_{\perp}}-\mathbf{R_{\perp}}^{\prime})\sim\delta(\mathbf{R_{\perp}}-\mathbf{R_{\perp}}^{\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_{\perp}})P(Z)$, with $P(Z)=\int\mathrm{d}\mathbf{R_{\perp}}I(\mathbf{R})$ the power.

For a monocromatic field $\mathbf{E}(\mathbf{R},\tau)=\hat{\mathbf{E}}_{0}A(\mathbf{R})e^{-\imath\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]:

[TABLE]

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, 19, 20, 21, 29, 30, 26, 33, 27] 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 [63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 12] 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}=kW_{0}^{2}$ the diffraction length, one obtains the normalized nonlocal NLSE

[TABLE]

with $\epsilon=\frac{L_{nl}}{L}=\sqrt{\frac{L_{nl}}{L_{d}}}$ a small quantity in strongly nonlinear (or weakly diffracting) regime, as the one we are considering, $\chi=\frac{n_{2}}{|n_{2}|}$, $\theta=\left|\frac{k\Delta nL_{nl}}{n_{0}}\right|$, $\alpha=\frac{L}{L_{loss}}$. From Eq. (11)

[TABLE]

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 TEM00) 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[\imath\frac{\phi(r,z)}{\epsilon}\right]$. For $D=2$ the transverse dimensionality and $u=\partial_{r}\phi$ the phase chirp, from Eqs. (15, 16) one obtains

[TABLE]

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

[TABLE]

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 $\epsilon^{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\epsilon$, 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 Figs. 1b,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 Sec. 2.3. The analysis made for random optical waves in Sec. 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 Fig. 1d (and, in the following, at the intensity profiles in Fig. 4a, Fig. 5, Figs. 6D,E, Figs. 8a,b and Figs. 16b,c), and comparing this to the fast oscillations in Fig. 1c [or Fig. 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.

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 $\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}$:

[TABLE]

with $\kappa(x)=W_{0}\tilde{K}(xW_{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 [82, 55]

[TABLE]

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 $\imath\partial_{z}\psi=\hat{H}\psi$, with the Hamiltonian $\hat{H}=\tfrac{1}{2}\hat{p}^{2}+p\kappa(x)$ $\left(\hat{p}=-\imath\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]:

[TABLE]

If $\psi=\exp\left(-\imath\kappa_{0}^{2}pz\right)\phi$, then $\imath\partial_{z}\phi=\hat{H}_{RHO}\phi$.

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 a orthonormal basis (both a orthogonality and a completeness relations hold):

[TABLE]

with $H_{n}(x)=(-1)^{n}x^{2}\frac{d^{n}}{dx^{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, 83, 73] through the transformation $\omega\rightarrow\imath\gamma$, $\hat{x}\rightarrow e^{-\imath\frac{\pi}{4}}\hat{x}$, $\hat{p}\rightarrow e^{\imath\frac{\pi}{4}}\hat{p}$ [70, 73]. The new stationary Schrödinger equation is $\hat{H}_{RHO}\mathfrak{f}^{\pm}(x)=\imath\frac{\Gamma}{2}\mathfrak{f}^{\pm}(x)$, solved by the spectrum $\frac{\Gamma_{n}}{2}=\gamma\left(n+\frac{1}{2}\right)$ and the non normalizable eigenfunctions

[TABLE]

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]:

[TABLE]

with

[TABLE]

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

[TABLE]

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{\Gamma_{n}}\mathfrak{f}_{n}^{+}$ gives the probability $p_{n}(z)$ of finding the system in a decaying GV

[TABLE]

which gives the $z-$dependent weight of the $n$-order GV. Initial weights $p_{n}(0)$ are reported in Fig. 3c as functions of $\gamma^{2}$. Since a Gaussian beam $\psi(x,0)=\varphi(x)=\exp(-x^{2}/2)/\sqrt[4]{\pi}$ is an even input, all the odd terms in Eq. (26) vanish due 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 Figs. 4b,c. While dotted profiles are numerical results from Eq. (27), continuous lines result from the general projection definition $p_{n}(z)=\Gamma_{n}|\langle\mathfrak{f}_{n}^{+}|\psi(x,z)\rangle|^{2}$, with $\psi(x,z)$ numerical solution of Eq. (19).

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 [86, 87, 74, 88, 89, 75, 33, 76, 39, 41]. 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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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 $\epsilon_{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 expansion with respect to $\epsilon_{c}$:

[TABLE]

for the local spectrum $n_{\mathbf{k}}(\mathbf{r},z)=\int\mathrm{d}\mathbf{\xi}B(\mathbf{r},\mathbf{\xi},z)\exp({-\imath\mathbf{k}\cdot\mathbf{\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 trasform of the response function $\tilde{U}(\mathbf{k})=\int\mathrm{d}\mathbf{r}U(\mathbf{r},z)\exp({-\imath\mathbf{k}\cdot\mathbf{r}})$, for the LRTE

[TABLE]

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

[TABLE]

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 an 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 Sec. 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\Delta k^{-1}$ decreases at shock front ($\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.6$mM. Measurements of the shock profiles are in Fig. 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}$cm2W*-1*, loss length $L_{loss}^{-1}=62$cm*-1*. For water $D_{T}=1.5\times 10^{-7}$m2s*-1*, $\rho_{0}=10^{3}$kg m*-3*, $c_{p}=4\times 10^{3}$J kg*-1K-1*, $\left|\frac{\partial n}{\partial T}\right|_{0}=10^{-4}$K*-1*. The degree of nonlocality is estimated as $\sigma=0.3$.

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 $\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 Fig. 6A. Samples are prepared by dispersing $0.1$mM 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}$m2W*-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.

Figures 6B,C report 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 Figs. 6D,E. These are signatures of DSWs in nonlinear media at high power.

Decay rates in Fig. 7 are detected by slicing the intensity profile $I(X,Z)$ at $X\simeq 0.1$mm (yellow line in Fig. 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 obey the quantized spectrum scaling $\Gamma_{2}/\Gamma_{0}=5$ at all investigated power levels, as shown in Fig. 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 TEM00 $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.

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}\rightarrow\mathfrak{p}$ and $\hat{x}\rightarrow\imath\partial_{\mathfrak{p}}$)

[TABLE]

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:

[TABLE]

where

[TABLE]

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):

[TABLE]

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 forth 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 Figs. 8a,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 Figs. 8c,d.

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 Sec. 2.2.

Undular bores of DSWs were analyzed and exhibited in Fig. 8f, while the field intensity undulation period $T$ versus $P$ is shown in Fig. 8g. In order to demonstrate univocally that $T\propto\sqrt[4]{P}$, inset in Fig. 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 [91, 54, 56]. 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 combines 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.1$mM 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 $\Delta n[|A|^{2}]=n_{2}|A|^{2}+\Delta n_{R}(X,Y,Z)$ and $L_{loss}\sim\infty$ (no loss). Through the same scaling of Sec. 2.1 and approximation to cylindrical symmetry $\partial_{y}\sim 0$, we obtain

[TABLE]

with $U_{R}(x,y,z)=\frac{\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 $\epsilon\sim 0$, the phase chirp behaves like a moving unitary mass particle [32]:

[TABLE]

with $U=\exp\left(-x^{2}/2\right)$ the deterministic potential for a Gaussiam TEM00 given by the nonlinearity, and $\eta_{R}=-\frac{\mbox{d}U_{R}}{\mbox{d}x}$ a Langevin force with Gaussian distribution, such that $\langle\eta_{R}(z)\eta_{R}(z^{\prime})\rangle=\eta^{2}\delta(z-z^{\prime})$ and $\eta=\sqrt{\left\langle\left(\frac{\mbox{d}U_{R}}{\mbox{d}x}\right)^{2}\right\rangle}\simeq\sqrt{\langle\left(\Delta n_{R}\right)^{2}\rangle}(|n_{2}|I_{0})^{-1}$ the disorder strength. Brackets $\langle,\rangle$ 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)$ (Figs. 9a,b) and phase space $(x,v)$ (Figs. 9c,d), where $v=\frac{\mbox{d}x}{\mbox{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 (Figs. 9a,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 presence of disorder, Figs. 9b,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.

We report here the experiments in RhB with silica spheres 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 Fig. 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 $1$mm$\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], 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 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 Fig. 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.

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$/$cm3 and refractive index $n_{0}=1.074$. Experimental setup is very similar to the previous ones (Fig. 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\mbox{mm})$ by the CCD camera. Observations are shown in Fig. 11. Images in the second and third rows of Fig. 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.

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 [95, 96, 58, 36, 97, 98, 43], 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 [49, 99, 29, 100, 32, 34, 101, 102, 35, 58, 38, 40], 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 [97, 43]. 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 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, [58]’s authors 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\%$ 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 literature. Authors generated the DSWs in Fig. 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\%$ 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 occured.

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, 104, 105, 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].

The red blood cell refractive index is mainly determined by Hb, which is the largest part of the erythrocyte dry content by weight [103]. Fig. 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 Figs. 13b-e, at variance of Hb concentration and input power. DSWs occur at high power (Figs. 13c,e), more visible in high Hb concentration regime (Fig. 13c).

3.4 Incoherent Dispersive Shock Waves

In this section, we report experimental evidence of incoherent DSWs theoretically introduced in Sec. 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 Fig. 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 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 (Fig. 15a for experiments, Fig. 15e for NLSE simulations) propagation is linear, so dominated by diffraction. As the input power increases (Fig. 15b for experiments, Fig. 15f for NLSE simulations), each speckle develops its own DSW, with a regular undular pattern. These is more evident in experimental (Figs. 15c,d) and numerical (Figs. 15g,h) zooms of shocklets undular bores.

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 (Figs. 16a,d for experiments, Fig. 16g for NLSE simulations), nor in the spectrogram (Fig. 16j for experiments, Fig. 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 (Figs. 16b,c,e,f for experiments, Figs. 16h,i for NLSE simulations). The corresponding spectrograms show the linear behavior (Figs. 16j,m) at low power and Z-shape distortion at higher powers. The latter phenomenon is a signature of incoherent DSWs (see Sec. 2.3), and is observable in Figs. 16k,l and Figs. 16n,o, respectevely in experiments and simulations.

All the shock phenomena reported here, as DSWs in general, do not arise in the weakly nonlinear (weak turbolence) 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 Fig. 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 Sec. 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 explaination [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 thoughful 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 731473 (Project QUOMPLEX), H2020 project grant number 820392, Sapienza Ateneo, PRIN 2015 NEMO, PRIN 2017 PELM, Joint Bilateral Scientic Cooperation CNR-Italy/RFBR-Russia 2018-2020, and the National Key R&D Program of China (2017YFA0303800).