# Irreversible thermalization vs reversible dynamics mediated by anomalous correlators: Wave turbulence theory and experiments in optical fibers

Theo Torres, J. Garnier, L. Zanaglia, M. Ferraro, C. Michel, V. Doya, J. Fatome, B. Kibler, S. Wabnitz, Antonio Picozzi, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Theo Torres, J. Garnier, L. Zanaglia, M. Ferraro, C. Michel, et al., Irreversible thermalization vs reversible dynamics mediated by anomalous correlators: Wave turbulence theory and experiments in optical fibers. 2025. [hal-05434515](https://hal.science/hal-05434515v1)

HAL Id: hal-05434515

https://hal.science/hal-05434515v1

Preprint submitted on 29 Dec 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Irreversible thermalization vs reversible dynamics mediated by anomalous correlators:
Wave turbulence theory and experiments in optical fibers

T. Torres¹, J. Garnier², L. Zanaglia³, M. Ferraro¹,⁴, C. Michel³, V. Doya³,

J. Fatome¹, B. Kibler¹, S. Wabnitz⁵, A. Picozzi¹, G. Millot¹

¹ Université Bourgogne Europe, CNRS, Laboratoire Interdisciplinaire Carnot de Bourgogne ICB UMR 6303, 21000 Dijon, France

² CMAP, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

³ Université Côte d'Azur, CNRS, Institut de Physique de Nice, Nice, France

⁴ Department of Physics, University of Calabria, Via P. Bucci, Rende, 87036, (CS), Italy and

⁵ Department of Information Engineering, Electronics, and Telecommunications, Sapienza University of Rome, Via Eudossiana 18, Rome, 00184, Italy

We theoretically and experimentally investigate spontaneous self-organization in a conservative (Hamiltonian) turbulent wave system, operating far from thermodynamic equilibrium. Our system is governed by two coherently coupled nonlinear Schrödinger equations, describing the polarization evolution of light in a dispersive nonlinear optical fiber. The analysis reveals the emergence of two fundamentally distinct turbulent regimes. In a first regime, the waves undergo a slow, irreversible thermalization process, which is accurately described by the wave turbulence kinetic equation and the associated *H*-theorem of entropy growth. In stark contrast with this expected irreversible process, we identify a second different regime, where strong phase-correlations spontaneously emerge, giving rise to a fast reversible oscillatory dynamics of the normal correlator and anomalous phase-correlator. Experimental observations confirm the occurrence of both irreversible thermalization and reversible dynamics mediated by the anomalous correlated fluctuations.

*Introduction.* Non-integrable Hamiltonian systems of random waves generally undergo thermalization – an irreversible evolution toward thermodynamic equilibrium, characterized by a state of maximum entropy. In the weakly nonlinear regime, this process is accurately described by the wave turbulence (WT) theory [1–6], which has been successfully applied to various physical systems [7–12], including optical waves [13–18]. The WT kinetic equation (KE) describes the actual irreversible evolution to the Rayleigh-Jeans equilibrium distribution, which is expressed by a *H*-theorem of entropy growth. Interestingly, light thermalization has recently been studied theoretically [19–21] and experimentally [21–24] in multimode fibers, thus spurring the emerging field of optical thermodynamics [19, 25–35].

Although rigorous derivations of the WT KE have been accomplished for certain partial differential equations (PDEs) [36–38], the broader validity of the WT theory across a wide range of nonlinear wave-interacting PDEs remains unclear. A central requirement in deriving the KE is the demonstration that *phase-correlations* among the random waves can be neglected.

Phase-correlations are known to underlie large-scale coherent structures, such as solitons or condensates, which emerge from strongly nonlinear turbulent regimes [2–8, 13, 39–43]. These coherent structures are, by their nature, phase-correlated entities. In contrast, we will focus here on phase-correlations in the *weakly nonlinear regime* of random dispersive waves.

Phase-correlations – also referred to as ‘anomalous correlations’ – were central to the original Bardeen–Cooper–Schrieffer theory of superconductivity [44]. They have since been extensively studied within the framework of *S*-theory [45–47], where anomalous

correlations arise from an external forcing (pumping) that parametrically excites a dissipative magnetic system [47]. From a different perspective, recent theoretical advances revealed the presence of phase-correlations in conservative (unforced) wave systems. Anomalous phase-correlations have been identified in Fermi–Pasta–Ulam–Tsingou chains [48], or in the coupled nonlinear Schrödinger equation (NLSE), arising from strong convection between wave components [49]. In *scalar* wave systems, the mechanism underlying the emergence of phase-correlations between different frequencies has been recently elucidated, by showing that they may only arise from Hamiltonian terms that break phase invariance, whereas systems that preserve phase invariance preclude such anomalous correlations [50].

In this Letter, we investigate a *vector* system of coherently coupled NLSEs, which describe the polarization evolution of temporally incoherent light waves in optical fibers. We identify two fundamentally different regimes: an irreversible thermalization, and a reversible turbulent dynamics driven by phase-correlations. Unlike recent 2D spatial studies in multimode optical fibers [21–24], we examine here the 1D temporal dynamics of incoherent waves propagating in a single-mode fiber [14–18]. By using the standard WT theory (neglecting phase-correlations), we describe nonequilibrium thermalization accurately. In contrast to this expected irreversible process, we show that phase-correlations can grow exponentially from fully uncorrelated initial random waves, leading to a reversible oscillatory turbulent dynamics. These anomalous phase-correlations emerge among different wave components in a phase-invariant vector system, hence they are of different nature from those investigated in Ref. [50]. Experimental results provide evidence
of both irreversible thermalization and reversible phase-
correlated dynamics. From a broader perspective, this
work advances our understanding of far-from-equilibrium
self-organization processes in closed (Hamiltonian) tur-
bulent wave systems.

*Model.* We consider a generic model for vector phenomena based on coherently coupled NLSEs. It describes spinor Bose–Einstein condensates with field-induced spin-flip coupling [51], allowing studies of quantum phase transitions [52], or bubble nucleation in ferromagnetic superfluids [53]. The coherently coupled NLSEs also describe the polarization dynamics of an optical wave propagating in a weakly birefringent fiber [9, 54]:

$$
i\partial_z \mathbf{u} = -\beta \partial_{tt} \mathbf{u} + \alpha \boldsymbol{\sigma} \mathbf{u} + \gamma (\kappa \mathbf{u}^\dagger \mathbf{u} \mathbf{u} + \rho \mathbf{u}^T \mathbf{u} \mathbf{u}^*), \quad (1)
$$

where **u**(t, z) = (u<sub>x</sub>, u<sub>y</sub>)<sup>T</sup> is the vector field in the linear
polarization basis, and the superscripts (T, †) stand for
the transpose, and conjugate transpose, operations. As
usual in optics, the distance z of propagation plays the
role of an evolution ‘time’ variable, while t denotes the
retarded time in a reference frame moving with the waves.
γ is the non-linear parameter, β is the group-velocity
dispersion, while κ and ρ are dimensionless interaction
coefficients; we will consider the case κ = 2ρ relevant
to our experiments. The coherent coupling parameter α
originates from the weak birefringence of the optical fiber,
with the matrix σ = diag(+1, −1). We assume α > 0
without loss of generality. The vector NLSE (1) conserves
the power (particle number) N = ∑<sub>μ</sub> N<sub>μ</sub>, where N<sub>μ</sub>(z) =
1/T ∫<sub>0</sub><sup>T</sup> |u<sub>μ</sub>|<sup>2</sup>dt, with T the size of the numerical window,
and μ = x, y. It also conserves the Hamiltonian H = E+
U, with a linear E, and a nonlinear U, contribution [56].

Let us introduce the usual nor-
mal correlators
$$
n_{\mu}(\omega, t, z) = \int \langle u_{\mu}(t + \tau/2, z) u_{\mu}^{*}(t - \tau/2, z) \rangle \exp(-i\omega\tau) d\tau \quad (\mu = x, y).
$$
The dependence of $n_{\mu}$ on the ('spatial') variable $t$
accounts for possible statistical inhomogeneities in
the random waves, while the angle brackets denote
an average over the random initial conditions. Fur-
thermore, the theorem in Ref.[50] for phase-invariant
scalar systems can be extended to vector systems such
as the vector NLSE (1), which implies the absence
of anomalous correlations among distinct frequencies.
However, as we show below, phase-correlations may
emerge between different wave (i.e., polarization)
components at the same frequency: these will be
characterized by the anomalous correlator $m(\omega, t, z) = \int \langle u_y(t + \tau/2, z) u_x^*(t - \tau/2, z) \rangle \exp(-i\omega\tau) d\tau$.

Standard WT theory: Irreversible thermalization. We consider the weakly nonlinear regime, where linear effects dominate over nonlinear effects |E/U| ≫ 1. Using the standard WT theory [1–6, 13, 18], we derive the KE under the usual assumptions: (i) statistical homogeneity, so that n_μ are t-independent, and (ii) absence of phase-correlations between u_x and u_y, i.e., the anomalous correlator is zero at any ‘time’ z, m(ω, t, z) = 0. The KE governing the evolution of the averaged spectra

<figure><img src="image_3.png" /><figcaption>FIG. 1: Irreversible thermalization. Evolution during the propagation of N<sub>x,y</sub>(z) (a), and corresponding degree of polarization P(z) [Eq.(3)] (b), obtained from the simulation of the NLSE (1) (5 realizations, dark and light blue lines), and the simulation of the WT KE (2) (red lines). Parameters: αL<sub>nl</sub> = 1, ΔN<sup>0</sup> = 0, στ<sub>0</sub> = 0.8π, κ = 2/3, with L<sub>nl</sub> = 1/(γN) the nonlinear length and τ<sub>0</sub> = √|β|L<sub>nl</sub> the ‘healing time’ (|E/U| ≃ (στ<sub>0</sub>)<sup>2</sup>).</figcaption></figure>

$\mathbf{n}(\omega, z) = (n_x, n_y)^T$ then takes the usual form

$$
\partial_z \mathbf{n}(\omega, z) = \kappa^2 \mathrm{Coll}_i[\mathbf{n}] + \rho^2 \mathrm{Coll}_c[\mathbf{n}]. \quad (2)
$$

The collision terms are rather cumbersome, see [56]. They are cubic nonlinear terms Coll<sub>i,c</sub>[n] ~ n<sup>3</sup>, describing the four-wave interaction as a collisional gas of ‘particles’. The KE (2) conserves N = &Sigma;<sub>&mu;</sub> N<sub>&mu;</sub>, with N<sub>&mu;</sub>(z) = 1/2&pi; &int; n<sub>&mu;</sub>(&omega;, z)d&omega; (&mu; = x, y), and the linear energy E [56]. At variance with the NLSE (1), the KE (2) is irreversible, as expressed by a H-theorem of entropy growth, dS/dz &ge; 0, where S = &Sigma;<sub>&mu;</sub> S<sub>&mu;</sub>, and S<sub>&mu;</sub>(z) = 1/2&pi; &int; log (n<sub>&mu;</sub>(&omega;, z))d&omega; is the nonequilibrium entropy of the &mu;-th component (&mu; = x, y).

We have performed numerical simulations of NLSE (1). As initial condition, we consider two uncorrelated random waves $u_{x,y}(t, z = 0)$ with zero mean, Gaussian statistics and Gaussian-shaped initial spectrum $n_x^0(\omega) = n_y^0(\omega)$, i.e., $N_x^0 = N_y^0$, with $n_\mu^0(\omega) = \sqrt{2\pi N_\mu^0} \exp(-\omega^2/(2\sigma^2))/\sigma$. Let us introduce the degree of polarization of the optical wave [57]:

$$
\mathcal{P}(z) = \sqrt{\Delta N(z)^2 + 4|M(z)|^2}/N, \quad (3)
$$

where $M(z) = \frac{1}{2\pi} \int m(\omega, z) d\omega = \langle u_y(t, z) u_x^*(t, z) \rangle$ denotes the integrated anomalous correlator, and $\Delta N(z) = N_y(z) - N_x(z)$ the unbalanced power distribution. Note that $\mathcal{P}$ is bounded by 0 and 1, which correspond to unpolarized, and fully polarized, waves, respectively.

We have also performed numerical simulations of the
WT KE (2): as shown in Fig. 1, a good agreement (with-
out adjustable parameters) is obtained with direct sim-
ulations of the NLSE (1). Recalling that the KE ne-
glects phase-correlations, such an agreement means that
anomalous correlations do not grow, i.e., M(z) ≃ 0, and
P(z) ≃ ΔN(z)/N. Accordingly, the thermalization pro-
cess is characterized by an irreversible transfer of power
from N<sub>x</sub> to N<sub>y</sub> due to their opposite contributions to the
linear energy, thus leading to the emergence of a non-
vanishing degree of polarization [56].
Reversible dynamics mediated by phase-correlations. We now show that the thermalization process can be inhibited by the spontaneous emergence of phase-correlations. To clarify this regime in a general framework, we do not implicitly assume the waves to obey homogeneous statistics,

i.e., we leave $n_μ(ω, t, z)$ and $m(ω, t, z)$ to depend on a slow non-homogeneous variation in the $t$ variable. Starting from the NLSE (1), we derive in the Supplementary Material the equations governing the coupled evolution of the normal and anomalous correlators [56]:

$$ \partial_z n_x(\omega, t) = -2\beta\omega\partial_t n_x + \gamma\partial_t(2N_x(t) + \kappa N_y(t))\partial_\omega n_x + 2\gamma\kappa\partial_t M^r(t)\partial_\omega m^r + 4\gamma\kappa M^r(t)m^i, \quad (4) $$

$$ \begin{aligned} \partial_z m(\omega, t) = & -2\beta\omega\partial_t m + \gamma(1 + \kappa/2)\partial_t(N_x(t) + N_y(t))\partial_\omega m + \gamma\kappa\partial_t M^r(t)\partial_\omega(n_x + n_y) \\ & - im(\gamma(2 - \kappa)(N_y(t) - N_x(t)) - 2\alpha) - 2i\gamma\kappa(n_x - n_y)M^r(t), \end{aligned} \quad (5) $$

with $M(t, z) = \frac{1}{2\pi} \int m(\omega, t, z)d\omega$, $N_μ(t, z) = \frac{1}{2\pi} \int n_μ(\omega, t, z)d\omega$, while the superscripts $m^{r,i}$ ($M^{r,i}$) denote the real and imaginary parts of $m$ ($M$). The equation for $n_y$ follows from Eq.(4) by exchanging $x \leftrightarrow y$ and reversing the sign of the last term.

The first two terms on the right-hand side of Eq.(4) correspond to the well-established Vlasov coupling between two uncorrelated random waves. Here, the generalized Vlasov-like Eqs.(4-5) provide an original extension accounting for phase-correlations, $m(\omega, t)$. In contrast to the KE (2), the system (4-5) is formally reversible in 'time' $z$ [56]. To avoid confusion with the usual KE (2), we shall refer to Eqs.(4-5) (and Eq.(8) below) as the anomalous-correlator kinetic equation (AC-KE).

We consider an initial condition of two uncorrelated random waves, of zero mean, with homogeneous statistics and spectra $n_μ(\omega, z = 0) = n_μ^0(\omega)$, so that the initial anomalous correlator is zero, $m(\omega, t, z = 0) = 0$. We carry out a linear stability analysis of Eqs.(4-5) around this state. Using the Laplace-Fourier transform $\hat{m}(\omega, \Omega, \lambda) = \int_0^z dz \int m(\omega, t, z) \exp(-\lambda z - i\Omega t)dt$, we obtain the dispersion relation $\lambda(\Omega)$ of the anomalous correlator [56]:

$$ \frac{2\pi}{\gamma\kappa} = \sum_s \int \frac{n_x^0(\omega) - n_y^0(\omega) - s\Omega\partial_\omega(n_x^0(\omega) + n_y^0(\omega))/2}{is\lambda - 2s\beta\omega\Omega - \gamma(2 - \kappa)\Delta N^0 + 2\alpha} d\omega, $$

where $\Delta N^0 = N_y^0 - N_x^0$, and $s = \pm 1$. Assuming the initial spectra Gaussian-shaped, we compute the corresponding growth-rate Re[$\lambda(\Omega)$]. The analysis reveals that, in general, the homogeneous mode with $\Omega = 0$ is the most unstable, with the maximum growth rate of the anomalous correlator:

$$ \lambda(\Omega = 0) = 2\sqrt{\alpha(2\gamma\Delta N^0/3 - \alpha)}, \quad (6) $$

where we have considered the case $\kappa = 2/3$ relevant to our experiments. The instability criterion then reads

$$ \alpha L_{nl} < (2/3)\Delta N^0/N, \quad (7) $$

revealing that the emergence of phase-correlations requires an unbalanced initial population, $N_y^0 > N_x^0$. The

<figure><img src="image_4.png" /><figcaption>FIG. 2: Reversible turbulent dynamics. (a) Exponential growth of the anomalous correlator starting from initially uncorrelated random waves $|M(z = 0)| \simeq 0$: Simulations of the NLSE (1) (100 realizations, black lines), and corresponding average (green line). The red line reports the theoretical prediction Eq.(6). (b) Nonlinear dynamics on the Poincaré sphere: Simulation of the NLSE (1) (red line), and corresponding theoretical prediction from the simulation of AC-KE (8) (black line), showing the periodic exchange among the normal and anomalous correlators. The blue line reports the homoclinic orbit emanating from the unstable fixed point $S_1 = +S_0$, i.e., $|M| = 0$. Parameters: $\alpha L_{nl} = 0.2$, $\sigma\tau_0 = 4\pi$, $\Delta N^0/N = 0.6$, $\kappa = 2/3$, $L_{nl} = 1/(\gamma N)$, $\tau_0 = \sqrt{|\beta|L_{nl}}$ ($|E/U| \simeq (\sigma\tau_0)^2$).</figcaption></figure>

simulations of the NLSE (1) confirm the theoretical prediction (6), see Fig. 2(a).

The fact that the homogeneous mode $\Omega = 0$ is the most unstable, indicates that the subsequent nonlinear dynamics described by Eqs.(4-5) preserves the homogeneous statistics, i.e., $n_μ(\omega, t)$ and $m(\omega, t)$ are $t$-independent. In this limit, the dynamics for the normal and anomalous correlators (4-5) can be recast in the form of a Stokes vector $S(z) = (\Delta N, 2M^r, 2M^i)^T$ that rotates on the surface of the Poincaré sphere [56]:

$$ \partial_z S(\tau, z) = R(z) \times S(z), \quad (8) $$

with $R(z) = (2\alpha - (2 - \kappa)S_1(z), -2\kappa S_2(z), 0)^T$. This formalism differs substantially from that typically used to describe polarization effects for fully coherent stationary waves [9], as originally introduced in [58]. Here, in
particular, the radius of the Poincaré sphere depends on the initial condition through the degree of polarization $S_0 = (\sum_{j=1}^3 S_j^2)^{1/2} = PN$, which is thus conserved during evolution, $P(z) = \text{const.}$ Consequently, the nonlinear dynamics described by the AC-KE (8) is essentially periodic, featuring oscillations that involve a reversible exchange between the normal correlator $\Delta N(z)$ and the anomalous correlator $M(z)$, while conserving

$$ \mathcal{P}^2 N^2 = \Delta N^2(z) + 4|M(z)|^2 = \text{const.} \quad (9) $$

This reversible dynamics is in agreement with the simulations of the NLSE (1), see Fig. 2(b). Note that the fast, reversible dynamics, governed by the quadratic nonlinearities in the AC-KE (4-5) and (8) in Fig. 2, stands in contrast to the slow, irreversible thermalization driven by the cubic nonlinear KE (2) in Fig. 1.

**Experiments.** We have experimentally investigated the two different turbulent regimes that we have previously described. We use 100-ps long pulses made of temporally incoherent optical waves with a near Gaussian beam shape and the following spectro-temporal features: 558-nm central wavelength, full width at half maximum (FWHM) of about 1.93 THz [56]. The temporally incoherent pulse is then divided into two linear polarization states with relative tunable power and temporal delay (much larger than their coherence time). The incoherent waves are then injected into a 6.2-m-long weakly birefringent silica fiber ($\alpha \approx 0.565 \text{ m}^{-1}$ at 558 nm) whose propagation is governed by Eq.(1) [9]. The complex anomalous correlator $M$, as well as the optical powers $N_{x,y}$ and spectra, along $x$ and $y$ axes, are measured at the input and output of the fiber using an optical spectrum analyzer and a polarimeter.

We begin by investigating the irreversible thermalization process. This is obtained when injecting the same power along both axes $x$ and $y$, i.e., $\Delta N^0 = 0$. Fig. 3(a) shows the measurements of the power difference $\Delta N/N$ and degree of polarization $P$ at the fiber output, as a function of the input power. We observe that the degree of polarization is effectively determined by the power imbalance $P \approx \Delta N/N$ (see Eq.(3)), since the measured anomalous correlator is negligible ($|M| \approx 0$), as expected for the thermalization regime under investigation. In Fig. 3(a) we can see that, as the input power increases, the degree of polarization $P$ at the fiber output grows larger. This behavior is found to be in agreement with the numerical simulations of NLSE (1), performed by faithfully reproducing the details of the experimental conditions [56]. As predicted by Eq.(7), a sufficiently large power difference may render the anomalous correlator unstable. It is worth emphasizing, however, that even for the largest power difference recorded at the fiber output in Fig. 3(a), the anomalous correlator remains stable (which is consistent with Eq.(7), as $\alpha L_{nl} \approx 0.19 > (2/3)\Delta N/N \approx 0.12$).

Figure 3(b) shows the measured input and output spectra along the $x$-axis of the fiber (similar results are obtained

<figure><img src="image_5.png" /><figcaption>FIG. 3: Observation of optical repolarization. (a) Experimental measurements of the normalized power difference $\Delta N/N$ (blue dots), and degree of polarization $P$ (red crosses), vs input power: In the thermalization regime, the anomalous correlator $M$ is negligible, and $P \approx \Delta N/N$, see Eq.(3). The measured repolarization is in good agreement with the simulations of NLSE (1) (dashed black lines). (b) Measured power spectrum at the input (blue line), and output (orange line), of the fiber. The experimental spectra agree well with the simulation of NLSE (1) (black lines), and exhibit a power-law signature of light thermalization $\sim \omega^{-2}$ in the spectral tail.</figcaption></figure>

along the $y$-axis), together with the spectrum obtained from NLSE simulations. The experimental and numerical results show the formation of a spectral tail decaying with the power law $\sim \omega^{-2}$, which is a signature of light thermalization reflecting energy equipartition among the modes [1-5, 18] (also see Fig. S4 in [56]).

We now turn our attention to the second regime characterized by a reversible oscillatory turbulent dynamics. Constraints on the injection conditions prevent simultaneously launching a strong unbalanced population $S_1$ with a small phase-correlation $|M| \approx 0$, so that we are unable to test the criterion (7). However, aside from the instability of the fixed point $S_1 = +S_0$, the AC-KE (8) predicts a general nonlinear coupling between normal and anomalous correlators for arbitrary initial conditions on the sphere $S(z=0)$. In our experiments, we probe this dynamics starting from $S_1(z=0) \approx 0$ with a large phase-correlation $|M|$. We thus restore a large anomalous correlator by means of a polarizer, placed before injection of the laser beam into the fiber. Figs. 4(a-b) report the experimental measurements of the power difference $\Delta N/N$ and complex anomalous correlator $M/N$ at the fiber output, as functions of the input power. Their non-monotonic evolutions agree well with the simulations of the NLSE (1). Such an agreement has been obtained by using a single adjustable parameter, namely the initial relative phase between $u_x$ and $u_y$, which accounts for a minor, uncontrolled elliptical polarization of the injected incoherent optical beam. The analysis of the evolution of $\Delta N/N$ and $M/N$ along the fiber length $L = 6.2 \text{m}$ uncovers the underlying oscillatory turbulent dynamics. This is illustrated in the insets of Fig. 4(a-b) for an input power of 50W, which corresponds to a strongly nonlinear regime with $|E/U| \approx 1$. The blue and red curves refer to simulations of the NLSE (1) and the AC-KE (8), respectively, whose corresponding evolution is also reported on
<figure><img src="image_6.png" /><figcaption>FIG. 4: **Observation of the reversible turbulent regime.** Measurements, at the fiber output and as a function of the input power, of (a) the power difference $\Delta N/N$ and (b) the real and imaginary parts of the normalized anomalous correlator $M/N$. The insets in (a)-(b) show the dynamics along the fiber length $L = 6.2\,\text{m}$ for $\Delta N/N$ and $M^t/N$ with an input power of $50\,\text{W}$ ($L_{nl} \approx 0.83\,\text{m}$) in the strongly nonlinear regime $|E/U| \approx 1$: The blue and red lines report the simulations of NLSE (1) and AC-KE (8), respectively. (c) Corresponding field evolution on the Poincaré sphere for the NLSE (blue) and the AC-KE (8) (red). Note in the insets (a)-(b) and in panel (c) that the AC-KE (8) (valid in the weakly nonlinear regime $|E/U| \gg 1$) still provides a qualitative description of the oscillatory dynamics even in the strongly nonlinear regime $|E/U| \approx 1$.</figcaption></figure>

the Stokes sphere, see Fig. 4(c). These results show that the AC-KE (8), which is rigorously valid in the weakly nonlinear regime $|E/U| \ll 1$ remains robust in describing qualitatively the oscillatory turbulent behavior even in the strongly nonlinear regime $|E/U| \simeq 1$.

**Conclusion and perspectives.** We have revealed, both theoretically and experimentally, that a phase-invariant Hamiltonian system of coupled waves may exhibit two distinct turbulent regimes: irreversible thermalization and reversible dynamics mediated by phase-correlations. The criterion (7) determines which regime emerges when the initial anomalous correlator vanishes, while a non-vanishing initial anomalous correlator drives the system into the reversible turbulent regime. The experiments provide the first evidence of such contrasting turbulent behaviors in conservative (Hamiltonian) systems of dispersive nonlinear waves.

Our study establishes a basis for future advances in understanding the interplay of these fundamentally different regimes. Until now, turbulent regimes dominated by slow, irreversible thermalization and those governed by fast, correlation-driven reversible dynamics have been

separately explored, leaving open the fundamental question about the role of anomalous correlators on thermalization, and the reciprocal impact of irreversible processes on the reversible dynamics. A key challenge for future research is therefore the development of a generalized theory that unifies these two aspects – reversible dynamics associated to anomalous correlators and irreversible thermalization – fully accounting for phase-correlations.

**Acknowledgement.** Funding was provided by Agence Nationale de la Recherche (Grants No. ANR-23-CE30-0021, No. ANR-21-ESRE-0040, No. ANR-21-CE42-0026), by the Marie Skłodowska-Curie Actions project BESCLING (Grant No. 101209943), and by the Italian Ministry of University and Research (MUR) through the Italian Science Fund (FIS2), project SOFT (H53C24001530001). Calculations were performed using HPC resources from Université Côte d'Azur's Center Azzura, and DNUM-CCUB (Université de Bourgogne Europe). S.W. acknowledges support from the European Innovation Council, project MULTISCOPE (101185664).

[1] V.E. Zakharov, V.S. L'vov, G. Falkovich, *Kolmogorov Spectra of Turbulence I* (Springer, Berlin, 1992).

[2] V. Zakharov, F. Dias, and A. Pushkarev, One-Dimensional Wave Turbulence, Phys. Rep. **398**, 1 (2004).

[3] S. Nazarenko, *Wave Turbulence* (Springer, Lectures Notes in Physics, 2011).

[4] A.C. Newell, B. Rumpf, Wave Turbulence, Annu. Rev. Fluid Mech. **43**, 59 (2011).

[5] S. Galtier, *Physics of Wave Turbulence* (Cambridge University Press, 2022).

[6] M. Onorato, Y.V. Lvov, G. Dematteis, S. Chibbaro, Wave Turbulence and thermalization in one-dimensional chains, Physics Reports **1040**, 1 (2023).

[7] V.S. L'vov, Y.V. Lvov, A.C. Newell, and V.E. Zakharov, Statistical Description of Acoustic Turbulence, Phys. Rev. E **56**, 390 (1997).

[8] B. Rumpf, A.C. Newell, and V.E. Zakharov, Turbulent Transfer of Energy by Radiating Pulses, Phys. Rev. Lett. **103**, 074502 (2009).

[9] V. Labarre, G. Krstulovic, S. Nazarenko, Wave-Kinetic Dynamics of Forced-Dissipated Turbulent Internal Gravity Waves, Phys. Rev. Lett. **135**, 014101 (2025).

[10] S. Galtier, S. Nazarenko, Turbulence of Weak Gravitational Waves in the Early Universe, Phys. Rev. Lett. **119**, 221101 (2017).

[11] G. Dematteis, L. Rondoni, D. Proment, F. De Vita, M.
Onorato, Coexistence of Ballistic and Fourier Regimes in the β-Fermi-Pasta-Ulam-Tsingou Lattice, Phys. Rev. Lett. 125, 024101 (2020).

[12] M. Onorato, L. Vozella, D. Proment, Y.V. Lvov, Route to thermalization in the α-Fermi-Pasta-Ulam system, Proc. Natl. Acad. Sci. (PNAS) **112**, 4208 (2015).

[13] J. Laurie, U. Bortolozzo, S. Nazarenko, S. Residori, One-dimensional optical wave turbulence: experiment and theory, Physics Reports **514**, 121-175 (2012).

[14] S.A. Babin, D.V. Churkin, A.E. Ismagulov, S.I. Kablukov, E.V. Podivilov, Four-wave-mixing-induced turbulent spectral broadening in a long Raman fiber laser, J. Opt. Soc. Am. B **24**, 1729-1738 (2007).

[15] S.K. Turitsyn, S.A. Babin, E.G. Turitsyna, G.E. Falkovich, E. Podivilov, D. Churkin, Optical Wave Turbulence, in Advances in Wave Turbulence (World Scientific, Singapore, 2013).

[16] D. Churkin, I. Kolokolov, E. Podivilov, I. Vatnik, S. Vergeles, I. Terekhov, V. Lebedev, G. Falkovich, M. Nikulin, S. Babin, S. Turitsyn, Wave kinetics of a random fibre laser, Nature Communications **2**, 6214 (2015).

[17] P. Suret, S. Randoux, H. Jauslin, A. Picozzi, Anomalous thermalization of nonlinear wave systems, Phys. Rev. Lett. **104**, 054101 (2010).

[18] A. Picozzi, J. Garnier, T. Hansson, P. Suret, S. Randoux, G. Millot, D.N. Christodoulides, Optical wave turbulence: Toward a unified nonequilibrium thermodynamic formulation of statistical nonlinear optics, Physics Reports **542**, 1-132 (2014).

[19] F.O. Wu, A.U. Hassan, D.N. Christodoulides, Thermo-dynamic theory of highly multimoded nonlinear optical systems, Nature Photonics **13**, 776 (2019).

[20] A. Fusaro, J. Garnier, K. Krupa, G. Millot, A. Picozzi, Dramatic acceleration of wave condensation mediated by disorder in multimode fibers, Phys. Rev. Lett. **122**, 123902 (2019).

[21] E.V. Podivilov, F. Mangini, O.S. Sidelnikov, M. Ferraro, M. Gervaziev, D.S. Kharenko, M. Zitelli, M.P. Fedoruk, S.A. Babin, S. Wabnitz, Thermalization of orbital angular momentum beams in multimode optical fibers, Phys. Rev. Lett. **128**, 243901 (2022).

[22] K. Baudin, A. Fusaro, K. Krupa, J. Garnier, S. Rica, G. Millot, A. Picozzi, Classical Rayleigh-Jeans condensation of light waves: Observation and thermodynamic characterization, Phys. Rev. Lett. **125**, 244101 (2020).

[23] H. Pourbeyram, P. Sidorenko, F. Wu, L. Wright, D. Christodoulides, F. Wise, Direct observations of thermalization to a Rayleigh-Jeans distribution in multimode optical fibres, Nature Physics **18**, 685 (2022).

[24] F. Mangini, M. Gervaziev, M. Ferraro, D. S. Kharenko, M. Zitelli, Y. Sun, V. Couderc, E.V. Podivilov, S.A. Babin, and S. Wabnitz, Statistical mechanics of beam self-cleaning in GRIN multimode optical fibers, Opt. Exp. **30**, 10850 (2022).

[25] M. Ferraro, F. Mangini, F.O. Wu, M. Zitelli, D.N. Christodoulides, S. Wabnitz, Calorimetry of photon gases in nonlinear multimode optical fibers, Phys. Rev. X **14**, 021020 (2024).

[26] A.L. Muniz, F. Wu, P. Jung, M. Khajavikhan, D. Christodoulides, U. Peschel, Observation of photon-photon thermodynamic processes under negative optical temperature conditions, Science **379**, 1019 (2023).

[27] M.S. Kirsch, G.G. Pyrialakos, R. Altenkirch, M.A. Selim, J. Beck, T.A.W. Wolterink, H. Ren, P.S. Jung, M. Kha-

javikhan, A. Szameit, M. Heinrich, D.N. Christodoulides, Observation of Joule-Thomson photon-gas expansion, Nat. Phys. **21**, 214-220 (2025).

[28] N. Berti, K. Baudin, A. Fusaro, G. Millot, A. Picozzi, J. Garnier, Interplay of thermalization and strong disorder: Wave turbulence theory, numerical simulations, and experiments in multimode optical fibers, Phys. Rev. Lett. **129**, 063901 (2022).

[29] K. Baudin, J. Garnier, A. Fusaro, N. Berti, C. Michel, K. Krupa, G. Millot, A. Picozzi, Observation of light thermalization to negative temperature Rayleigh-Jeans equilibrium states in multimode optical fibers, Phys. Rev. Lett. **130**, 063801 (2023).

[30] M. Lian, Y. Geng, Y.-J. Chen, Y. Chen, J.-T. Lü, Coupled Thermal and Power Transport of Optical Waveguide Arrays: Photonic Wiedemann-Franz Law and Rectification Effect, Phys. Rev. Lett. **133**, 116303 (2024).

[31] A. Ramos, L. Fernández-Alcázar, T. Kottos, B. Shapiro, Optical Phase Transitions in Photonic Networks: a Spin-System Formulation Phys. Rev. X **10**, 031024 (2020).

[32] A.Y. Ramos, C. Shi, L.J. Fernandez-Alcazar, D.N. Christodoulides, T. Kottos, Theory of localization-hindered thermalization in nonlinear multimode photonics, Communications Physics **6**, 189 (2023).

[33] A. Kurnosov, L.J. Fernández-Alcázar, A. Ramos, B. Shapiro, T. Kottos, Optical Kinetic Theory of Nonlinear Multimode Photonic Networks, Phys. Rev. Lett. **132**, 193802 (2024).

[34] G. Yang, D. Bongiovanni, D. Song, R. Morandotti, Z. Chen, N.K. Efremidis, Thermalization dynamics in photonic lattices of different geometries, APL Photonics **9**, 066115 (2024).

[35] M. Ferraro, K. Baudin, M. Gervaziev, A. Fusaro, A. Picozzi, J. Garnier, G. Millot, D. Kharenko, E. Podivilov, S. Babin, F. Mangini, S. Wabnitz, Wave turbulence, thermalization and multimode locking in optical fibers, Physica D **481**, 134758 (2025).

[36] Y. Deng, Z. Hani, Full derivation of the wave kinetic equation, Invent. math. **233**, 543-724 (2023).

[37] Y. Deng and Z. Hani, Long time justification of wave turbulence theory, arXiv:2311.10082 (2023).

[38] G. Staffilani and M.-B. Tran, On the wave turbulence theory for stochastic and random multidimensional KdV type equations, arXiv:2106.09819 (2021).

[39] E. Turitsyna, G. Falkovich, A. El-Taher, X. Shu, P. Harper, and S. Turitsyn, Optical turbulence and spectral condensate in long fibre lasers, Proc. R. Soc. London Ser. A **468**, 2145 (2012).

[40] E. Turitsyna, S. Smirnov, S. Sugavanam, N. Tarasov, X. Shu, S. Babin, E. Podivilov, D. Churkin, G. Falkovich, and S. Turitsyn, The Laminar-Turbulent Transition in a Fibre Laser, Nat. Photonics **7**, 783 (2013).

[41] C. Rodda, C. Savaro, V. Bouillaut, P. Augier, J. Sommeria, T. Valran, S. Viboud, N. Mordant, From Internal Waves to Turbulence in a Stably Stratified Fluid, Phys. Rev. Lett. **131**, 264101 (2023).

[42] D. Pierangeli, F. Di Mei, G. Di Domenico, A. J. Agranat, C. Conti, and E. DelRe, Turbulent Transitions in Optical Wave Propagation, Phys. Rev. Lett. **117**, 183902 (2016).

[43] L. Dieli, D. Pierangeli, E. DelRe, C. Conti, Observation of Two-Dimensional Dam Break Flow and a Gaseous Phase of Solitons in a Photon Fluid, Phys. Rev. Lett. **133**, 183801 (2024).

[44] J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Theory of
superconductivity, Phys. Rev. **108**, 1175 (1957).

[45] V.E. Zakharov, V.S. L'Vov, and S.S. Starobinets, Stationary nonlinear theory of parametric excitation of waves, Sovient Physics JETP **32**, 656 (1971).

[46] V.E. Zakharov, V.S. Lvov, and S.S. Starobinets, Spin-wave turbulence beyond the parametric excitation threshold, Sov. Phys. Usp. **17**, 896 (1975).

[47] V.S. L'vov, Wave Turbulence Under Parametric Excitation – Application to Magnets (Springer, Berlin, 1994).

[48] J. Zaleski, M. Onorato, Y.V. Lvov, Anomalous Correlators in Nonlinear Dispersive Wave Systems, Phys. Rev. X **10**, 021043 (2020).

[49] M. Guasoni, J. Garnier, B. Rumpf, D. Sugny, J. Fatome, F. Amrani, G. Millot, A. Picozzi, Incoherent Fermi-Pasta-Ulam Recurrences and Unconstrained Thermalization Mediated by Strong Phase Correlations, Phys. Rev. X **7**, 011025 (2017).

[50] A. Villois, G. Dematteis, Y.V. Lvov, M. Onorato, and J. Shatah, Anomalous correlators, negative frequencies and non-phase-invariant Hamiltonians in random waves, Phys. Rev. Lett. **135**, 167201 (2025). Considering a scalar wave equation, it was proven in this work that phase invariant systems cannot develop anomalous correlations. The theorem can be extended to vector equations such as NLSE (1), for which phase invariance implies that, for given initial data $u_0$, the solution $u(z) = F(z, u_0)$, verifies $e^{i\theta}F(z, u_0) = F(z, e^{i\theta}u_0)$, $\theta \in \mathbb{R}$. Consequently, neither the $u_x$ component nor the $u_y$ component can individually develop anomalous correlations in the NLSE (1).

[51] L. Pitaevskii and S. Stringari, *Bose-Einstein Condensation and Superfluidity* (Oxford University Press, 2016).

[52] R. Cominotti, A. Berti, C. Dulin, C. Rogora, G. Lamporesi, I. Carusotto, A. Recati, A. Zenesini, and G. Ferrari, Ferromagnetism in an Extended Coherently Coupled Atomic Superfluid, Phys. Rev. X **13**, 021037 (2023).

[53] A. Zenesini, A. Berti, R. Cominotti, C. Rogora, I.G. Moss, T.P. Billam, I. Carusotto, G. Lamporesi, A. Recati, G. Ferrari, False vacuum decay via bubble formation in ferromagnetic superfluids, Nature Physics **20**, 558 (2024).

[54] G. Millot, E. Seve, S. Wabnitz Polarization Symmetry Breaking and Pulse Train Generation from the Modulation of Light Waves, Phys. Rev. Lett. **79**, 661 (1997).

[55] G. Agrawal, Nonlinear Fiber Optics (sixth Ed., Academic Press, New York, 2019).

[56] The Supplementary Material presents the derivation of a reduced form of the WT KE (2), as well as the derivation of the AC-KE (4-5)-(8) and the corresponding stability analysis, Eq.(6). It also provides technical details regarding the experimental procedures and measurements, as well as additional numerical simulations.

[57] The definition of the degree of polarization in Eq.(3) is equivalent to the usual one: $P(z) = (1 - 4\det(J)/\text{tr}(J)^2)^{1/2}$, where $J_{ij} = \langle u_i u_j^* \rangle$ is the 2×2 coherence matrix ($j = x, y$), with ‘det’ and ‘tr’ the determinant and trace operations, see, e.g., L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, New York, 1995).

[58] B. Daino, G. Gregori, S. Wabnitz, Stability analysis of nonlinear coherent coupling, Journal of Applied Physics **58**, 4512 (1985).
Supplementary information on the article:

I. THERMALIZATION: DERIVATION OF A REDUCED WAVE TURBULENCE KINETIC EQUATION FOR SIMULATIONS

In Fig. 1 (main text), we report numerical simulations of the WT KE (2) (main text). The one-dimensional interaction allows the KE (2) (main text) to be written in a simplified form, which greatly improves the efficiency of its numerical integration. In the following, we derive the simplified version of KE (2) (main text).

We start from the coherently coupled NLSE (1) (main text):

$$i\partial_z u_x = -\beta\partial_{tt}u_x + \alpha u_x + \gamma(|u_x|^2 + \kappa|u_y|^2)u_x + \gamma\rho u_x^* u_y^*, \quad (10)$$

$$i\partial_z u_y = -\beta\partial_{tt}u_y - \alpha u_y + \gamma(|u_y|^2 + \kappa|u_x|^2)u_y + \gamma\rho u_y^* u_x^*. \quad (11)$$

The corresponding dispersion relations read $k_x(\omega) = \beta\omega^2 + \alpha$ and $k_y(\omega) = \beta\omega^2 - \alpha$. The NLSE conserves the power (particle number) $N = N_x + N_y$, with $N_\mu = \frac{1}{T}\int_0^T |u_\mu|^2 dt$ ($\mu = x, y$), where $T$ denotes the size of the numerical window. It also conserves the Hamiltonian $H = E + U$, which has a linear and a nonlinear contribution. The linear contribution $E = E_k + E_c$ can be split into a kinetic contribution and a resonant coupling contribution:

$$E_k = \frac{1}{T} \int_0^T \beta (|\partial_t u_x|^2 + |\partial_t u_y|^2) dt, \quad E_c = \frac{1}{T} \int_0^T \alpha (|u_x|^2 - |u_y|^2) dt. \quad (12)$$

The nonlinear contribution reads

$$U = \frac{1}{T} \int_0^T \frac{\gamma}{2} (|u_x|^4 + |u_y|^4) + \gamma \kappa |u_x|^2 |u_y|^2 + \frac{\gamma \rho}{2} (u_x^{*2} u_y^2 + u_x^2 u_y^{*2}) dt. \quad (13)$$

Following the standard wave turbulence theory [1–3], we derive from the vector NLSE (10-11) the following kinetic equations

$$\partial_z n_x(\omega, z) = \kappa^2 \mathrm{Coll}_i[n_x, n_y] + \rho^2 \mathrm{Coll}_c[n_x, n_y], \quad (14)$$

$$\partial_z n_y(\omega, z) = \kappa^2 \mathrm{Coll}_i[n_y, n_x] + \rho^2 \mathrm{Coll}_c[n_y, n_x], \quad (15)$$

where $\langle \tilde{u}_\mu(\omega_1, z) \tilde{u}_\mu^*(\omega_2, z) \rangle = 2\pi n_\mu(\omega_1, z) \delta(\omega_1 - \omega_2)$ under the assumption of homogeneous statistics, where $\tilde{u}_\mu(\omega, z) = \int u_\mu(t, z) \exp(-i\omega t) dt$ is the Fourier transform, for $\mu = x, y$. The collision term for the non-resonant interaction reads:

$$\mathrm{Coll}_i[n_x, n_y] = \frac{\gamma^2}{2\pi} \iiint d\omega_{1-3} n_x(\omega_1) n_y(\omega_2) n_y(\omega_3) n_x(\omega) (n_x^{-1}(\omega) + n_y^{-1}(\omega) - n_y^{-1}(\omega_2) - n_x^{-1}(\omega_1)) \\ \times \delta(k_x(\omega_1) + k_y(\omega_2) - k_y(\omega_3) - k_x(\omega)) \delta(\omega_1 + \omega_2 - \omega_3 - \omega), \quad (16)$$

$$\mathrm{Coll}_i[n_y, n_x] = \frac{\gamma^2}{2\pi} \iiint d\omega_{1-3} n_y(\omega_1) n_x(\omega_2) n_x(\omega_3) n_y(\omega) (n_y^{-1}(\omega) + n_x^{-1}(\omega) - n_x^{-1}(\omega_2) - n_y^{-1}(\omega_1)) \\ \times \delta(k_y(\omega_1) + k_x(\omega_2) - k_x(\omega_3) - k_y(\omega)) \delta(\omega_1 + \omega_2 - \omega_3 - \omega), \quad (17)$$

with $d\omega_{1-3} = d\omega_1 d\omega_2 d\omega_3$. The collision term for the resonant interaction reads:

$$\mathrm{Coll}_c[n_x, n_y] = \frac{2\gamma^2}{2\pi} \iiint d\omega_{1-3} n_y(\omega_1) n_y(\omega_2) n_x(\omega_3) n_x(\omega) (n_x^{-1}(\omega) + n_x^{-1}(\omega_3) - n_y^{-1}(\omega_2) - n_y^{-1}(\omega_1)) \\ \times \delta(k_y(\omega_1) + k_y(\omega_2) - k_x(\omega_3) - k_x(\omega)) \delta(\omega_1 + \omega_2 - \omega_3 - \omega), \quad (18)$$

$$\mathrm{Coll}_c[n_y, n_x] = \frac{2\gamma^2}{2\pi} \iiint d\omega_{1-3} n_x(\omega_1) n_x(\omega_2) n_y(\omega_3) n_y(\omega) (n_y^{-1}(\omega) + n_y^{-1}(\omega_3) - n_x^{-1}(\omega_2) - n_x^{-1}(\omega_1)) \\ \times \delta(k_x(\omega_1) + k_x(\omega_2) - k_y(\omega_3) - k_y(\omega)) \delta(\omega_1 + \omega_2 - \omega_3 - \omega). \quad (19)$$

Note that, because of the degenerate resonances inherent to the purely one-dimensional problem considered here, the collision term describing the self-interaction (first nonlinear term in Eqs.(10-11)) vanishes identically, see [1–3]. We focus on the analysis of the resonant collision term. Two integrals can be computed owing to the Dirac $\delta$-functions,
which gives:

$$
\begin{equation}
\begin{split}
\partial_z n_x(\omega, z) = \kappa^2 \mathrm{Coll}_i[n_x, n_y] &+ \frac{2\rho^2\gamma^2}{4\pi|\beta|} \int n_y(\omega + \omega') n_y(\omega - \nu/\omega') n_x(\omega + \omega' - \nu/\omega') n_x(\omega) \\
& \times (n_x^{-1}(\omega) + n_x^{-1}(\omega + \omega' - \nu/\omega') - n_y^{-1}(\omega + \omega') - n_y^{-1}(\omega - \nu/\omega')) \frac{d\omega'}{|\omega'|}
\end{split}
\tag{20}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\partial_z n_y(\omega, z) = \kappa^2 \mathrm{Coll}_i[n_y, n_x] &+ \frac{2\rho^2\gamma^2}{4\pi|\beta|} \int n_x(\omega + \omega') n_x(\omega + \nu/\omega') n_y(\omega + \omega' + \nu/\omega') n_y(\omega) \\
& \times (n_y^{-1}(\omega) + n_y^{-1}(\omega + \omega' + \nu/\omega') - n_x^{-1}(\omega + \omega') - n_x^{-1}(\omega + \nu/\omega')) \frac{d\omega'}{|\omega'|}
\end{split}
\tag{21}
\end{equation}
$$

where $\nu = 2\alpha/\beta$.

$$
\begin{equation}
\begin{split}
\partial_z n_x(\omega, z) ={}& \frac{\kappa^2 \gamma^2}{4\pi|\beta|} \int n_x(\omega') n_y(\omega) n_y(\omega') n_x(\omega) (n_x^{-1}(\omega) + n_y^{-1}(\omega') - n_y^{-1}(\omega) - n_x^{-1}(\omega')) \frac{d\omega'}{|\omega' - \omega|} \\
& + \frac{2\rho^2 \gamma^2}{4\pi|\beta|} \int n_y(\omega + \sqrt{|\nu|\eta}) n_y(\omega - \sqrt{|\nu|\sigma/\eta}) n_x(\omega + \sqrt{|\nu|\eta} - \sqrt{|\nu|\sigma/\eta}) n_x(\omega) \\
& \times (n_x^{-1}(\omega) + n_x^{-1}(\omega + \sqrt{|\nu|\eta} - \sqrt{|\nu|\sigma/\eta}) - n_y^{-1}(\omega + \sqrt{|\nu|\eta} - n_y^{-1}(\omega - \sqrt{|\nu|\sigma/\eta})) \frac{d\eta}{|\eta|},
\end{split}
\tag{22}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\partial_z n_y(\omega, z) ={}& \frac{\kappa^2 \gamma^2}{4\pi|\beta|} \int n_x(\omega') n_y(\omega) n_y(\omega') n_x(\omega) (n_x^{-1}(\omega) + n_y^{-1}(\omega') - n_y^{-1}(\omega) - n_x^{-1}(\omega')) \frac{d\omega'}{|\omega' - \omega|} \\
& + \frac{2\rho^2 \gamma^2}{4\pi|\beta|} \int n_x(\omega + \sqrt{|\nu|\eta}) n_x(\omega + \sqrt{|\nu|\sigma/\eta}) n_y(\omega + \sqrt{|\nu|\eta} + \sqrt{|\nu|\sigma/\eta}) n_y(\omega) \\
& \times (n_y^{-1}(\omega) + n_y^{-1}(\omega + \sqrt{|\nu|\eta} + \sqrt{|\nu|\sigma/\eta}) - n_x^{-1}(\omega + \sqrt{|\nu|\eta} - n_x^{-1}(\omega + \sqrt{|\nu|\sigma/\eta})) \frac{d\eta}{|\eta|},
\end{split}
\tag{23}
\end{equation}
$$

with $\sigma = \operatorname{sgn}(\nu)$.

The singularity in the resonant collision terms (proportional to $\rho^2$) can be rewritten in a convenient form by the
change of variable $\eta \to se^\theta$, $s = \pm 1$, so that the kinetic equations take the form

$$
\begin{equation}
\begin{split}
\partial_z n_x(\omega, z) ={}& \frac{\kappa^2 \gamma^2}{4\pi|\beta|} \int n_x(\omega') n_y(\omega) n_y(\omega') n_x(\omega) (n_x^{-1}(\omega) + n_y^{-1}(\omega') - n_y^{-1}(\omega) - n_x^{-1}(\omega')) \frac{d\omega'}{|\omega' - \omega|} \\
& + \frac{2\rho^2 \gamma^2}{4\pi|\beta|} \sum_s \int n_y(\omega + \sqrt{|\nu|}se^\theta) n_y(\omega - \sqrt{|\nu|}\sigma se^{-\theta}) n_x(\omega + \sqrt{|\nu|}se^\theta - \sqrt{|\nu|}\sigma se^{-\theta}) n_x(\omega) \\
& \times (n_x^{-1}(\omega) + n_x^{-1}(\omega + \sqrt{|\nu|}se^\theta - \sqrt{|\nu|}\sigma se^{-\theta}) - n_y^{-1}(\omega + \sqrt{|\nu|}se^\theta - n_y^{-1}(\omega - \sqrt{|\nu|}\sigma se^{-\theta})) d\theta,
\end{split}
\tag{24}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\partial_z n_y(\omega, z) ={}& \frac{\kappa^2 \gamma^2}{4\pi|\beta|} \int n_y(\omega') n_x(\omega) n_x(\omega') n_y(\omega) (n_y^{-1}(\omega) + n_x^{-1}(\omega') - n_x^{-1}(\omega) - n_y^{-1}(\omega')) \frac{d\omega'}{|\omega' - \omega|} \\
& + \frac{2\rho^2 \gamma^2}{4\pi|\beta|} \sum_s \int n_x(\omega + \sqrt{|\nu|}se^\theta) n_x(\omega + \sqrt{|\nu|}\sigma se^{-\theta}) n_y(\omega + \sqrt{|\nu|}se^\theta + \sqrt{|\nu|}\sigma se^{-\theta}) n_y(\omega) \\
& \times (n_y^{-1}(\omega) + n_y^{-1}(\omega + \sqrt{|\nu|}se^\theta + \sqrt{|\nu|}\sigma se^{-\theta}) - n_x^{-1}(\omega + \sqrt{|\nu|}se^\theta - n_x^{-1}(\omega + \sqrt{|\nu|}\sigma se^{-\theta})) d\theta.
\end{split}
\tag{25}
\end{equation}
$$

The kinetic equations conserve the total power (density) $N = \frac{1}{2\pi} \sum_{\mu} \int n_{\mu}(\omega, z) d\omega$, with $\mu = x, y$. Note that if $\sigma = 1$, then

$$
\begin{equation}
\begin{split}
\partial_z N_x(z) ={}& \frac{2\rho^2\gamma^2}{2\pi|\beta|} \int d\omega \int d\theta n_y(\omega + \sqrt{|\nu|} \cosh\theta) n_y(\omega - \sqrt{|\nu|} \cosh\theta) n_x(\omega + \sqrt{|\nu|} \sinh\theta) n_x(\omega - \sqrt{|\nu|} \sinh\theta) \\
& \times (n_x^{-1}(\omega + \sqrt{|\nu|} \sinh\theta) + n_x^{-1}(\omega - \sqrt{|\nu|} \sinh\theta) - n_y^{-1}(\omega + \sqrt{|\nu|} \cosh\theta) - n_y^{-1}(\omega - \sqrt{|\nu|} \cosh\theta)) ,
\end{split}
\tag{26}
\end{equation}
$$

$$
\partial_z N_y(z) = -\partial_z N_x(z). \tag{27}
$$

If $\sigma = -1$, the same result is obtained by changing $\cosh(\theta)$ and $\sinh(\theta)$. The kinetic equation also conserves the linear energy $E = \frac{1}{2\pi} \sum_{\mu} \int k_{\mu}(\omega) n_{\mu}(\omega, z) d\omega$, with $\mu = x, y$. The reduced form of the wave turbulence KE (24-25)
<figure><img src="image_7.png" /><figcaption>FIG. 5: Irreversible thermalization. Numerical simulations of the NLSE (10-11) (blue and light-blue lines), and the wave turbulence KE (24-25) (red lines): (a) Evolution of the kinetic $E_k(z)$, and coherent $E_c(z)$, contribution to the energy $E = E_k(z) + E_c(z)$ in Eq.(12), the inset shows the corresponding evolution of the entropy $S(z)$. Spectra of the waves $n_x(\omega)$ (b), and $n_y(\omega)$ (c) at $z = 3000L_{nl}$. The thermalization process is characterized by a growth of the kinetic energy $E_k(z)$, which entails a decrease of $E_c(z)$: The irreversible transfer of power (particles) from $N_x$ to $N_y$ entails the repolarization effect (increase of $P(z)$), see Fig. 1 (main text). An average over 5 realizations of the initial condition (i.e., initial random spectral phases) has been considered for the NLSE simulations.</figcaption></figure>

has been numerically integrated. The results have been compared with the simulations of the NLSE (10-11). A good agreement is obtained without using adjustable parameters, see Fig. 1 (main text), and Fig. 5.

Unlike conventional nonlinear systems where the kinetic energy $E_k$ is the sole contribution to the linear energy $E$, our system includes an additional contribution arising from the coherent coupling, $E = E_k + E_c$ in Eq.(12). Consequently, the growth in $E_k$ is no longer constrained, as it can be offset by a corresponding decrease in $E_c$, see Fig. 5(a) [4]. Recalling that $E_c(z) = \alpha(N_x - N_y)$, thermalization manifests itself by an irreversible transfer of power (particles) from $N_x$ to $N_y$, thus leading to the emergence of a nonvanishing degree of polarization.

## II. ANOMALOUS CORRELATOR KINETIC EQUATION

### A. Derivation of Eqs.(4-5) (main text)

We start from the NLSE (10-11) and define the normal correlator $U_\mu(t_1, t_2, z) = \langle u_\mu(t_1, z)u_\mu(t_2, z)^* \rangle$ for $\mu = x, y$, and the anomalous correlator describing a correlation among the waves $u_x$ and $u_y$: $W(t_1, t_2, z) = \langle u_y(t_1, z)u_x^*(t_2, z) \rangle$. Following the usual procedure, we derive the equations for the correlators $U_\mu$ and $W$, which involve fourth-order moments of the waves that can be expanded into products of second-order moments in the weakly nonlinear regime of Gaussian statistics. Here, caution should be exercised in keeping all anomalous correlations between the waves. We obtain

$$ 
\begin{aligned}
i \partial_z U_x^{12} &= -\beta(\partial_{t_1}^2 - \partial_{t_2}^2)U_x^{12} + \gamma U_x^{12}[2(U_x^{11} - U_x^{22}) + \kappa(U_y^{11} - U_y^{22})] + \gamma \kappa [W^{12}(W^{11} + W^{11*}) - W^{21*}(W^{22} + W^{22*})], \\
i \partial_z U_y^{12} &= -\beta(\partial_{t_1}^2 - \partial_{t_2}^2)U_y^{12} + \gamma U_y^{12}[2(U_y^{11} - U_y^{22}) + \kappa(U_x^{11} - U_x^{22})] + \gamma \kappa [W^{21*}(W^{11} + W^{11*}) - W^{12}(W^{22} + W^{22*})], \\
i \partial_z W^{12} &= -\beta(\partial_{t_1}^2 - \partial_{t_2}^2)W^{12} - 2\alpha W^{12} + \gamma W^{12}[2(U_y^{11} - U_x^{22}) + \kappa(U_x^{11} - U_y^{22})] \\
&\quad + \gamma \kappa [U_x^{12}(W^{11} + W^{11*}) - U_y^{12}(W^{22} + W^{22*})],
\end{aligned}
$$

where $U_\mu^{ij} = U_\mu(t_i, t_j, z)$ for $\mu = x, y$, and $W^{ij} = W(t_i, t_j, z)$. Introducing the variables $t = (t_1 + t_2)/2$, $\tau = t_1 - t_2$, this also reads

$$ i \partial_z U_x = -2\beta \partial_{t\tau} U_x + \gamma \mathcal{L}^{U_x}, \qquad (28) $$

$$ i \partial_z U_y = -2\beta \partial_{t\tau} U_y + \gamma \mathcal{L}^{U_y}, \qquad (29) $$

$$ i \partial_z W = -2\beta \partial_{t\tau} W - 2\alpha W + \gamma \mathcal{L}^{W}, \qquad (30) $$
with the three nonlinear terms

$$ \mathcal{L}^{U_x}(t, \tau) = U_x(t, \tau)[(2U_x + \kappa U_y)(t + \tau/2, 0) - (2U_x + \kappa U_y)(t - \tau/2, 0)] \\ + 2\kappa [W(t, \tau)W^r(t + \tau/2, 0) - W(t, -\tau)^*W^r(t - \tau/2, 0)], \quad (31) $$

$$ \mathcal{L}^{U_y}(t, \tau) = U_y(t, \tau)[(2U_y + \kappa U_x)(t + \tau/2, 0) - (2U_y + \kappa U_x)(t - \tau/2, 0)] \\ + 2\kappa [W(t, -\tau)^*W^r(t + \tau/2, 0) - W(t, \tau)W^r(t - \tau/2, 0)], \quad (32) $$

$$ \mathcal{L}^{W}(t, \tau) = W(t, \tau)[(2U_y + \kappa U_x)(t + \tau/2, 0) - (2U_x + \kappa U_y)(t - \tau/2, 0)] \\ + 2\kappa [U_x(t, \tau)W^r(t + \tau/2) - U_y(t, \tau)W^r(t - \tau/2, 0)], \quad (33) $$

where $W^r = \text{Re}(W)$. We can make use of the assumption of quasi-stationary statistics, that is to say the time correlation $t_c$ (i.e. the microscopic time scale of random fluctuations) is much smaller than the macroscopic time scale $t_s$ of stationary (homogeneous) statistics. In the regime $\varepsilon = t_c/t_s \ll 1$, the three nonlinear terms can be expanded up to first order as

$$ \mathcal{L}^{U_x}(t, \tau) = U_x(t, \tau)\tau\partial_t(2U_x + \kappa U_y)(t, 0) + 2\kappa W^r(t, 0)[W(t, \tau) - W(t, -\tau)^*] + \kappa\tau\partial_t W^r(t, 0)[W(t, \tau) + W(t, -\tau)^*], \quad (34) $$

$$ \mathcal{L}^{U_y}(t, \tau) = U_y(t, \tau)\tau\partial_t(2U_y + \kappa U_x)(t, 0) - 2\kappa W^r(t, 0)[W(t, \tau) - W(t, -\tau)^*] + \kappa\tau\partial_t W^r(t, 0)[W(t, \tau) + W(t, -\tau)^*], \quad (35) $$

$$ \mathcal{L}^{W}(t, \tau) = (2 - \kappa)W(t, \tau)(U_y - U_x)(t, 0) + (2 + \kappa)W(t, \tau)\frac{\tau}{2}\partial_t(U_y + U_x)(t, 0) \\ + 2\kappa W^r(t, 0)(U_x - U_y)(t, \tau) + \kappa\tau\partial_t W^r(t, 0)(U_x + U_y)(t, \tau). \quad (36) $$

This gives

$$ i\partial_z U_x(t, \tau, z) = -2\beta\partial_{t\tau}U_x + \gamma\tau U_x\partial_t(2N_x(t) + \kappa N_y(t)) - 2\gamma\kappa W_0^r(t)(W^*(t, -\tau) - W(t, \tau)) \\ + \gamma\kappa\tau\partial_t W_0^r(t)[W(t, \tau) + W(t, -\tau)^*], \quad (37) $$

$$ i\partial_z U_y(t, \tau, z) = -2\beta\partial_{t\tau}U_y + \gamma\tau U_y\partial_t(2N_y(t) + \kappa N_x(t)) - 2\gamma\kappa W_0^r(t)(W(t, \tau) - W^*(t, -\tau)) \\ + \gamma\kappa\tau\partial_t W_0^r(t)[W(t, \tau) + W(t, -\tau)^*], \quad (38) $$

$$ i\partial_z W(t, \tau, z) = -2\beta\partial_{t\tau}W + W(\gamma(2-\kappa)(N_y(t) - N_x(t)) - 2\alpha) - 2\gamma\kappa(U_y - U_x)W_0^r(t) \\ + \gamma(1+\kappa/2)\tau W\partial_t(N_y(t) + N_x(t)) + \gamma\kappa\tau(U_y + U_x)\partial_t W_0^r(t), \quad (39) $$

where $W_0(t, z) = W(t, \tau = 0, z)$, and $N_\mu(t, z) = U_\mu(t, \tau = 0, z)$.

The equations for the normal correlator $n_\mu(t, \omega, z) = \int U_\mu(t, \tau, z) \exp(-i\omega\tau)d\tau$, and the anomalous correlator $m(t, \omega, z) = \int W(t, \tau, z) \exp(-i\omega\tau)d\tau$, are obtained by taking the Fourier transform of Eqs.(37-39), leading to the set of Eqs.(4-5) (main text):

$$ \partial_z n_x(\omega, t) = -2\beta\omega\partial_t n_x + \gamma\partial_t(2N_x(t) + \kappa N_y(t))\partial_\omega n_x + 2\gamma\kappa\partial_t M^r(t)\partial_\omega m^r(\omega, t) + 4\gamma\kappa M^r(t)m^i(\omega, t), \quad (40) $$

$$ \partial_z n_y(\omega, t) = -2\beta\omega\partial_t n_y + \gamma\partial_t(2N_y(t) + \kappa N_x(t))\partial_\omega n_y + 2\gamma\kappa\partial_t M^r(t)\partial_\omega m^r(\omega, t) - 4\gamma\kappa M^r(t)m^i(\omega, t), \quad (41) $$

$$ \partial_z m(\omega, t) = -2\beta\omega\partial_t m + \gamma(1 + \kappa/2)\partial_t(N_x(t) + N_y(t))\partial_\omega m + \gamma\kappa\partial_t M^r(t)\partial_\omega(n_x + n_y) \\ - im(\gamma(2 - \kappa)(N_y(t) - N_x(t)) - 2\alpha) - 2i\gamma\kappa(n_x - n_y)M^r(t), \quad (42) $$

where $M(t, z) = W_0(t, z) = \frac{1}{2\pi} \int m(t, \omega, z)d\omega$, and $N_\mu(t, z) = \frac{1}{2\pi} \int n_\mu(t, \omega, z)d\omega$, while $m^{r,i}$ ($M^{r,i}$) denote the real and imaginary parts of $m$ ($M$).

Note that the system of Eqs.(40-42) is formally reversible in the 'time' variable $z$, i.e., it is invariant under the transformation: ($z \to -z, \omega \to -\omega, m \to m^*$).

## B. Stability analysis for the emergence of phase-correlations: Derivation of Eq.(6) (main text)

We assume that the initial random waves $u_\mu(t, z = 0)$ have a homogeneous statistics with spectra $n_\mu^0(\omega)$, and are uncorrelated with each other, so that the initial anomalous correlator is zero, $m(\omega, t, z = 0) = 0$. We perform a linear stability analysis of Eqs.(40-42) around this solution with

$$ n_\mu(t, \omega, z) = n_\mu^0(\omega) + \delta n_\mu(t, \omega, z), \quad \mu = x, y, \quad (43) $$

$$ m(t, \omega, z) = \delta m(t, \omega, z), \quad (44) $$
with $\delta n_\mu(t, \omega, z) \ll n_\mu^0(\omega)$ and $|\delta m(t, \omega, z)| \ll n_\mu^0(\omega)$. Introducing (43-44) into Eqs.(40-42), we obtain after linearization

$$ \partial_z \delta n_x(t, \omega) = -2\beta\omega \partial_t \delta n_x + \gamma \partial_\omega n_x^0 \partial_t (2\delta N_x(t) + \kappa \delta N_y(t)), \quad (45) $$

$$ \partial_z \delta n_y(t, \omega) = -2\beta\omega \partial_t \delta n_y + \gamma \partial_\omega n_y^0 \partial_t (2\delta N_y(t) + \kappa \delta N_x(t)), \quad (46) $$

$$ \begin{aligned} \partial_z \delta m(t, \omega) = & -2\beta\omega \partial_t \delta m - i\delta m (\gamma(2-\kappa)(N_y^0 - N_x^0) - 2\alpha) - 2i\gamma\kappa(n_x^0 - n_y^0)\delta M^r(t) \\ & + \gamma\kappa \partial_\omega (n_x^0 + n_y^0) \partial_t \delta M^r(t), \end{aligned} \quad (47) $$

where $N_\mu^0 = \frac{1}{2\pi} \int n_\mu^0(\omega) d\omega$, $\delta N_\mu(t, z) = \frac{1}{2\pi} \int \delta n_\mu(t, \omega, z) d\omega$, $\delta M(t, z) = \frac{1}{2\pi} \int \delta m(t, \omega, z) d\omega$. Note that the evolution of $\delta n_\mu$ in (45-46) is decoupled from the evolution $\delta m$ in (47). In the following we study separately the stability analysis for the anomalous correlator, and the normal correlators.

1. Stability analysis of the anomalous correlator, Eq.(47)

We perform the stability analysis of the anomalous correlator by introducing $\delta \tilde{m}(\Omega, \omega, z) = \int \delta m(t, \omega, z) \exp(-i\Omega t) dt$, we obtain from Eq.(42)

$$ i \partial_z \delta \tilde{m}(\Omega, \omega, z) = (2\beta\omega\Omega + \gamma(2-\kappa)(N_y^0 - N_x^0) - 2\alpha)\delta \tilde{m} \\ + \frac{\gamma\kappa}{2\pi}(n_x^0(\omega) - n_y^0(\omega) - \frac{\Omega}{2}\partial_\omega n_y^0(\omega) - \frac{\Omega}{2}\partial_\omega n_x^0(\omega)) \int (\delta \tilde{m}(\Omega, \omega_1, z) + \delta \tilde{m}^*(-\Omega, \omega_1, z)) d\omega_1. $$

To derive the dispersion relation, it proves convenient to introduce the function $\delta \tilde{p}(\Omega, \omega, z) = \delta \tilde{m}^*(\Omega, -\omega, z)$, which gives the system

$$ \begin{aligned} i \partial_z \delta \tilde{m}(\Omega, \omega, z) = & (2\beta\omega\Omega + \gamma(2-\kappa)(N_y^0 - N_x^0) - 2\alpha)\delta \tilde{m} \\ & + \frac{\gamma\kappa}{2\pi}(n_x^0(\omega) - n_y^0(\omega) - \frac{\Omega}{2}\partial_\omega n_y^0(\omega) - \frac{\Omega}{2}\partial_\omega n_x^0(\omega)) \int (\delta \tilde{m}(\Omega, \omega_1, z) + \delta \tilde{p}(\Omega, \omega_1, z)) d\omega_1, \\ -i \partial_z \delta \tilde{p}(\Omega, \omega, z) = & (-2\beta\omega\Omega + \gamma(2-\kappa)(N_y^0 - N_x^0) - 2\alpha)\delta \tilde{m} \\ & + \frac{\gamma\kappa}{2\pi}(n_x^0(\omega) - n_y^0(\omega) + \frac{\Omega}{2}\partial_\omega n_y^0(\omega) + \frac{\Omega}{2}\partial_\omega n_x^0(\omega)) \int (\delta \tilde{m}(\Omega, \omega_1, z) + \delta \tilde{p}(\Omega, \omega_1, z)) d\omega_1. \end{aligned} $$

Introducing the Laplace transforms of the anomalous correlator $\delta \hat{m}(\Omega, \omega, \lambda) = \int_0^z ds \delta \tilde{m}(\Omega, \omega, s) \exp(-\lambda s)$, $\delta \hat{p}(\Omega, \omega, \lambda) = \int_0^z ds \delta \tilde{p}(\Omega, \omega, s) \exp(-\lambda s)$, we obtain the dispersion relation for the instability growth-rate of the anomalous correlator Re[$\lambda(\Omega)$]:

$$ \frac{2\pi}{\gamma\kappa} = \int \frac{n_x^0(\omega) - n_y^0(\omega) - (\Omega/2)\partial_\omega(n_x^0(\omega) + n_y^0(\omega))}{i\lambda - 2\beta\omega\Omega - \gamma(2-\kappa)(N_y^0 - N_x^0) + 2\alpha} + \frac{n_x^0(\omega) - n_y^0(\omega) + (\Omega/2)\partial_\omega(n_x^0(\omega) + n_y^0(\omega))}{-i\lambda + 2\beta\omega\Omega - \gamma(2-\kappa)(N_y^0 - N_x^0) + 2\alpha} d\omega. \quad (48) $$

We have considered initial Gaussian spectra of the form: $n_\mu^0(\omega) = \frac{\sqrt{2\pi}}{\sigma} N_\mu^0 \exp(-\omega^2/(2\sigma^2))$. The analysis shows that, in general, Re[$\lambda(\Omega)$] is peaked on $\Omega = 0$, i.e., the homogeneous mode is the most unstable. Note in this respect that the growth-rate of the homogeneous mode does not depend on the initial spectrum: Regardless of the particular form of $n_\mu^0(\omega)$, the growth-rate for the homogeneous mode reads:

$$ \lambda(\Omega = 0) = \sqrt{(\gamma(2-\kappa)\Delta N^0 - 2\alpha)(\gamma(3\kappa - 2)\Delta N^0 + 2\alpha)}. \quad (49) $$

Note that this expression recovers Eq.(6) (main text) for $\kappa = 2/3$.

Fig. 6 shows the influence of the various parameters on Re[$\lambda$]. We can see that in most cases, the maximum growth, i.e. the maximum of Re[$\lambda$], is achieved for $\Omega = 0$. Increasing the spectral width, $\sigma$, reduces the range of frequencies $\Omega$ that are unstable without changing the maximum instability rate. Decreasing the initial power difference, $\Delta N^0$, reduces both the range of unstable modes $\Omega$ and the instability growth-rate. The influence of $\alpha$ is more subtle and reducing $\alpha$ can lead to a shift of the maximum instability rate from $\Omega = 0$ to $\Omega \neq 0$.
<figure><img src="image_8.png" /><figcaption>FIG. 6: Impact of the various parameters on Re[λ]. Panel (a) shows the influence of varying αLnl while keeping ΔN⁰/N = 1 and στ₀ = 2π. Panel (b) shows the influence of varying ΔN⁰/N while keeping αLnl = 1/3 and στ₀ = 2π. Panel (c) shows the influence of varying στ₀ while keeping αLnl = 1/3 and ΔN⁰/N = 1.</figcaption></figure>

## 2. Stability analysis of the normal correlators, Eqs. (45-46)

Following a procedure similar to that for the anomalous correlator, we now carry out the stability analysis for the normal correlators whose evolution is governed by the linearized Eqs. (45-46). Defining δ̂n_μ(Ω, ω, λ_n) = ∫_0^z ds ∫ dt δn_μ(t, ω, z) exp(-iΩt - λ_n s), we obtain the set of coupled equations:

$$ \delta \hat{n}_x(\Omega, \omega, \lambda_n) = \frac{i\gamma\Omega\partial_\omega n_x^0(\omega)}{2\pi(\lambda_n + 2i\beta\omega\Omega)} \int 2\delta\hat{n}_x(\Omega, \omega_1, \lambda_n) + \kappa\delta\hat{n}_y(\Omega, \omega_1, \lambda_n)d\omega_1, \quad (50) $$

$$ \delta \hat{n}_y(\Omega, \omega, \lambda_n) = \frac{i\gamma\Omega\partial_\omega n_y^0(\omega)}{2\pi(\lambda_n + 2i\beta\omega\Omega)} \int 2\delta\hat{n}_y(\Omega, \omega_1, \lambda_n) + \kappa\delta\hat{n}_x(\Omega, \omega_1, \lambda_n)d\omega_1, \quad (51) $$

The dispersion relation λ_n(Ω) for the normal correlator is obtained by solving the integral equation:

$$ (1 - 2\mathcal{I}_x(\Omega, \lambda_n))(1 - 2\mathcal{I}_y(\Omega, \lambda_n)) = \kappa^2 \mathcal{I}_x(\Omega, \lambda_n)\mathcal{I}_y(\Omega, \lambda_n), \quad (52) $$

where the integrals I_μ(Ω, λ_n) are given by

$$ \mathcal{I}_{\mu}(\Omega, \lambda_n) = \frac{i\gamma\Omega}{2\pi} \int \frac{\partial_{\omega} n_{\mu}^{0}(\omega)}{\lambda_n + 2i\beta\omega\Omega} d\omega = -\frac{\gamma\beta\Omega^2}{\pi} \int \frac{n_{\mu}^{0}(\omega)}{(\lambda_n + 2i\beta\omega\Omega)^2} d\omega, \quad \mu = x, y. \quad (53) $$

We consider initial Gaussian spectra of the form: $n_μ^0(\omega) = \frac{\sqrt{2\pi}}{\sigma} N_μ^0 \exp(-\omega^2/(2\sigma^2))$. We have solved numerically Eq. (52) to compute the dispersion relation λ_n(Ω). In the weakly nonlinear regime verifying $L_{lin} \sim (\beta\sigma^2)^{-1} \ll L_{nl} \sim (\gamma N)^{-1}$, the numerical analysis shows that the normal correlator is always stable, Re(λ_n(Ω)) < 0.

## C. Nonlinear dynamics with homogeneous statistics: Derivation of Eq. (8) (main text)

The analysis of the dispersion relation has revealed that, in general, the homogeneous mode Ω = 0 is the most unstable mode, see Fig. 6. Here, we analyze the nonlinear dynamics of Eq. (40-42) in the limit of stationary statistics:

$$ \partial_z n_x(\omega, z) = +4\gamma\kappa M^r m^i(\omega), \quad (54) $$

$$ \partial_z n_y(\omega, z) = -4\gamma\kappa M^r m^i(\omega), \quad (55) $$

$$ \partial_z m(\omega, z) = -i(\gamma(2-\kappa)(N_y - N_x) - 2\alpha)m - 2i\gamma\kappa(n_x - n_y)M^r. \quad (56) $$

Note that, $\langle \tilde{u}_y(\omega_1, z) \tilde{u}_x^*(\omega_2, z) \rangle = 2\pi m(\omega_1, z) \delta(\omega_1 - \omega_2)$ under the assumption of homogeneous statistics ($\tilde{u}_\mu(\omega, z) = \int u_\mu(t, z) \exp(-i\omega t) dt$ being the Fourier transform of the field), so that anomalous correlations between $u_x$ and $u_y$ arise only at the same frequency.

By integration over ω, we get the following system for the normal and anomalous correlators:

$$ \partial_z \Delta N = -8\gamma\kappa M^r M^i, \quad (57) $$

$$ \partial_z M^r = (\gamma(2-\kappa)\Delta N - 2\alpha) M^i, \quad (58) $$

$$ \partial_z M^i = -(\gamma(2-3\kappa)\Delta N - 2\alpha) M^r. \quad (59) $$
where $\Delta N(z) = N_y(z) - N_x(z)$, and we recall that $N_\mu(z) = \frac{1}{2\pi} \int n_\mu(\omega, z) d\omega$, $M(z) = \frac{1}{2\pi} \int m(\omega, z) d\omega$. Introducing the Stokes vector $S(z) = (\Delta N(z), 2M^r(z), 2M^i(z))^T$, the dynamics (57-59) can be recast as a rotation of this vector on the Poincaré sphere (see Eq.(8) in the main text):

$$ \partial_z S(\tau, z) = R(z) \times S(z), \qquad (60) $$

$$ R(z) = (2\alpha - \gamma(2 - \kappa)S_1(z), -2\gamma\kappa S_2(z), 0)^T. \qquad (61) $$

Because $S \cdot \partial_z S = 0$, the vector $S(z)$ evolves on the surface of the Poincaré sphere of constant radius

$$ S_0 = \left( \sum_{j=1}^{3} S_j^2 \right)^{1/2} = \sqrt{\Delta N(z)^2 + 4|M(z)|^2} = \text{const.} \qquad (62) $$

Consequently, the degree of polarization is conserved during evolution, $P(z) = S_0/N = \text{const}$.

Considering the case $\kappa = 2/3$ relevant to our experiments: The fixed point $S_1 = -S_0$ is stable. Conversely, as discussed in the main text through Eq.(8), the fixed point $S_1 = +S_0$ becomes unstable for $\alpha < \alpha_c = \frac{2}{3}\gamma S_1$, giving rise to two new stable fixed points located at $S_1 = 3\alpha/(2\gamma)$, $S_2 = 0$, $S_3 = \pm \sqrt{S_0^2 - S_1^2}$.

The Poincaré-Stokes formalism employed here differs fundamentally from the one commonly used to describe the polarization dynamics of fully coherent stationary waves [9]. First, the Stokes vector represents the correlation functions of an incoherent random wave, rather than the coherent field amplitudes. Second, in the coherent stationary case, the radius of the Poincaré sphere $S_0$, is determined solely by the total power $N$; consequently, all solutions with fixed power $N$ evolve on the surface of a single sphere. In contrast, in the present formulation the sphere radius also depends on the degree of polarization, $S_0 = PN$. As a result, $S_0$ varies with the initial conditions even when the total power is held fixed ($N = \text{const}$): The full set of solutions is not described by a single sphere, but instead by an ensemble of spheres with different radii.

### III. EXPERIMENTAL SET-UP

The experimental set-up is presented in Fig. 7.

<figure><img src="image_9.png" /><figcaption>FIG. 7: **Experimental set-up.** Laser source (Nd:YAG), halfwave plate (λ/2), polarizer (P), microscope objectives (MO), multimode fiber (MMF), spectral filter (SF), polarization beam splitter (PBS), retroreflector (RR), optical spectrum analyzer (OSA) polarimeter (PM).</figcaption></figure>

The source is a Nd:YAG laser delivering subnanosecond pulses (400ps) at $\lambda = 532$nm. Using a microscope objective (×20), the beam is first injected into a standard GRIN multimode fiber [8]. During the nonlinear propagation, a temporally incoherent pulse is generated via the Raman beam clean-up effect which allows to produce Stokes beam with good spatial beam quality [10, 11]. At the output of the multimode fiber (MMF), we spectrally filter the second Raman Stokes, centered around $\lambda = 558$ nm with FWHM of 2 nm and characterized by a nearly Gaussian beam shape. The beam is then passed through a half-wave plate to control the power on each polarization axis before a polarization beam splitter (PBS). The half-wave plate is oriented such that there is an equal amount of power on the two orthogonal polarization axes. Each polarization propagates for a different time $\tau_0$ and $\tau_1$ using free-space delay lines before being recombined through the PBS. We introduce a time delay, $\Delta\tau = \tau_1 - \tau_0 \approx 15$ ps which is larger than the coherence time of the pulse (< 1 ps), hence removing the mutual coherence between the two polarization axes, i.e. $|M| \approx 0$. The beam is then passed through another half-wave plate in order to align the two polarizations with
the principal axes of the weakly birefringent fiber. The weak birefringence is obtained by carefully winding a spun fiber onto a spool with a diameter of 14.5 cm. Spun fibers are made by rapidly spinning a preform during the drawing process. This results in a fast rotation of the slow and fast axis of the fiber, leading the averaged net birefringence to be zero. The mechanical torsion from the winding introduces a weak, controlled linear birefringence $dn \approx 10^{-7}$. The axes of the weakly birefringent fiber were identified by measuring the polarization modulation instability of the coherent Nd:YAG beam directly injected [5]. Microscope objectives (×40) were used to inject and collect light in and out of the fiber. At the output of the fiber, we measure the optical spectrum and the power on each polarization axis. To characterize optical spectra, we used the OSA YOKOGAWA-AQ6374. Furthermore, the temporal profile of the pulse was monitored before the injection into the weakly birefringent fiber using a high-speed oscilloscope combined with 12-GHz photodiode. To measure the complex anomalous correlator, we used a homemade polarimeter made of a combination of quarter-wave, half-wave plates, polarizers, and powermeter as commonly applied in standard polarization measurement.

## IV. NUMERICAL SIMULATIONS OF THE EXPERIMENTAL REGIME

To interpret the experimental observations, we have numerically integrated the coupled NLS Eqs. (10-11), where $\beta = \beta_2/2$, with $\beta_2$ the group-velocity dispersion coefficient [9]. The coherent coupling parameter $\alpha = \Delta\beta/2$ is related to the fiber birefringence $\Delta\beta = (2\pi/\lambda_0)dn$, with $dn = |n_x - n_y|$ the refractive index difference among the axes and $\lambda_0$ the pulse central wavelength. Furthermore, in the numerical simulations of the experiments we have included the minor correction provided by the group-velocity difference of the waves $u_x$ and $u_y$, which was found to have a negligible impact. The fiber used in the experiment was the same as in Ref. [5, 6], with parameters $dn \approx 1 \times 10^{-7}$, effective area $A_{\text{eff}} \approx 15 \ \mu\text{m}^2$ and $\beta_2 \approx 60 \ \text{fs}^2/\text{mm}$ at $\lambda_0 \approx 558$ nm, nonlinear refractive index $n_2 = 3.2 \times 10^{-20} \text{m}^2/\text{W}$, and corresponding nonlinear parameter $\gamma = 2\pi n_2/(\lambda_0 A_{\text{eff}}) \approx 0.024 \ \text{W}^{-1}\text{m}^{-1}$.

The spectrum and temporal profile of the initial pulses were recorded experimentally before injection into the fiber and served as initial conditions in the numerical simulations (spectrum with a FWHM of $\approx 2$-nm and 100-ps long pulses). A random phase was added on each spectral component to generate an incoherent temporal wave.

### A. Thermalization in the experimental regime

In the main text, we have reported the repolarization effect. Due to experimental limitations, we have been able to investigate the early stage of the thermalization process. As discussed in [7], the optical field approaches at each z a local quasi-equilibrium within a limited spectral window. As z increases, this spectral window grows and the spectral distribution follows, locally, an equilibrium state, characterized by the power-law decay $\sim \omega^{-2}$ in the spectral tail. This is illustrated in Fig. 8, which shows the evolution of the spectrum as z increases. We can see that the spectral window of local equilibrium increases as the field propagates over large propagation lengths, exhibiting the thermalization process in the experimental configuration. This confirms that the experimental results reported through Fig. 3 (main text) provide a signature of the thermalization process.

<figure><img src="image_10.png" /><figcaption>FIG. 8: Local quasi-equilibrium in the experimental regime. Numerical simulation of the NLSE in the experimental configuration for a peak power of $\approx 18$ W. The spectral window in which the field is in local quasi-equilibrium increases as z increases.</figcaption></figure>
[1] V.E. Zakharov, V.S. L'vov, G. Falkovich, *Kolmogorov Spectra of Turbulence I* (Springer, Berlin, 1992).

[2] V. Zakharov, F. Dias, A. Pushkarev, Phys. Rep. **398**, 1 (2004).

[3] S. Nazarenko, *Wave Turbulence* (Springer, Lectures Notes in Physics, 2011).

[4] A. Picozzi, Opt. Express **16**, 17171 (2008).

[5] G. Millot, E. Seve, S. Wabnitz, M. Haelterman, J. Opt. Soc. Am. B **15**, 1266 (1998).

[6] G. Millot, E. Seve, S. Wabnitz, M. Haelterman, Opt. Lett **23**, 511 (1998).

[7] L. Zanaglia, J. Garnier, I. Carusotto, V. Doya, C. Michel,

A. Picozzi, arXiv:2506.23536

[8] K. Krupa, A. Tonello, A. Barthélémy, V. Couderc, B.M. Shalaby, A. Bendahmane, G. Millot, S. Wabnitz, Phys. Rev. Lett. **116**, 183901 (2016).

[9] G. Agrawal, *Nonlinear Fiber Optics* (sixth Ed., Academic Press, New York, 2019).

[10] N. B. Terry, T. G. Alley, and T. H. Russell, Optics express, 15(26), 17509-17519 (2007).

[11] H. Pourbeyram, G. P. Agrawal, A. Mafi, Appl. Phys. Lett. 102 (20):201107 (2013).