Optical Spatial Shock Waves in Nonlocal Nonlinear Media

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.

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 inten-sity. 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][47][48][49][50][51][52][53][54][55][56][57][58][59][60]. In these materials, light propagation is affected by a highly nonlocal Kerr nonlinearity, ruled by nonlocal NLSE.
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 [39,41,76].
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].

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 with E the electric field and P the medium nonlinear polarization [78]. The relation between P and E depends on the material properties. Including all the nonlinear terms, we have where 1 + χ (1) = n 2 0 , n 0 is the medium refractive index, χ (2) and χ (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 radiationmatter interaction is more properly represented by the following non instantaneous superposition of linear and nonlinear polarization [78]: where * is the convolution product If we have a third-order isotropic and centrosymmetric material, the nonlinear polarization is (4) and the related dielectric tensor changes as where E·E = 1 2 |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: with I = |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]: with ∂n ∂T 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): where γ = (L loss ρ 0 c P D T ) −1 , L loss is the loss characteristic length, ρ 0 is the material density, c P is the specific heat at constant pressure, D T is the thermal diffusivity and R = (X, Y, Z) = (R ⊥ , Z). The solution can be written as with G(R ⊥ ) 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 radiationmatter interaction: it works only in a neighborhood of the sample midpointẐ = Z 0 /2, not in proximity of the borders. Mathematically, this is deciphered in a longitudinal parabolic approximation with characteristic width L nloc = |n2| γ| ∂n From Eqs. (8,10) we obtain the 2D heat equation with In low absorption regime (Z 0 << L loss ) ∆T (R) ∼ ∆T ⊥ (R ⊥ ) and ∂ Z I(R) ∼ 0 (intensity longitudinal changes are negligible as for solitary wave packets), therefore we attain n[I](R) = n 0 + ∆n[I](R ⊥ ), with the refractive index nonlocal perturbation and n 2 K(R ⊥ ) = ∂n ∂T 0 G ⊥ (R ⊥ ). 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 W 0 (local approximation), i.e., K(R ⊥ − R ⊥ ) ∼ δ(R ⊥ − R ⊥ ), and the opposite case L nloc W 0 , that is, the highly nonlocal approximation (HNA), where K * I(R) ∼ K(R ⊥ )P (Z), with P (Z) = dR ⊥ I(R) the power.
For a monocromatic field E(R, τ ) =Ê 0 A(R)e −ıωτ in a third-order thermal medium, in paraxial and slowly varying envelope approximations, introducing the delayed time τ = t − n0 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]: with k = 2πn0 λ = ωn0 c the wavenumber.
Starting from Eq. (14), through the scaling , with I 0 the intensity peak, L = √ L nl L d , L nl = n0 k|n2|I0 the nonlinear length scale associated to a local Kerr effect, L d = kW 2 0 the diffraction length, one obtains the normalized nonlocal NLSE Ld a small quantity in strongly nonlinear (or weakly diffracting) regime, as the one we are considering, where σ = Lnloc W0 is the nonlocality degree, which expresses the nature of the Kerr effect through the limits we have previously discussed: if σ 1 we are considering the local limit θ ∼ |ψ| 2 , instead, if σ 1, by HNA θ ∼ κ(x, y)p(z), with p(z) = dxdy |ψ(x, y, z)| 2 = P (Z) The fundamental laser mode (Gaussian TEM 00 ) is described by an axisymmetric Gaussian input ψ 0 (r) = exp(−r 2 ), with r = x 2 + y 2 , which evolves in the WKB approximation [81] as ψ(r, z) = ρ(r, z) exp ı φ(r,z) . For D = 2 the transverse dimensionality and u = ∂ r φ the phase chirp, from Eqs. (15,16) one obtains Figure 1 reports phase chirp and field amplitude for D = 1, so for ∂ y ∼ 0 and r → x, in a defocusing medium (χ = −1) without losses (α = 0). The local case (σ = 0) is illustrated in Figs. 1a,c and follows from system (17): Eqs. (18) are equivalent to Euler and continuity equations, respectively, for a fluid of speed u, mass density ρ and pressure proportional to ρ 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 ∼ , which appear simultaneously in phase chirp u and intensity ρ. For large values of σ, the normalized refractive index variation, here expressed by θ(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 ρ 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 [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.

High Nonlocality and Time Asymmetric Quantum Mechanics
Let us consider the nonlocal NLSE in Eq. (14) with a medium response function , through the approximation ∂ Y ∼ 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 : with κ(x) = W 0K (xW 0 ) = exp (−|x|/σ) /(2σ). We take into account a medium where the nonlocality length is much larger than the beam waist. By HNA we have [55,82] where κ is a function no more depending on |ψ| 2 . In a system without loss, that is, α = 0, the normalized power p is conserved and the NLSE is mapped into a linear Schrödinger equation ı∂ z ψ =Ĥψ, with the HamiltonianĤ = 1 2p 2 + pκ(x) (p = −ı∂ x ). When we express the even function κ as its second order expansion, that is, we obtain the reversed harmonic oscillator (RHO) Hamiltonian [68,70,73]: If ψ = exp −ıκ 2 0 pz φ, then ı∂ z φ =Ĥ RHO φ. 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):Ĥ with H n (x) = (−1) n x 2 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,73,83] through the transformation ω → ıγ,x → e −ı π 4x ,p → e ı π 4p [70,73]. The new stationary Schrödinger equation isĤ RHO , solved by the spectrum Γn 2 = γ n + 1 2 and the non normalizable eigenfunctions 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]: Figure 3 shows the GV square norms (Fig. 3a) and phase chirps (Fig. 3b). The evolution of the normalized field ψ presents a Gamow part resulting as a superposition of exponential decays with quantized decay rates [73]: 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 √ Γ n f + n gives the probability p n (z) of finding the system in a decaying GV which gives the z−dependent weight of the n-order GV. Initial weights p n (0) are reported in Fig. 3c as functions of γ 2 . Since a Gaussian beam ψ(x, 0) = ϕ(x) = exp(−x 2 /2)/ 4 √ π 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

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][75][76][86][87][88][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 longrange 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), 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 The statistics is said to be homogeneous if B depends only on the distance |r 1 − r 2 |. Following Eq. (28) with P (r, ξ, z) = B(r, ξ, z) dr U (r ) N (r − r + 1 2 ξ, z) − N (r − r − 1 2 ξ, z) , In the last equation, N is the averaged field power, also depending on r because of the inhomogeneity of the statistics. By defining the length scale of random fluctuations as the coherence length λ c and W 0 the incoherent beam waist, we can assume that the statistics is quasi-homogeneous if c = λ c /W 0 << 1. For a Gaussian response function U (r) = 1 2πσ 2 exp − |r| 2 2σ 2 , we get the SRVE if σ << 1 and the LRVE if σ >> 1 through two different multiscale expansion with respect to c : for the local spectrum n k (r, z) = dξB(r, ξ, z) exp(−ık · ξ), with averaged power N (r, z) = 1 (2π) 2 dkn k (r, z). The nonlocal features of Eq. (32) are traced by the generalized dispersion relation. Once defined the Fourier trasform of the response function Ũ (k) = drU (r, z) exp(−ık · r), for the LRTẼ One can also compute the momentum both from the NLSE and the LRVE: 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 λ c ∼ ∆k −1 decreases at shock front (∆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.

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 1mm filled with an aqueous solution of rhodamine B (RhB), with a concentration of 0.6mM. Measurements of the shock profiles are in Fig. 5. A Gaussian CW laser beam of intensity waist W 0 = 20µm, at wavelength λ = 532nm, propagates in a material with linear refractive index n 0 = 1.3, defocusing Kerr coefficient n 2 = −7 × 10 −7 cm 2 W −1 , loss length L −1 loss = 62cm −1 . For water D T = 1.5 × 10 −7 m 2 s −1 , ρ 0 = 10 3 kg m −3 , c p = 4 × 10 3 J kg −1 K −1 , ∂n ∂T 0 = 10 −4 K −1 . The degree of nonlocality is estimated as σ = 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.

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 Γ 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.1mM of RhB in water. The solution is placed in a cuvette 1mm thick in the propagation direction. The measured defocusing Kerr coefficient is |n 2 | = 2×10 −12 m 2 W −1 and the absorption length is L loss 1.6mm at the laser wavelength 532nm [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 10cm, and a focus spot size of 10µ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 = 20cm 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.0mm in the Y direction and 35µ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.0mm. 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 Mshaped 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 0.1mm (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 Γ 2 /Γ 0 = 5 at all investigated power levels, as shown in Fig. 7D. This demonstrates that authors of [38] excited the fundamental state f − 0 and the first excited state 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.  [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.9mm (blue). Reprinted by permission from Macmillan Publishers Ltd. from [38]. Copyright 2015.  Fig. 6B). B Numerically calculated decays in the conditions of panel A. C Peak region of the experimental curve at P = 450mW. 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, Γ 0 (filled circles) and the excited state, Γ 2 , (triangles). Reprinted by permission from Macmillan Publishers Ltd. from [38]. Copyright 2015.
tum basis (p → p andx → ı∂ p ) 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: where with W is the finite size of the physical system. During the evolution, the Gamow ground state has the lowest decay rate, i.e., γ/2. This allows to consider, in the long term evolution, only the fundamental GV, and so only the Fourier transform F of the fundamental state of Eq. (36): 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 γ, 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 γ, being γ = p √ πσ 2 . The Gaussian beam profile in the Y direction, that is, p ∝ exp(−y 2 ), provides the link between Y , P and γ. 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 ∝ 4 √ P , inset in Fig. 8g reports the period T as function of 4 √ P . The resulting linear behavior confirms the theoretical results.

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 disorderinduced 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 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µm diameter in 0.1mM aqueous solution of RhB. The second one is a 1mm×1mm×8.5mm parallelepiped of silica aerogel. Despite observations of DSWs √ 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.
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 ∆n[|A| 2 ] = n 2 |A| 2 + ∆n R (X, Y, Z) and L loss ∼ ∞ (no loss). Through the same scaling of Sec. 2.1 and approximation to cylindrical symmetry ∂ y ∼ 0, we obtain with U R (x, y, z) = ∆nR(X,Y,Z) n2I0 taken as a random dielectric noise mainly acting on the phase [32]. In the hydrodynamic limit ∼ 0, the phase chirp behaves like a moving unitary mass particle [32]: with U = exp −x 2 /2 the deterministic potential for a Gaussiam TEM 00 given by the nonlinearity, and η R = − dUR dx a Langevin force with Gaussian distribution, such that η R (z)η R (z ) = η 2 δ(z − z ) and η = dUR dx 2 (∆n R ) 2 (|n 2 |I 0 ) −1 the disorder strength. Brackets , denote the statistical average, and the dependence of η R on x, y is neglected for stochastic independence and cylindrical symmetry, respectively, thus η R η R (z). Figure 9 shows trajectories x(z) (Figs. 9a,b) and phase space (x, v) (Figs. 9c,d), where v = dx dz , respectively without (η = 0) and with (η = 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.005w/w to 0.03w/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 10µm. The aqueous solutions are put in 1mm×1cm×3cm glass cells with propagation along the 1mm 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.2mm to 1.6mm (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. Figure 10. Transverse intensity patterns for different input power P and silica spheres concentration c: (a) P = 5mW , c = 0w/w, (b)P = 400mW , c = 0w/w, (c) P = 5mW , c = 0.017w/w, (d) P = 400mW , c = 0.017w/w, (e) P = 5mW , c = 0.030w/w, (f) P = 400mW , c = 0.030w/w. White 1D curves show the measured section of the intensity profiles vs X. FIG. 1 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 indepth details are given in [35]. It turns out that the sample used in the experiment has mass density ρ = 0.215g/cm 3 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.5mm) 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. Figure 11. Far field intensity profiles at the output of the silica aerogel for P in ranging from 1mW to 1W, and input beam waist w 0 ranging from 43µm to 1.4mm. 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.

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][96][97][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][100][101][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 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 2mm 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 × 10 −8 cm 2 /W, for a nylon mass concentration of 3.5% it is −1.6 × 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 532nm) to 20 ± 1µm onto the surface of M-Cresol/Nylon solution of 3.5% nylon concentration. The input laser power was varied ranging from 2µW to 20mW and, when it reached 5mW, the wave-breaking occured. 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][104][105][106]. In normal conditions, red blood cells are disc-shaped malleable cells, averagely with 8µm of diameter and 2µ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 532nm focused through a lens of 125mm focal length into a 3cm long glass cuvette filled with the red blood cell suspensions. In particular, the focused beam has initial waist W 0 = 28µm at the focal point, which was located at 1cm 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).

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 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. Fig. 14. A CW laser (λ = 532nm) is made incoherent by passing through a 4-f telescope with a ground-glass plate in the middle. The initial coherence length λ 0 c ∼ 200µ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, λ = 532nm, is sent through a 4-f telescope. A groundglass 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.3mm, while the initial coherence length λ 0 c 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 (λ 0 c ∼ 250µ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 15. Experimental observation at short-range regime. a,b Experimental beam profiles of the output intensity recorded at power a P = 0.05W, b P = 2.50W. 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 (λ 0 c ∼ 200µ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, 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. 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.05W, (b,e,h,k,n) P = 1.25W, (c,f,i,l,o) P = 2.50W. Reprinted by permission from Macmillan Publishers Ltd. from [39]. Copyright 2015.

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.