Nonlinear dispersion relation in integrable turbulence

We investigate numerically and experimentally the concept of nonlinear dispersion relation (NDR) in the context of partially coherent waves propagating in a one-dimensional water tank. The nonlinear random waves have a narrow-bandwidth Fourier spectrum and are described at leading order by the one-dimensional nonlinear Schrödinger equation. The problem is considered in the framework of integrable turbulence in which solitons play a key role. By using a limited number of wave gauges, we accurately measure the NDR of the slowly varying envelope of the deep-water waves. This enables the precise characterization of the frequency shift and the broadening of the NDR while also revealing the presence of solitons. Moreover, our analysis shows that the shape and the broadening of the NDR provides signatures of the deviation from integrable turbulence that is induced by high order effects in experiments. We also compare our experimental observations with numerical simulations of Dysthe and of Euler equations.

www.nature.com/scientificreports/ experiments described at leading order by 1DNLSE. In particular the (Fourier) spectral signature of solitons has not been reported in this context. In this manuscript, we first use numerical simulations of 1DNLSE to investigate the NDR and the power spectral density -PSD-(in the space-time (ω, k) Fourier plane) in IT. We show that the PSD and the NDR present clear signature of the growing influence of solitons when the nonlinearity strength is increased. We then report on an experiment in which partially coherent (random) deep water surface gravity waves propagate along a long one-dimensional flume. We demonstrate that removing the carrier wave enables the accurate measurement of the NDR of wavefields having a narrow Fourier spectrum by using a very limited number of gauges. In experiments, the NDR is found to reveal both the existence of solitons embedded in the random field and a deviation from IT induced by higher-order nonlinear effects not taken into account in the 1DNLSE model. Simulations of Euler and Dysthe equations confirm the influence of high order terms and provide a deeper insight into the mechanisms behind the formation of coherent structures observed in the experiment.

Results
Numerical simulations (nonlinear Schrödinger equation). Considering unidirectional deep water gravity waves having a narrow spectrum, the surface elevation η(τ , z) is: where ψ r is the slowly-varying complex envelope, f 0 = ω 0 /(2π) and k 0 = ω 2 0 /g are the frequency and the modulus of the wavevector of the carrier wave respectively, z is the propagation distance and τ is the time measured in the laboratory frame. For deep water gravity waves, the dynamics of ψ(t, z) = ψ r (τ , z) is described at leading order by the focusing 1DNLSE 29 : where t = τ − z/c g , c g = g/(2ω 0 ) is the group velocity evaluated at the frequency f 0 , γ = k 3 0 and g is the gravity acceleration. In simulations and in experiments, the initial conditions are partially coherent waves produced from the linear superposition of numerous independent Fourier components having a Gaussian spectrum (see "Methods" and 19,21,30 ). It is useful to introduce the degree of nonlinearity of the wave propagation, which is given by the parameter Ŵ: where f ≪ f 0 is the initial spectral bandwidth evaluated at z = 0 and where ... denotes the averaging over time and/or realizations. Note that �|ψ| 2 � = 2�η 2 � and that Ŵ = BFI 2 where BFI is the Benjamin-Feir Index 31,32 (see "Methods"). It has been shown that large values of Ŵ lead to the formation of rogue waves characterized by heavy-tailed statistical distributions of the surface elevation 16,19,30,31,[33][34][35] .
We have first performed numerical simulations of Eq. (2) for three values of Ŵ (see "Methods" for details). Fig. 1a, represents the typical spatio-temporal dynamics of IT developing from partially coherent waves in the focusing regime of 1DNLSE. Remarkably, the number of isolated pulses emerging in the random wave field increases with the nonlinearity. Note that our numerical simulations reveal the presence of elastic collisions (see for the example the white circle in Fig. 1), a signature of solitons in integrable systems 13 .
We define the space-time double Fourier transform as: The power spectrum (PSD) | ψ(ω, k)| 2 of one realization of ψ(z, t) is plotted in blue in Fig. 1b. As pointed out in 25 and reported below in "Methods", the straight lines in the k − ω space are signature of the solitons observed in the spatio-temporal dynamics (Fig. 1a). Note that the theoretical description of the statistical distribution of the slopes of the straight lines-associated to the velocities of the solitons-is a theoretical open question. We also plot the spectrum �| ψ(ω, k)| 2 � averaged over several realizations in Fig. 1c. The linear dispersion of 1DNLSE reads k(ω) = ω 2 /g 36 and is plotted in dashed black lines in Fig. 1b,c. The NDR can be defined as the peak wavevectors k of the PSD for each value of ω in the (ω, k) plane 25 . For moderate nonlinearities and narrow Fourier spectra, the nonlinear NDR gets shifted from the linear one and reads 25,36 : As expected from Eq. (5), the NDR computed from numerical simulations shifts toward negative values of k (see red dashed lines in Fig. 1b,c). The Fig. 1c,d show that, the PSD broadens around the NDR when Ŵ increases because of the energy exchange among Fourier modes. This broadening phenomenon is well-known in standard WT with resonant interactions 1-3,37 but, to the best of our knowledge, it has not been reported in IT where resonances are forbidden (see "Methods"). The k−spectrum at ω = 0 , i.e. n k = �| ψ(0, k)| 2 � is plotted in Fig. 1d. For small values of Ŵ , the maximum of this curve coincides with the value predicted by the weakly nonlinear theory (i.e. Eq. (5), red dashed line in Fig 1d). Note that, in this weakly nonlinear regime, n k can be empirically fitted by a Lorentzian distribution (purple dashed line in Fig. 1d). To the best of our knowledge, this remarkable fact, not known yet in IT, has www.nature.com/scientificreports/ not been described theoretically. At higher nonlinearities, the number of straight lines associated with solitons (individually observed in Fig. 1a,b for Ŵ = 0.65 ) increases. The NDR of a single soliton follows the straight line where c s is the group velocity of the soliton in the (t, z) plane, ω s is the central frequency of the soliton and k s = k(ω s ) = ω 2 s /g follows the linear dispersion (see 25 and "Methods"). As a consequence, the nonlinear phase shift acquired by solitons during their propagation places the solitonic lines well below the dispersion parabola described by Eq. (5) (see Fig. 1b).
Consequently, the strong and asymmetric broadening of the PSD profile n k toward negative values of k (Fig. 1d) can be considered as a signature of the increasing number of solitons in the strongly nonlinear regime of IT.

Experiments.
In order to investigate experimentally the NDR described above, we have used the setup described in 38 and schematically shown in Fig. 2a. Regimes close to IT can be achieved in unidirectional deep water waves having a narrow Fourier spectrum. Such waves are generated at one end of a 148 m long, 5 m wide and 3 m deep wave flume by a computer-assisted flap-type wavemaker (see Fig. 2a). The flume is equipped with an absorbing device strongly reducing wave reflection at the opposite end 38 . The surface elevation η(τ , z) is measured by using 20 equally spaced resistive wave gauges that are installed along the water tank at distances z j = 6 j m, j = 1, 2, ...20 from the wavemaker located at z = 0 m. This provides an effective measuring range of 120 m and a resolution of the k−spectrum of 2π/120 rad m −1 (see "Methods"). The envelope ψ(t, z = 0) of the www.nature.com/scientificreports/ surface elevation has a central frequency f 0 = 1.15 Hz and it is designed using the same procedure as the one used in our numerical simulations, i.e. by performing the linear superposition of a large number of independent random Fourier modes. The degree of nonlinearity is varied by changing either the averaged amplitudes or the initial spectral width f of the waves generated in the water tank.
A typical temporal evolution of the surface elevation experimentally recorded at the first gauge ( z = 6 m) is plotted in the Fig. 2b. The slowly varying amplitude ψ(z, t) is determined by using Hilbert Transform following techniques described, e.g., in Ref. 29 . Typical spatio-temporal evolution of |ψ| is plotted in Fig. 3a where we use the retarded time t = τ − z/c g . When Ŵ increases, we observe the emergence of pulses localized in space and time. The emerging pulses become narrower and most of them achieve a negative speed in the (t, z) diagram (see below).
As expected from our numerical simulations, the experimental PSD broadens and shifts toward the negative values of k (Fig. 3b). However, at high nonlinearity, the measured NDR deviates significantly from the numerical simulations and becomes asymmetric with ω . This phenomenon is the well-known "frequency downshift" of surface gravity waves induced by high order nonlinearities responsible for the negatives speeds observed in the (t, z) diagram shown in Fig. 3a 39 . The numerical simulations of Dysthe equation (see below) and of the Euler equations (see Supplementary Material) confirm that this shape of the NDR is indeed induced by effects (not included in 1DNLSE) which break integrability of the wave equation.
Simulations of the Dysthe equation. The 1DNLSE (Eq. 2) is established under the assumptions of weak nonlinearity of the wave field as well as the narrowbandedness of its energy content. As it could be expected, at high values of the nonlinear parameter Ŵ , the measured NDR deviates from the one computed by performing numerical simulations of the 1DNLSE. The Dysthe equation (a higher-order nonlinear and not integrable generalized version of the 1DNLSE) provides a simple model that reproduces qualitatively the ω asymmetry of the NDR observed in experiments. Dysthe equation can be expressed as follows 40 : where H stands for the Hilbert transform defined as follows:F (H(f (t))) = −i sign(ω)F (f (t)) , where F represents the Fourier transform and sign is the signum function. The extra terms labelled 'a' , 'b' and ' c' in Eq. (6) represent higher order terms in a perturbative approach of Euler equation where the small parameter is the spectral width f . If f → 0 , the derivatives terms 'a' , 'b' and ' c' vanish whereas the dominant cubic term ik 3 0 |ψ| 2 ψ remains unchanged. In order to investigate the influence of additional terms present in Dysthe equation, we have simulated the nonlinear propagation of identical initial conditions used in Fig. 1a, Ŵ = 0.65, including Dysthe terms labeled 'a' , 'b' , and ' c' separately as shown in Fig. 4. Parameters used in the numerical simulations are the same as those used for the numerical integration of the 1DNLE reported above.
As one can see in Fig. 4, the derivatives included in the terms 'a' and 'b' contribute to the nonlinear spectral shift, leading to the change of solitons' group velocity (asymmetry in ω ). The term ' c' reduces the shift of solitons  www.nature.com/scientificreports/ leading to a picture similar to 1DNLSE but with a smaller amplitudes of the localized structures which can be seen in the corresponding NDR plots. Importantly, we have also run numerical simulations Euler equations by using high-order spectral (HOS) method (see Supplementary Material). The results of these simulations are close to those obtained with the Dysthe equation.

Nonlinear shift and broadening of the nonlinear dispersion relation.
Our data analysis reveals that higher-order effects discussed above reduce the number of solitons embedded in random waves: indeed, we found that for Ŵ ≤ 1 , the straight lines in the (ω, k) space-signatures of solitons-appear much less frequently than in 1DNLSE simulations. Nevertheless, we have also observed these spectral signatures of solitons at extremely high nonlinearity (which is achieved with small values of the f , see the fourth column of Fig. 3).
The shapes of the PSD and of the NDR provide another signature of the deviation from IT: while the NDR predicted from the 1DNLSE becomes asymmetric at high nonlinearity, the experimental PSDs profiles n k = �| ψ(0, k)| 2 � in Fig. 3c coincide with a Lorentzian fit for all values of Ŵ , a result also observed for HOS simulations.
The influence of the nonlinearity and of high order effects can be quantified with the help of the position k = k M (Ŵ) of the maximum and of the full width at half maximum �k(Ŵ) of n k (see Fig. 5a,b respectively). Figure 5a shows that, both in experiments and simulations, k M (Ŵ) evolves more slowly than predicted by the weakly nonlinear theory-see Eq. (5). Moreover, because of the smaller number of solitons, the shift of the NDR toward negative values of k is weaker in experiments and HOS simulations than in IT (1DNLSE simulations).

Discussion
Our work provides new insight into the spectral properties of unidirectional nonlinear random waves. Our numerical simulations of the 1DNLSE show how the number of emerging solitons is directly related to the strength of nonlinearity in IT. Moreover, while a previous study has focused on the position of the maxima of the NDR 25 , we also investigate the broadening of the NDR. Simulations reveal that the NDR provides a spectral signature of the solitons embedded in the random waves (asymmetric broadening in the k direction at high nonlinearity). Further investigations are needed to establish the theoretical link between the nonlinear spectra computed in the framework of IST and the broadening of the NDR (Fourier spectra).
From the experimental point of view, by focusing our analysis on the slowly-varying amplitude ψ , we were able to measure very accurately slight deviations from the predicted nonlinear dispersion relation. The results reported here provide new insights into an old fundamental problem of hydrodynamics: the measurement of the dispersion relation of random surface gravity waves (see for example 5   www.nature.com/scientificreports/ experimental works have been devoted to this question, see e.g. 5,36,[41][42][43][44][45][46][47] . For water tanks equipped with transparent side walls, cameras may be used to record the full spatio-temporal dynamics 23 . Otherwise, a large number of gauges is needed 5 . The central point of our strategy is to remove the carrier wave frequency in order to retrieve the NDR of the slowly varying amplitude of waves having a narrow spectrum. Our simple experimental technique can be easily implemented in further investigations of the NDR of unidirectional water waves by using a very limited number of probes (20 probes here instead of 384 probes in 5 for example, see "Methods").
Our work contributes to the understanding of the NDR in nonlinear waves systems that involve solitons. Recently, a non trivial NDR has been established for finite gap solutions in the context of NLS soliton gases 48 . The possible relationship between this NDR derived for soliton gases and the Fourier NDR studied here is an open fundamental question. On another side, the concept of NDR is receiving interest in the photonics community, where it has been for example recently used to explain the effect of dissipative soliton hopping in a photonic dimer 49 or to describe noise properties of a soliton frequency comb in a synchronously pumped cavity 50 . In this article, we have demonstrated that the NDR measured in deep water waves experiments is close to the one predicted by using the 1DNLSE for low nonlinearity. For high nonlinearity, the measured NDR exhibits a signature of solitons but experiments deviates from IT because perturbative effects not taken into account in the 1DNLSE limit the emergence of solitons. It has been demonstrated that optical fiber experiments can be very close to integrability for high nonlinearity 19,30 . The measurement of the NDR in optical fibers is extremely challenging but it has been recently demonstrated in a double loop fiber devices 51 . We hope that our work will also stimulate further investigation of the NDR of IT in photonics.

Nonlinear Schrödinger equation (1DNLSE). Strength the nonlinearity: Benjamin-Feir index. In order
to compare nonlinearity and group velocity dispersion in the framework of the 1DNLSE (Eq. 2), it is useful to introduce a linear and a nonlinear propagation length as follows: where f is a typical initial spectral bandwidth and P 0 = �|ψ(t, z = 0)| 2 � where ... is the averaging over time and/or realizations. The degree of nonlinearity of the wave propagation is given by the parameter Ŵ = z lin /z nlin (see Eq. 3).
Note that in the context of ocean waves, Ŵ = BFI 2 where BFI is the Benjamin-Feir Index 31,32 . BFI index can also be expressed as follows: where ε = k 0 √ 2σ is the wave steepness with σ 2 = �η 2 � = �|ψ| 2 �/2 and f and f 0 are the average spectral width and the central frequency of the initial wave packets respectively. Along the propagation in the focusing regime of 1DNLSE, the statistics deviates from Gaussianity and the probability density function of the wave amplitude becomes heavy-tailed. Note that the statistical characteristics of partially coherent waves are very different than plane waves initially perturbed by noise which is also investigated in 25 . For a comparison between the two cases, refer e.g. to 52 .
Numerical simulations of Eq. (2) are realized using step-adaptive high order Runge-Kutta method. We construct different initial conditions using the random phase approach where a uniformly distributed phase is added to every Fourier component of a Gaussian spectrum with f = 0.2 Hz. Typical temporal windows used in the numerical simulations correspond to 100 seconds. Three parameters of simulations depend on the value of the steepness: the number of points N, the length of propagation L max and the number of realizations N sample . We separate the numerical studies into three ranges of steepness ε (see Table 1).
(7) z lin = g (2π�f ) 2 and z nlin = 1 γ P 0 , www.nature.com/scientificreports/ In order to reconstruct numerically NDR, we multiply the spatiotemporal diagram by a Super-Gaussian window with power 15 along z direction avoiding thereby undesirable effects related to the Fourier analysis of non-periodic signals.
Nonlinear dispersion of an isolated soliton. The fundamental soliton solution of the 1DNLSE (Eq. 2) reads 56 : where the duration of the soliton is τ = 2/(γ g|ψ 0 | 2 ) . k s = k(ω s ) obeys the linear dispersion relation, c s = c(ω s ) = dk/dω = 2ω s /g is the group velocity of the soliton in the (z, t) plane and A = γ |ψ 0 | 2 /2 . The double Fourier transform of ψ S (t, z) is given by: As a consequence, the NDR of a single soliton follows the straight line k = c s ω + (k s − c s ω s − 1 2 γ |ψ 0 | 2 ) 25 .

Measurement of the nonlinear dispersion relation
Resolution of the measurement of k. The key point in our approach is to remove the carrier wave before computing the NDR. This allows us to reveal the details of the NDR of the slowly varying envelop. In the water tank, the gauges are separated by 6 m and the maximum measurable wavevector is around 1.05 m −1 . As a consequence, contrary to Taklo et al. who use 384 probes, we do not resolve the wavevectors of the carrier wave 2π/ 0 ≃ 5.3 m −1 and of the harmonics. Our strategy enables the accurate measurement of the NLDR of the slow varying envelop of the wave by using only 20 probes. This provides an effective measuring range of 120 m with a resolution of the measurement of the k−spectrum of �k min = 2π/120 rad m −1 . Note finally that, in 5 , the accuracy of the measurement of k given by the length of the water tank is �k min /k 0 = 0.014 while our setup enables an accuracy of �k min /k 0 < 0.01. In order to measure the averaged spectra and NDR, we use 3, 6 and, 3 experimental runs with a duration of 512 s for Ŵ = 0.12 , 0.33, 0.65, respectively. One run of 128 s have been used for Ŵ = 6.18.
Note that in NLS and HOS (high-order spectral, see Supplementary Material) simulations, the chosen lengths of propagation depend on the parameters and vary typically from 300 to 500 m. The uncertainty of measurement of k M and k is therefore significantly lower in simulations than in experiments.
Evaluation of the full width at half maximum k of the position of the maximum k M . In Fig. 4, we report the evaluation of the full width at half maximum k and of the position of the maximum k M of the function f (k) = | ψ(k, ω = 0)| 2 in 1DNLSE, HOS simulations and in experiments. The accuracy of the measurement of k and k M is limited both by the discretization of k (see above the uncertainty k min ) and by the random fluctuations of f(k). In order to overcome these difficulties, when it is appropriate, we evaluate k and k M by using best fitting procedure with Lorentzian function.

Data availability
The datasets generated and/or analysed during the current study are available in the Figure data: Nonlinear dispersion relation in integrable turbulence repository, https:// zenodo. org/ record/ 65954 29.