From coherent shocklets to giant collective incoherent shock waves in nonlocal turbulent flows

Understanding turbulent flows arising from random dispersive waves that interact strongly through nonlinearities is a challenging issue in physics. Here we report the observation of a characteristic transition: strengthening the nonlocal character of the nonlinear response drives the system from a fully turbulent regime, featuring a sea of coherent small-scale dispersive shock waves (shocklets) towards the unexpected emergence of a giant collective incoherent shock wave. The front of such global incoherent shock carries most of the stochastic fluctuations and is responsible for a peculiar folding of the local spectrum. Nonlinear optics experiments performed in a solution of graphene nano-flakes clearly highlight this remarkable transition. Our observations shed new light on the role of long-range interactions in strongly nonlinear wave systems operating far from thermodynamic equilibrium, which reveals analogies with, for example, gravitational systems, and establishes a new scenario that can be common to many turbulent flows in photonic quantum fluids, hydrodynamics and Bose–Einstein condensates.

T he statistical description of nonequilibrium behaviour of random dispersive waves is well developed in the weak nonlinear limit, for which wave turbulence theory provides a powerful tool to interpret observations arising in contexts as different as ocean waves, quantum fluids, plasmas and nonlinear optics to name a few [1][2][3][4][5][6][7] . However, when nonlinearities are strong, such an approach breaks down and no general theory exists. Therefore, novel scenarios that can be predicted starting from universal models and observed in real systems are of paramount importance. In this context, we explore how shock waves (at variance with other coherent structures such as vortices, quasi-solitons, collapsing wavepackets or rogue waves 1-3,6-16 ) affect the turbulent flow of a conservative system of random optical waves (that is, a photon fluid) interacting through defocusing (or repulsive) nonlinearities. As a distinguishing feature, we consider the nonlocal character of the nonlinearity, which is indeed common to a large variety of nonlinear wave systems [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] . In such defocusing media, the generic signature of the strongly nonlinear regime is the formation of expanding undulatory structures known as dispersive shock waves (DSWs), which originate from the dispersion (diffraction) acting over steep fronts 33 developing via the nonlinearity. These fascinating DSW structures 34 have been observed from the small scales as in optics 24,[35][36][37][38] or in Bose-Einstein condensates [39][40][41][42][43] to intermediate hydrodynamic scales [44][45][46] , or even long astrophysical spatial scales 47 .
In this paper, we report a remarkable transition of the turbulent flow that occurs in our photon fluid system when the characteristic range of nonlocality is increased. In the quasilocal regime, the field evolution is ruled by the stochastic formation of small-scale DSW structures that we naturally denote as dispersive 'shocklets' (the term shocklet was introduced to designate sporadic steep fronts reported in high-speed compressible turbulence 48 and as flank formations associated with planetary-scale shocks in space 49 ). In this 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. In marked contrast with such regime, a remarkable self-organized regime emerges when the nonlinearity becomes highly nonlocal (long-range interaction). In this case, the turbulent flow follows an unexpected global collective behaviour that manifests itself in the formation of a giant shock singularity that emerges from the fluctuating field as a whole. Such a phenomenon is characterized by a strong non-homogeneous redistribution of the spatial fluctuations, whose description is provided in terms of a hydrodynamic-like model derived from singular solutions of a long-range Vlasov equation (LRVE). Our analysis reveals that (i) in marked contrast with coherent DSWs in conservative (Hamiltonian) systems, the regularization of the global incoherent shock does not require the formation of a regular oscillating DSW structure; and (ii) the self-organization of the turbulent wave ensemble into such a global shock structure is intimately related to the long-range character of the interaction, thus revealing interesting links with gravitationallike systems 50,51 . We emphasize that this regime is fundamentally different from other forms of turbulence that are dominated by inverse or direct cascades 1,2,[6][7][8]12 , vortex dynamics 13 , acoustic turbulence 52,53 or from the incoherent undulatory shocks that manifest solely in frequency domain in the weak Langmuir turbulent regime 54 . We directly observe the emergence of the giant collective incoherent shock phenomenon in an optical experiment that involves a graphene-based thermal medium, which allows for a tunable degree of nonlocality in the system.

Results
Nonlinear Schrödinger model. A nonlocal nonlinearity is found in several systems (dipolar Bose-Einstein condensates 26 , roton excitations in superfluids 28 , atomic vapours 23 , nematic liquid crystals 12,18 , thermal media 24 , glasses [19][20][21]25 , plasmas 27 or plasmonic metamaterials 32 ) and can be modelled by the following universal two-dimensional (2D) nonlinear Schrödinger equation (NLSE): where the dynamics occurs in the transverse plane r ¼ (x, y) and z plays the role of time. The parameters b ¼ 1=k L and g refer to the linear (dispersive) and nonlinear coefficients, k L being the laser wave number. We denote by s the spatial extension of the nonlocal response function U(r), that is, the range of nonlocal interaction. Equation (1) conserves two important quantities during the propagation, the 'power' N ¼ R j c j 2 dr, and the 'energy' (Hamiltonian) H ¼ E þ U ¼ const, where EðzÞ and UðzÞ denote the evolutions of the linear and nonlinear energy contributions (see Methods section). The dynamics is ruled by the comparison of s with the coherence length of the field l c , and the 'healing length' L ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi b=ð2 j g j rÞ p , where r is the intensity. The healing length denotes the typical length scale for which linear and nonlinear effects are of the same order, for example, the typical size of a soliton or a vortex. In the following, we focus the presentation on the more interesting 2D defocusing regime, g40, while the analysis can easily be transposed to the focusing regime or to one spatial dimension (Supplementary Note 3).
Local versus nonlocal regimes. The formation of shock waves is known to require a strong nonlinear interaction 24,[33][34][35][36][37][38][39][40][41][42][43][44][45] , that is, the initial random wave, c 0 ðrÞ, is such that U 0 ) E 0 , where U 0 and E 0 refer to the initial values of the energies at z ¼ 0 (Fig. 1a). Note that this condition is analogous to l 0 c ) L, l 0 c being the initial coherence length of the random wave, c 0 ðrÞ. Considering this strong nonlinear regime, we now contrast the case of a quasi-local (short-range) interaction, s $ L, with a highly nonlocal (longrange) interaction, s ) L. For s $ L, the incoherent wave leads to the formation of several coherent DSWs, which develop within each individual fluctuation of c 0 ðrÞ as shown in Fig. 1c,d,f. Since the range of the nonlocal response is smaller than the coherence length s $ L ( l 0 c À Á , every individual fluctuation evolves independently of each other. Hence, in this quasi-local (strong) turbulent regime, the incoherent wave develops singularities that are in essence coherent DSWs. These dispersive shocklets can be regarded as the conservative counterpart of viscous shocklets considered in high-speed turbulent flows 48 . Note that the development of a sea of DSWs manifests itself by a spectral filamentation process in frequency k-space (Fig. 1e). This is associated with the strong interaction among DSWs that emanate from different fluctuations and whose compression against each other leads to polygon-like patterns featured by effective one-dimensional (1D) DSW sides (Fig. 1d,f). Also note that although this shocklets regime occurs in the presence of wave dispersion (oðkÞ ¼ bk 2 from the linearized NLSE), it exhibits some interesting connections with a long-standing challenging issue of weakly dispersive acoustic-like wave turbulence 52,53 (see Discussion section and Supplementary Note 4).
This physical picture changes in a dramatic way in the highly nonlocal regime, s ) L. This is illustrated in Fig. 1g-l, which shows the evolution of the field starting from the same initial condition c 0 ðrÞ ð Þ as in the quasi-local regime. In the highly nonlocal regime, the fluctuations of the incoherent wave exhibit a global collective behaviour, which is responsible for the formation of a large-scale incoherent shock wave of a fundamental different nature. It is now the incoherent wave as a whole that develops a shock: the momentum of the speckled beam is radially outgoing (p NLS ðr; zÞ $ Iðc Ã rcÞ oriented along r as illustrated by the arrows in Fig. 1k) and exhibits a shock-like singularity, while the envelope of the intensity of the beam experiences an annular (ring-shaped) collapse-like behaviour ( Fig. 1g-l). 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 (l c decreases with z), while its internal region exhibits a significant coherence enhancement (l c increases with z) (Fig. 1g-l).
We study this peculiar phenomenon through the analysis of the 'local spectrum' (or spectrogram) of the random wave, that is, a spectrum that depends on the spatial position (r) because the random wave is characterized by fluctuations that are inhomogeneous in space. As illustrated in Fig. 1m-o, the evolution of the spectrogram exhibits a peculiar Z-shaped distortion, which contrasts with the regular deformation observed for E 0 $ U 0 (Supplementary Note 1 and Supplementary Fig. 2). In the strong nonlinear regime, U 0 ) E 0 , the spectrogram exhibits a pronounced self-steepening that leads, nearby the overtaking point, to a dramatic spectral broadening on the annular boundary of the incoherent beam, where the coherence length decreases down to the healing length, l c $ L ðE $ UÞ.
Here the increase in the linear energy E up to U shown in Fig. 1p is due indeed to such peculiar process of coherence degradation. This is in marked contrast to coherent DSWs in Hamiltonian systems in which such a balance between linear and nonlinear energies stems from the well-known formation of regular DSW oscillations 24,[33][34][35][36][38][39][40][41][42][43][44][45] .
Singular solutions. We describe theoretically the phenomenon of large-scale incoherent shock within the general framework of the Starting from the same initial condition (a), in the highly nonlocal regime s ¼ 100L, the random wave as a whole develops a giant collective incoherent shock (j-l); the corresponding intensity lineouts y ¼ 0 (g-i) indicate an annular collapse-like behaviour.
(m-o) Corresponding spectrogram evolutions of the incoherent shock: the Z-shaped distortion reveals a dramatic coherence degradation on the annular boundaries of the beam (the correlation length decreases at the shock front, l c;shock $ 1=Dk shock ), while a significant coherence enhancement occurs in the internal region of the beam (l c;inner $ 1=Dk inner increases). (p) Evolutions during the propagation of the linear (E, green) and nonlinear (U, red) contributions to the total energy (constant Hamiltonian, H, dashed blue): The incoherent shock develops in the strongly nonlinear regime U ) E ð Þ . The propagation length in the sample, z, is in units of the nonlinear length, L 0 ¼ 1=ðgrÞ.
wave turbulence formalism. The wave turbulence theory has been shown to provide a natural asymptotic closure of the hierarchy of moment equations for a system of weakly nonlinear dispersive waves [3][4][5] . On the basis of a generalized wave turbulence approach of the problem, we show that the appropriate statistical description of the large-scale incoherent shocks is provided by a collisionless LRVE 14,55 . Note in this respect that resonant fourwave interactions described by a collision term are still present in the long-term evolution of the system, although in the long-range regime their influence on the dynamics results of higher order with respect to the LRVE dynamics, see Methods section.
The long-range Vlasov formalism discussed here differs from the traditional Vlasov equation describing random waves in hydrodynamics 8,15,16 , plasmas 27 or optics, such as, for example, incoherent modulational instabilities 14,56 , or the collective clustering of incoherent solitons with inertial nonlinearities 57 . On the other hand, the structure of the LRVE is analogous to that describing systems of particles with long-range, for example, gravitational, interactions 50,51 . Contrary to conventional weak-turbulence approaches 1,2,6,7,58 , we show that the LRVE is valid beyond the weakly nonlinear regime of interaction: Simulations of incoherent shocks in the strong nonlinear regime U 0 ) E 0 ð Þ reveal a quantitative agreement between LRVE and NLSE models without adjustable parameters ( Fig. 2 and Supplementary Fig. 6).
We show in the Supplementary Note 2 that, if the initial condition is strongly nonlinear U 0 ) E 0 ð Þ , the momentum of the incoherent wave results to be radially outgoing, so that the 2D LRVE remarkably reduces to an effective 1D equation: where ň k (r, z) is the local averaged spectrum of the field at radial position r ¼ |r|, Vðr; zÞ ¼ g ð2pÞ 2 R 1 0Ũ ðr; r 0 ÞŇðr 0 ; zÞdr 0 , withŨ ðr; r 0 Þ ¼ R 2p 0 Uððr 2 þ r 02 À 2rr 0 cosyÞ 1 2 Þdy an effective radial response function andŇðr; zÞ ¼ R 1 0ň k ðr; zÞdk ¼ rNðr; zÞ, see Methods section. A simulation of equation (2) is reported in Fig. 2a-d, which confirms the peculiar folding of the spectrogram discussed above through NLSE simulations in Fig. 1m-o. Physical insight into the strongly nonlinear regime can be obtained through solutions of equation (2) with extremely narrow spectral distribution, ň k (r,z) ¼ Ň (r, z)d(k-K(r, z)), where the 'particle density' Ň (r, z) and 'momentum' K(r, z) satisfy the following hydrodynamic-like model.
We remark that, from a broader perspective, the general class of singular solution considered here are sometimes called 'mono-kinetic' (or single speed), in the sense that to each spatial position (r, z) corresponds a unique well-defined spectral (velocity) component (k, z) (ref. 51). This special form of singular solutions introduced in ref. 59 are important in that they make a direct correspondence with a generic hydrodynamic formalism, a property that found fundamental implications in a large variety of long-range interacting systems, for example, self-gravitating systems, 2D geophysical fluids, collisionless plasma or more generally wave-particle interactions 50,51 . Numerical simulations of equations (3) and (4) are found in quantitative agreement with those of NLSE and LRVE equations, without using any adjustable parameter, as remarkably illustrated in Fig. 2e-n. Starting from K(r,z ¼ 0) ¼ 0, the 'spectrogram' K(r, z) is first driven by the last nonlinear term in equation (3), while the Burgers-like (second) term of equation (3) subsequently leads to the gradient catastrophe of K(r, z). Quite remarkably, the finite 'time' (distance, z) shock singularity of K(r, z) is responsible The annular collapse singularity, max r (Ň (r, z)) (i), and shock singularity, max r (@ r K(r, z)) (n), are in agreement with the power law predicted by the theory B(z N À z) À 1 , where z N denotes the blow-up propagation length, see equation (5) (dark-dashed). The propagation length z is in units of L 0 . The NLSE simulation in e-n refers to the radial averaging of the results reported in Fig. 1j-l. The quantitative agreement between NLSE, LRVE and hydrodynamic models in e-n is obtained without using adjustable parameters (see Methods section).
for an annular blow-up singularity of the intensity envelope Ň (r, z) on the circular boundary of the incoherent beam. These singular behaviours can be described theoretically by solving equations (3) and (4) by the method of the characteristics, which shows that the gradient of the spectrogram and the intensity envelope along a characteristic R(z) exhibit a blow-up singularity at the finite 'time' z N (see Methods section): The tendency to diverge as Bz À 1 occurs regardless of the dimensionality of the problem that is, also in 1D (Supplementary Note 3 and Supplementary Figs 5 and 6), and has been confirmed by the simulations shown in Fig. 2i,n. Such dynamics becomes regularized by the NLSE and LRVE models, although an accurate description of the underlying regularization mechanism is hindered by a long-standing mathematical problem, namely, achieving a closure of the infinite hierarchy of k-moments equation for kinetic theories, see Methods section(Supplementary Note 3 and Supplementary Fig. 7). We emphasize that the 'hydrodynamic' model equations (3) and (4) recovers the 1D shallow-water equations under the substitution V(r, z)-Ň (r, z) (refs 33,40). Note that, in the optical context, this limit is relevant to the description of incoherent wave propagation in inertial nonlinear media 14,56,57 . However, shallow-water equations are hyperbolic equations that are known to not exhibit collapse singularities: the annular collapse singularity of the field intensity in equations (3) and (4) originates into the nontrivial nonlocal coupling between the effective potential,Ũðr; r 0 Þ and the 'particle density', Ň (r, z) (see the expression of V(r, z) below equation (2)). Note that, in one spatial dimension, such a coupling reduces to a convolution, and the system still develops finite-time collapse and shock singularities (Supplementary Note 3 and Supplementary Figs 5 and 6).
Experiments. We performed a series of experiments that provide evidence of both regimes of incoherent shocks and dispersive shocklets by varying the effective range of nonlocality in the nonlinear material. As illustrated in the experimental set-up in Fig. 3, the beam from a continuous wave laser (532 nm) is sent through a 4-f telescope. A ground-glass plate is placed in the focus of the first lens of the telescope, which creates a speckle pattern at the input of the sample. Notice that this initial speckle beam does not need to satisfy specific (for example, Gaussian) statistics because the incoherent shock phenomenon occurs irrespective of the wave statistics-a property justified by the fact that the LRVE formalism accurately describes the strongly nonlinear regime of interaction. The beam radius in the absence of the ground plate, and at the sample input window is 2.3 mm. The size of the speckles (that is, the coherence length of the random speckle pattern) at the sample input can be controlled by changing the beam size on the ground-glass plate, for example, by changing the beam size at the telescope input with an iris or, more simply, by shifting the position of the ground-glass plate along the beam propagation axis (Supplementary Note 1 and Supplementary Fig. 4). After the sample, the beam is imaged onto a CCD camera.
The liquid in the sample exhibits a thermal nonlocal nonlinearity that originates from a laser-induced heating of the liquid. The nonlinearity is intrinsically nonlocal due to heat diffusion and usually defocusing as a result of the negative thermo-optic coefficient (dn/dTo0). The nonlinear coefficient can be enhanced by including a highly absorbing dye (rhodamine 24 or iodine). In our experiments, we chose a dilute solution of graphene nanoscale flakes that provides optimal conversion of absorbed laser energy into heat because of the absence of any fluoresence mechanisms. The medium exhibits very low absorption, which entails a highly nonlocal (long range) interaction, a property that has been exploited to study the incoherent shock phenomenon. The measured absorption coefficient is a ¼ 0.013 cm À 1 , so that the nonlocal length is in the range sB400 mm (ref. 24), an estimate that should be considered as lower bound, since in extended media, such as ours, the nonlocal length may be significantly larger, sB900 mm (Supplementary Note 1; ref. 22). Note that, despite the presence of absorption, the nonlinear system behaves essentially as a conservative (Hamiltonian) system, since the characteristic absorption length (C76.9 cm) is much larger than the sample length (15 cm) and the nonlinear length (C0.04 cm). The corresponding measured refractive index change is DnC À 2 Â 10 À 4 at maximum input power (2.5 W) of the laser, which gives LB5 mm. Considering typical values of the correlation length of the speckle beam, l 0 c $ 200mm, we obtain the scale separation required for the observation of incoherent shocks, s4l 0 c ) L (see Methods section and Supplementary Note 1 for more details).
Figure 4a-f reports typically measured output intensity patterns of the optical beam. As the input power is increased, we clearly see an annular reshaping of the speckled beam, with high-frequency components piling up on the boundaries, whereas low-frequency components dominate the internal region of the beam. This is confirmed by the corresponding y ¼ 0 lineouts in Fig. 4a-c, which evidence the formation of two lateral high-intensity peaks, as predicted by the collapse behaviour. The corresponding spectrogram measurements reported in Fig. 4j-l were obtained by measuring the far-field of the optical beam for several different x-positions. The singular distortion of the spectrograms provides a clear experimental signature of the phenomenon of incoherent shock. We also performed experiments in other liquids (for example, ethanol) and obtained similar results. Moreover, reducing the initial coherence length l 0 c À Á leads to a substantial different dynamics with no evidence of incoherent shock formation ( Supplementary Figs 1 and 2).
The simulations of NLSE (1) have been performed by following the model developed in ref. 24, in which an effective 2D reduction of the three-dimensional (3D) heat equation is used to describe the thermal nonlinearity (Supplementary Note 1). Several works highlighted the role of remote boundary conditions of ARTICLE the thermal nonlinear sample in the dynamics of soliton propagation by considering the complete 3D heat equation [19][20][21] . In this respect, despite the fact that the nonlinear response function is in principle 'infinite range', the analysis revealed that in actual experiments the response function can be accurately characterized by an effective well-defined finite range 22 . The 2D reduction of the heat equation considered here precisely accounts for such an effective finite nonlocal range, s (refs 18,24). In this way, we model the experiment by the universal form of the NLSE (1), in which the nonlocal response function does not depend on the propagation (time) variable, z, in the present case, U(r) ¼ K 0 (r/s)/(2ps 2 ), where K 0 ( Á ) refers to the modified Bessel function of the second kind (Supplementary Note 1). It is interesting to note that, in spite of the simplicity of the 2D NLSE model (1) (for example, it does not account for heat-induced convection of the liquid, which is responsible for the up-down asymmetry in Figs 4 and 5), this model captures the essential phenomenological behaviour observed experimentally, as revealed by the qualitative agreement between the simulated intensity patterns (Fig. 4g-i and Fig. 5e-h) and spectrograms (Fig. 4m-o) with the experimental results.
In order to demonstrate that the emergence of such collective incoherent structure is ultimately linked to the long-range regime of interaction, we have made a second series of experiments to investigate the development of dispersive shocklets from a speckled beam. To this end, we significantly increased the concentration of graphene nano-flakes so as to increase the absorption (aC1.7 cm À 1 ), and thus reduce the nonlocal range of interaction by one order of magnitude, see Methods section. In order to inhibit long-range collective effects, the correlation length of the input speckle beam has to be chosen much larger than the nonlocal range, l 0 c ð $ 250mmÞ ) sð $ 70mmÞ. As illustrated in Fig. 5a-d, each individual speckle of the incoherent beam develops its own DSW, which leads to the formation of regular undular patterns, whose typical spatial period has been found in good agreement with the corresponding NLSE (1) simulations. As expected from the previous analysis, experiments realized with a smaller correlation length lead to an effective global behaviour characterized by a strong interaction among the different speckles, which strongly reduces the development of DSWs within each individual speckle of the input beam ( Supplementary Fig. 3). For more details on the experimental results, in both short-and long-range regimes, and the corresponding NLSE simulations, see Supplementary Note 1.

Discussion
We have reported a novel scenario of strongly nonlinear turbulent flows characterized by the emergence of a large-scale incoherent shock, which is inherently a collective phenomenon of the turbulent field as a whole. This strongly nonlinear turbulent regime is of a fundamental different nature than other forms of turbulence, which are dominated by nonequilibrium stationary cascades 1,2,6-8,12 , or deeply affected by localized nonlinear Spectrograms Collective incoherent shock waves   (Fig. 2e-n).
structures, such as solitonic turbulence 7,12 , or entangled vortices in superfluid turbulence 13 . By varying the range of the nonlocal thermal nonlinearity over one order of magnitude, a remarkable qualitative agreement has been obtained between the experimental results and the simulations of the NLSE model based on a 2D reduction of the heat equation. Despite such a qualitative agreement, it would be interesting to generalize the theory by considering the whole 3D Poisson equation describing heat diffusion also along the optical propagation axis 20,21 . This would lead to the unusual and interesting situation in which the long-range potential in the Vlasov formalism would depend itself on the 'temporal' z variable, a property that would introduce intriguing unexplored collective behaviours of the turbulent flow.
We also briefly comment an interesting connection between the coherent shocklets regime reported here and a long-standing challenging issue inherent to weakly dispersive turbulent systems 52,53 . At variance with strongly dispersive wave systems that can be accurately described by the wave turbulence theory 1-7 , a proper kinetic formulation of random nonlinear waves that exhibit weak acoustic-like dispersion constitutes a difficult problem that has been the subject of important developments in the past 1,2,52,53 . The first serious attempt of an analytical description of a 'semi-dispersive' wave turbulence was reported in ref. 52. Basically, the theory accurately describes how energy is shared along 'rays' in the frequency domain, that is, between wave vectors with the same orientation in k-space; however, it does not describe how energy is redistributed between distinct neighbouring rays. The foundation for an evaluation of the contribution to energy exchange that occurs at higher orders was formulated theoretically in ref. 53. This leaves unanswered an open challenging issue 2,53 : Given an initial anisotropic spectral distribution of the weakly dispersive random wave, one may wonder whether the nonlinear interactions of the next order lead to a spectral isotropic redistribution or to condensation along specific rays in k-space. In the light of this problem, we have performed NLSE simulations in the presence of weak dispersion (Bogoliubov regime), which reveal that, locally, the spectrum can exhibit some remarkable effects of spectral star formation and ray condensation in k-space, although their properties and underlying mechanisms depend on the nature of the random fluctuations as well as the dispersive properties of the waves (Supplementary Note 4 and Supplementary Figs 8-13). These preliminary simulations open an array of interesting questions, which will be considered in future works specifically devoted to this important issue of acoustic turbulence.
From a broader perspective, our photon fluid system can be viewed as a general platform for the experimental study of this acoustic turbulent regime. Furthermore, thanks to its widely tunable degree of nonlocality, it also opens the door to the experimental study of a variety of collective incoherent wave phenomena. We briefly comment here the applicability to extreme wave events, also called rogue, killer or freak waves 15,46,60 . Although these phenomena have been shown to emerge from a turbulent state, so far, the rogue wave itself has been always considered as being inherently a coherent localized entity. Here, by considering the shock and collapse singularities in the focusing regime go0 ð Þ, we may expect that the incoherent wave as a whole would lead to the formation of an extreme high-amplitude event, a feature that would be interpreted in analogy with an incoherent rogue wave phenomenon. Works are in progress in order to study the existence and the spontaneous emergence of this novel form of extreme events from a turbulent state of the field. The experimental study requires a highly nonlocal focusing nonlinearity, while usual pure liquids exhibit a defocusing thermal nonlinearity. Interestingly, however, a focusing thermal nonlinearity can be achieved by exploiting the Sôret effect in nanoparticle suspensions, which would thus couple the nonequilibrium behaviour of the turbulent field with the thermophoretic motion of the particles in the liquid 30,31 . Along these general lines, it would also be interesting to consider the relevance of the global collective phenomenological behaviours reported here within the general framework of opto-fluidics 61 . In this respect, while optical manipulation in strongly scattering densely packed (for example, biological) suspensions is usually considered impossible, a recent work demonstrated that some controlled manipulation of opaque suspensions can be achieved Simulations P = 0.05 W P = 2.5 W Zoom 1 for (f) Zoom 2 for (f)  (1) performed with the experimental parameters (l 0 c $ 250 mm, L ¼ 3 mm, s ¼ 59 mm), and corresponding zooms (g,h). As revealed by the zooms, a qualitative agreement is obtained between the experimental and numerical DSW periodic patterns. Note that simulations with s4150 mm do not reveal the formation of DSWs, which confirms the effective short-range regime of interaction in the experiment (Supplementary Note 1). NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9131 ARTICLE through the formation of a large-scale shock of nanoparticle concentration 62 . The long-range kinetic formalism developed here could be relevant to shed new light on the complex nonlinear mechanisms, which underlie this interesting example of highly nonlocal and diffusive opto-fluidic system. Furthermore, considering the universality of the NLSE, these collective nonequilibrium behaviours find applications in a multitude of disciplines in nonlinear science. Moreover, such large-scale incoherent singularities can also develop in the temporal domain, thanks to a noninstantaneous (instead of nonlocal) response function. At variance with spatial nonlocality, temporal nonlocality is characterized by a response function that is constrained by the causality condition 14,63 , which breaks the Hamiltonian structure, and thus introduces collective incoherent shocks and collapses phenomena of a different nature than those reported here. Then in addition to nonlocal wave systems discussed above through (refs 17-20,23-27), our work is also relevant to many physical systems involving a radiation-matter interaction featured by a finite-time nonlinear response.
In order to compare the simulations of the NLSE and LRVE models, we normalized the fields with respect to the size of the numerical window (L ¼ 1,200L in Figs 1 and 2). As discussed in the main text, a quantitative agreement has been obtained between NLSE and Vlasov simulations even in the strongly nonlinear regime. As far as we know, there is no rigorous proof of the validity of the LRVE in the strongly nonlinear regime. This fact may be qualitatively interpreted by using arguments similar to those of refs 55,64, in which it was shown that Gaussian statistics is preserved during the nonlinear evolution of the system.
Influence of the collision term on LRVE. The simulations reveal that in the highly nonlocal (that is, long range) regime of interaction, s=L ) 1, the dynamics of a random wave characterized by a non-homogeneous statistics is dominated by the LRVE, which indicates that the wave turbulence collision term is of higher order. This can be interpreted by considering the following collision term in the rhs of the LRVE, Coll½n k ¼ g 2 =b R R R Qðn k ; n k1 ; n k2 ; n k3 Þ T 2 k123 dðk 1 þ k 2 À k 3 À kÞdðk 2 1 þ k 2 2 À k 2 3 À k 2 Þdk 1 dk 2 dk 3 where Qðn k ; n k1 ; n k2 ; n k3 Þ ¼ n k1 n k2 n k3 n k ðn À 1 k þ n À 1 k3 À n À 1 k2 À n À 1 k1 Þ while the tensor T k123 accounts for the nonlocal interaction 1,6,14 , T k123 ¼ 1 4 ðŨ k1 þŨ k2 þŨ 31 þŨ 32 Þ, with U ij ¼Ũðk i À k j Þ, andŨðkÞ the Fourier transform of the response function U(r). A qualitative analysis of the LRVE reveals that the typical nonlinear propagation length scale is of the order l Vlas $ L 0 ðs=LÞ (ref. 65). On the other hand, the collision term slows down the dynamics by a factor of order equal to or larger than ðs=LÞ 2 : l nonloc Coll \l loc Coll ðs=LÞ 2 (ref. 14), where l loc Coll is the nonlinear propagation length scale of the local collision term (U(x)-d(x), and thus T k123 ¼ 1).
Considering a mode of the order kB1/L, we have typically l loc Coll $ L 0 , so that l nonloc Coll \l Vlas ðs=LÞ. Then in the long-range regime, s=L4 41, the collision term is of higher order with respect to LRVE (note that l loc Coll $ L 0 should be considered as a lower bound, since l loc Coll increases for higher modes k4 41/L required to verify the weakly nonlinear regime of interaction, o nl =o lin ¼ 1=ðk 2 l loc Coll Þo o1 (ref. 3)). In the highly nonlocal configuration of the experiment, we have s=LI2Â10 2 (Fig. 4), which indicates the irrelevance of the collision term to describe the experimental results in the long-range regime. Note that in the limit of a local interaction, the generalized collisional Vlasov equation recovers the form of the kinetic equation recently considered in ref. 8 to study the stability of Kolmogorov-Zakharov spectra of turbulence.
Regularization and closure. The finite 'time' singularities (5) described by the hydrodynamic model (3)(4) are regularized by the NLSE and LRVE equations. The spectrogramň k ðr; zÞ evolves in the 2D phase-space (r, k) and can thus become 'multi-valued' (Fig. 2 and Supplementary Fig. 7). However, the derivation of reduced equations describing the regularization of the singularities (5) is a difficult task fundamentally related to a long-standing mathematical problem, namely achieving a closure of the infinite hierarchy of equations that govern the evolutions of k-moments in transport kinetic equations. Note that the wave turbulence closure is justified in the weakly nonlinear regime 1,2,6,7 , while the closure considered here concerns the opposite strongly nonlinear regime. We address this problem of regularization in the Supplementary Note 3 through the analysis of higher-order truncations of the hierarchy (Supplementary Fig. 7). Our study reveals that all higher-order k-moments suddenly become of the same order of magnitude nearby the singularities (5), which prevents a closure of the infinite hierarchy and thus a reduced description of the dynamics beyond the incoherent shock point. Also note that in the context of self-gravitating systems, a weak diffusive effect has been introduced in a pure phenomenological way in order to regularize the wave breaking shock singularity described by the inviscid Burgers equation, thus leading to the so-called 'adhesion model' 66 . To our knowledge, no rigorous theory has been developed to justify such an heuristic approach 51 .
Experiments. For the incoherent shock experiment, Fig. 4, the sample is a 15-cmlong, 2-cm-diameter cylindrical tube, with anti-reflection (at 532 nm)-coated windows, and is filled with a solution of methanol and graphene nano-flakes. The graphene flakes have strongly subwavelength dimension (an average thickness of 7 nm, Graphene Supermarket) and thus do not induce significant scattering of light. They do, however, induce additional absorption and then release the absorbed energy to the surrounding liquid in the form of heat-less than 20% measured absorption over the whole sample length has been measured a ¼ 0:013cm À 1 ð Þ . For the dispersive shocklets experiment (Fig. 5), the concentration of graphene flakes has been increased (DnC À 4 Â 10 À 4 ) so as to increase absorption a ¼ 1:7cm À 1 ð Þ , and thus reduce the nonlocal range. Note that, due to such a strong absorption, the optical power transmitted through the sample is significantly reduced in the dispersive shocklets experiment, so that the sample length has been reduced to 1 cm. Although light absorption during propagation is no longer negligible, its impact on the dynamics is still perturbative with respect to nonlinear effects, L nl ð ' 0:02cmÞ ( L a ð ' 0:58cmÞ. According to the scaling s / 1= ffiffi ffi a p (ref. 24), the nonlocal range has been reduced by one order of magnitude in the coherent shocklets experiment as compared with the incoherent shock experiment. This important reduction of the nonlocal range has been confirmed by a direct comparison between the experiments and the NLSE simulations, in which the precise value of s has been adjusted to reproduce the experimental results: A minimum value of s4500mm is required in order to simulate the collective behaviour of the incoherent shock experiments, while accurate simulations of the dispersive shocklets experiments require a significant reduction of the nonlocal length, typically s480 mm (Supplementary Note 1 and Supplementary Figs 1-3). The measurement of the experimental spectrogram (x-k x ) in Fig. 4j-l has been performed by recording the far-field spectrum of the optical beam for several different x-positions and by keeping fixed the y-position and the pump power. For more details on the experimental results and corresponding numerical NLSE simulations, see Supplementary Note 1.