Dynamical Thermalization, Rayleigh-Jeans Condensate, Vortexes and Wave Collapse in Quantum Chaos Fibers and Fluid of Light

1. Introduction

In recent years a number of interesting physical phenomena have been established for laser beam propagation in multimode optical fibers (see e.g., overview [1]). Among them there is a phenomenon of self-cleaning of the laser beam when an injected laser power at high modes is transferred to low-energy modes. It was realized that this effect appears due to thermalization among linear fiber modes induced by nonlinear four-wave interactions between modes [2,3,4,5,6,7]. Thus, even thermalization with negative temperatures has recently been observed in such optical fibers [7]. It is argued that this thermalization emerges as a result of purely Hamiltonian dynamics without any external thermal bath.

The phenomenon of dynamical thermalization attracts a deep interest of the scientific community, starting from the famous Boltzmann-Loschmidt dispute on thermalization emerging from the reversible dynamical equations of particle motion [8,9,10] (see also [11]). In modern times the interest in the problem of dynamical thermalization was pushed forward by numerical experiments of Fermi, Pasta, Ulam (FPU) in 1955 [12] with the conclusion that “The results show very little, if any, tendency toward equipartition of energy between the degrees of freedom.” [12]. The absence of thermalization in the FPU problem of linear oscillators with nonlinear couplings was found to be linked to its proximity to exactly integrable models with solitons [13,14,15] and Toda lattice [16,17]. Also, at weak nonlinearity there is no overlapping of main system resonances, the system is below the chaos border, and the dynamics remains quasi-integrable [18,19,20]. Such a regime corresponds to a spirit of Kolmogorov-Arnold-Moser (KAM) integrability which takes place at weak perturbation of integrable system while chaos and thermalization may appear only at relatively strong nonlinear perturbation (see details in [21,22,23,24]). At present the FPU problem still remains an active area of research [25].

In fact an example of FPU problem shows that for investigations of dynamical thermalization in classical many-body systems it is important to choose a model corresponding to a generic situation of oscillators with nonlinear interactions. Such a generic model was proposed in [26] with unperturbed Hamiltonian of linear oscillators described by a random matrix with additional nonlinear perturbations. Thus the linear oscillators are described by the Random Matrix Theory (RMT) invented by Wigner for a description of complex quantum systems like nuclei, atoms and molecules [27]. At present RMT finds a variety of applications in quantum systems, mesoscopic physics [28,29] and systems of quantum chaos [30,31]. In [26] it was shown that a nonlinear perturbation of RMT leads to dynamical chaos and the Rayleigh-Jeans thermal distribution over linear oscillator modes if nonlinearity strength exceeds a certain chaos border. In presence of additional dissipation and energy pumping this nonlinear random matrix (NLIRM) model describes also the spectra of Kolmogorov-Zakharov turbulence [32]. A brief discussion of a model similar to NLIRM was also done in [33].

In practice it may be not so easy to reach an experimental realization of such NLIRM model but in [26] it was proposed that this system can be approximately presented by a multimode optical fiber with a fiber cross-section being a chaotic billiard, e.g., D-shape billiard (circle with a cut. see Figure 1). Indeed, it is known that classically chaotic billiards have many properties being the same as for RMT [30,31]. It was shown that the quantum chaos exists for variety of cuts of a circle billiard [34,35,36] even if for certain cut positions there can be present significant integrable islands. The advantage of D-shape billiard is that it can be relatively simply realized in multimode fibers. Indeed, observations of signatures of quantum chaos in D-shape billiard had been reported in [35,36,37]. However, dynamical thermalization induced by nonlinearity was not studied in these works. We argue that such quantum chaos fibers can be used for investigations of dynamical thermalization for a generic case of oscillators with nonlinear interactions. In principle similar effects can be studied for lasing dynamics in semiconductor chaotic microcavities realized in [38] but in this work we restrict our studies to the case of D-shape multimode fibers.

It is natural to ask what is the situation with dynamical thermalization observed with multimode fibers in [3,4,5,6,7]. In fact in these circular shape fibers the spectrum of linear modes corresponds to a spectrum of two oscillators with equal frequencies that leads to a significant spectrum degeneracy. In such a case the KAM theory for nonlinear perturbation of such a system is not applicable [21,22,23,24]. Thus in [39] it was shown that for three oscillators with equal frequencies, the measure of chaotic component remains on a level of 50% for arbitrary weak nonlinear perturbation (its strength only gives a reduction of the positive Lyapunov exponent but not the measure of chaos). Similar effects exist for other nonlinear systems [40] and the FPU $\alpha$-model [41]. Thus for fibers, which have spectrum of two oscillators with equal frequencies, as those in [3,4,5,6,7], we expect that arbitrary small four-wave nonlinear interaction, that couples these modes, should lead to a chaotic dynamics and eventual thermalization between degenerate modes. However, since the Lyapunov exponent decreases with the nonlinearity strength (as in [39,40,41]) one would need to have longer and longer fiber to reach such chaotization. It is possible that certain non homogeneous perturbations along a fiber act as a time-dependent noise and facilitate development of chaos. Due to these reasons we argue that the experiments [3,4,5,6,7] are done in a specific regime of degenerate mode spectrum when the KAM theory is not valid and thus it is of fundamental interest to study the generic case of D-shape fiber when the KAM theory works and dynamical thermalization appears only above a certain chaos border. We should note that the mathematical KAM theory results for the NSE in billiards are very difficult to obtain and the question about asymptotic time behavior of nonlinear spreading on high energy modes remains open for chaotic billiards and moderate nonlinearity (see e.g. [42]).

With this aim we present here the numerical and analytical study of dynamics of nonlinear waves in D-shape billiard described by the Nonlinear Schrödinger Equation (NSE), or the Gross-Pitaevskii equation (GPE) [43], which depicts the laser beam propagation in multimode fibers [1]. We note that the onset of dynamical thermalization for the NSE with defocusing nonlinearity $\beta >0$ in a chaotic Bunimovich billiard was studied in [44]. The results reported there show emergence of dynamical thermalization above a chaos border when nonlinear perturbation exceeds a typical energy level spacing in such a billiard. However, in [44] only relatively short times were reached in numerical simulations and the Rayleigh-Jeans distribution was not reached. Here we significantly increased the time evolution scale with emergence of the Rayleigh-Jeans thermalization in a D-shape fiber (billiard) with defocusing nonlinearity.

Another important feature is that usually for fibers a nonlinearity is focusing $\beta <0$ and thus, according to the so called Vlasov-Petrishchev-Talanov theorem [45,46], a wave collapse can take place for strong nonlinearity at finite times. Important implications of such a collapse for Langmuir waves was discussed in [47]. However, in [45,46,47] a collapse was discussed in a case of infinite system size. For a D-shape fiber chaos and thermalization in a finite size billiard may lead to an unusual interplay between collapse and dynamical thermalization so that there is a fundamental interest to study the properties of NSE in a D-shape focusing fiber. Here we present results of our studies of such a system considering both cases of negative (focusing) and positive (defocusing) signs of nonlinearity. We show that that the dynamical thermalization with the Rayleigh-Jeans distribution takes place at moderate strength of nonlinearity being above the chaos border. Below this border the dynamics is quasi-integraable corresponding to the KAM theory.

It is interesting to note that the GPE describes Bose-Einstein condensate (BEC) of cold boson atoms in optical trap [43]. Thus there is a temptation to expect [44] that the dynamical thermalization for GPE or NSE in a chaotic billiard would lead to the thermal Bose-Einstein (BE) distribution over eigenenergies $E_m$ of linear system [48]:

$${\rho }_m={\displaystyle \frac{1}{exp[(E_m-\mu )/T]-1}}\left(\mathrm{BE}\right).$$

(1)

However, NSE describes classical nonlinear fields without second quantization and thus in this case we should have the classical Rayleigh-Jeans (RJ) thermal energy equipartition with chemical potential $\mu$ due to an additional conservation integral of norm or number of particles, see e.g. [48,49]:

$${\rho }_m={\displaystyle \frac{T}{E_m-\mu }}\left(\mathrm{RJ}\right).$$

(2)

In the distributions (1), (2) T is a system temperature, $\mu$ is a chemical potential and the total norm of probabilities ${\rho }_m$ on N energy levels $E_m$ is preserved and normalized to unity (${{\sum }^N}_{m=1}{\rho }_m=1$; $1\le m\le N$). In optical fibers such classical thermal distribution (2) is usually called as Rayleigh-Jeans (RJ) distribution (see e.g. [3]). Indeed, in [26] for the NLIRM and related models it was shown that above the chaos border the dynamical chaos at moderate nonlinearity leads to the RJ thermal distribution and not to the BE one. For dynamical systems like NSE we assume that the nonlinear term is relatively weak or moderate and it leads to dynamical chaos for nonlinearity above the chaos border producing thermal distribution (2) over linear eigenenergies (or mode eigenfrequencies) that are not significantly affected by weak nonlinearity.

In fact the RJ distribution is a limiting case of the BE one in the limit of high temperature significantly exceeding typical energies of the linear system ($T\gg E_m$). Using this RJ description of BE distribution at high temperature Fröhlich showed in 1968 that for the RJ distribution there is a significant condensation of probability at the lowest system energy level that may be close to hundred percent of total probability norm or total number of bosons [50,51]. The model proposed by Fröhlich assumes a presence of external pumping of the system in presence of dissipative process leading to a certain steady-state. The total number of phonons in the system depends on a pumping strength leading to condensation in the ground state above a certain critical strength. Thus this model does not directly correspond to the RJ thermal distribution studied here, but there are still certain similarities due to the fact that the steady-state in [50,51] is described by the RJ distribution with a certain prefactor depending on pumping and dissipation.

On the basis of this result Fröhlich made a conjecture that such a condensate, appearing at high temperatures, concentrates energy in a single mode and thus it can play an important role in coherent properties of cell membranes and biomolecules. This conjecture attracted a significant interest of physicists and biologists with a large number of theoretical and experimental studies of this phenomenon called Fröhlich condensate, see e.g. [52,53,54,55,56,57].

In this work we show that in the RJ thermal distribution (2) there is also condensation at the ground state, which we call Rayleigh-Jeans (RJ) condensate. This RJ condensate captures almost all norm, or number of particles, it emerges at low temperatures ($T\ll E_m$) and high number of levels in the system. Since RJ condensate appears only at low temperatures it describes the thermalization of classical fields since for quantum systems the thermalization is described by the BE distribution (1) which at low temperatures is very different from RJ one (2). Thus the RJ condensate exists only at low temperatures ($T\ll E_m$) at norm conservation and absence of external pumping. Hence it is very different from the Fröhlich condensate which exists at high temperatures ($T\gg E_m$) in presence of external pumping.

We argue that the self-cleaning effect [1,2] is an example of the RJ condensate appearing due to chaos induced by nonlinear interaction of waves in multimode fibers. The emergence of such condensation was also reported for numerical simulations of NSE in two and three dimensional finite systems with quadratic spectrum of waves [58,59,60,61]. In these works the emergence of RJ distribution was attributed to the Kolmogorov-Zakharov spectra of turbulence [49]. However, as argued above, we attribute the appearance of RJ thermal distribution in these numerical simulations to the degeneracy of linear eigenmode frequencies due to which the KAM theory cannot be applied to these systems and thus chaos appears at very weak nonlinearity strength. For the case of focusing fibers we discuss an interesting interplay between the wave collapse and chaos. We also consider the problem of ultraviolet catastrophe at high energy modes in a chaotic regime.

Finally, we note that a laser beam propagation in a nonlinear Kerr media with atomic vapor [62] is also described by the NSE which we study here. The important advantage of this physical system is that the sign of NSE nonlinearity can be easily changed by a frequency detuning of a laser beam. In [62] it was experimentally and numerically demonstrated a formation of vortexes of such light fluid described by a dissipative NSE within a circular integrable billiard cross-section. Here we show that in the case of chaotic D-shape cross-section there are also long living vortexes described by the unitary NSE evolution. The NSE dynamics of laser beam propagation in fiber or other nonlinear media is also considered as a fluid of light [61,62]. In fact the NSE and GPE are closely related to the Ginzburg-Landau equation used for the description of superconductivity, superfluidity and quantum turbulence as it is described in the reviews [63,64,65,66].

The article is composed as follows: Section 2 presents the model and description of quantum chaos properties in absence of nonlinearity, numerical methods of integration of NSE are given in Section 3, results for simple models with RJ condensate are displayed in Section 4, the results for dynamical thermalization for defocusing fiber NSE are described in Section 5, properties of the wave collapse for focusing NSE are presented in Section 6, the dynamical thermalization for focusing fiber NSE is described in Section 7, properties of vortexes in defocusing fiber NSE are displayed in Section 8, NSE long time dynamics is discussed in Section 9, experimental parameters for a quantum chaos fiber are displayed in Section 10, discussion and conclusion are given in Section 11. Appendix presents some additional results.

2. Model Description

The D-shape fiber model is described by the NSE with the Dirichlet boundary conditions:

$$i\hslash {\displaystyle \frac{\partial \psi (\overrightarrow{r},t)}{\partial t}}=-{\displaystyle \frac{{\hslash }^2}{2m}}\Delta \psi (\overrightarrow{r},t)+\beta {\left|\psi (\overrightarrow{r},t)\right|}^2\psi (\overrightarrow{r},t)$$

(3)

where we consider $\hslash =1$, $m=0.5$, $\overrightarrow{r}$ is a vector in $(x,y)$ plane. The form of D-shape fiber billiard is shown in Figure 1. The time variable t corresponds to z direction of a beam propagation along fiber (see e.g. Eq.(1) in [1]), $\psi$ represents a field inside fiber and $\beta$ is a field strength coefficient of nonlinearity. This equation has two integrals of motion being the norm $\int |\psi (\overrightarrow{r},t){|}^{2}dxdy=1$ and energy $E=<\psi \left(t\right)|\widehat{H_0}{|\psi \left(t\right)>+\beta |\psi \left(t\right)|}^4/2$ where ${\widehat{H}}_0$ is the operator of linear Hamiltonian at $\beta =0$.

Here we study the D-shape billiard with $w=1+1/\sqrt{2}$ shown in Fig.Figure 1. In this case the Poincaré section [24] of classical ray dynamics is shown in the left panel of Figure 2. The phase space is fully chaotic even if very small islands of regular motion are not excluded. For other values of w it is possible to have large islands with regular dynamics as shown in the right panel of Figure 2. Our results show that in the range of $1.7\le w\le 1.95$ the dynamics is fully chaotic similar to the case of $w=1+1/\sqrt{2}$.

Millions of eigenvalues $E_m$ and eigenvectors of linear Hamiltonian $\widehat{H_0}$ can be found numerically with the efficient methods developed in the field of quantum chaos (see e.g. [67]). At $\beta =0$ we consider the eigenstates ${\psi }_m$ with even-odd parity in y ($\psi (x,y,t)=\pm \psi (x,-y,t)$). Usually multimode optical fiber may have up to about thousand of modes that determined our consideration with up to 1000 states.

Examples of several eigenstates ${\psi }_m(x,y)$ are shown in Figure 3 with eigeneneries $E_m$ of first 100 states. A chaotic structure of wave function eigenstates at high energies is evident. According to the Weyl theorem the average energy level spacing between adjacent levels $\Delta$ is independent of energy E being $\Delta \approx 4\pi /A\approx 4.4$ where A is the billiard area [31] ($A=3\pi /4+1/2\approx 2.856$). The level spacing statistical distribution $P\left(s\right)$ for the lowest 2000 energies $E_m$ for even and odd parity eigenstates is shown in Figure 4. It is in a good agreement with the Wigner result from RMT and other results for quantum chaos billiards [30,31].

3. Numerical Integration of NSE

For numerical integration of NSE (3) we use two different methods. The first one uses the same strategy as in [44]. The time evolution of (3) is integrated by small time steps of Trotter decomposition of linear and nonlinear terms with a step size going down to $\Delta t=4\times {10}^{-5}$. There are up to $N_l=1085$ linear eigenmodes ${\varphi }_m$ for linear part of time propagation; then the wave function is transferred in coordinate space $(x,y)$ with $N_p\approx 3N_l$ points inside the billiard. The change of basis from coordinates to energies (and vice versa) is given by a unitary matrix in double precision. At any step the wave function is expanded in the basis of linear modes ${\varphi }_m$ so that $\psi (x,y,t)={\sum }_m{C}_m\left(t\right){\varphi }_m(x,y)$. A special aliasing procedure [68] is used with an efficient suppression of nonlinear numerical instability at high modes. This integration scheme exactly conserves the probability norm providing the total energy conservation with an accuracy better than $2\%$ at large values $\beta =10$ and better than $1\%$ at lower $\beta$ values. But this approach allows to obtain numerical results only on relatively short times with $t\le 40$.

Due to the numerical restrictions of the first approach we perform all presented numerical simulations with the second approach using the public codes of FREEFEM package [69]. The number of net nodes of FREEFEM used in simulations is usually up to $N_F=136241$ with the integration time step $\Delta t=dt={10}^{-4}$ for the case of vortex dynamics at such high $\beta$ values as $\beta =45$. These parameters depend on nonlinearity and initial state usually composed from linear eigenmodes ${\psi }_m$. For typical cases with $\beta =10$ we use $\Delta =5\times {10}^{-4}$ and $N_F=34501$ for maximal time studied $t=1000$. For such typical cases the norm is conserved with accuracy of approximately 3% and total energy conservation is at 6% at $t=1000$ ($\beta =10$, initial state at $m=11$). Another typical example of computation is $\beta =-10$, with initial $m=11$, it has a CPU run time of 2 weeks on a workstation with 100 processors with $N_F=136241,\Delta t={10}^{-4}$ and reached time $t=800$ with accuracy of norm and energy relative conservation of 0.05% and 1% respectively. In certain cases with collapse we used even higher $N_F$ and smaller $\Delta t$ values (see below).

We checked that the obtained results are not sensitive to a variation of $N_F$ and $\Delta t$. At short times $t\le 40$ both methods give close results but such long times at $t=1000$ can be reached only with the FREEFEM codes with which all presented NSE results are obtained.

4. Model Systems with RJ Condensate

4.1. Simple Estimates

To understand the properties of RJ condensate we consider a simple model with N equidistant energy levels $E_m=m\Delta$ with level index $m=1,2,\dots ,N$ and $\Delta$ neing a level spacing which we suppose to be unity. The probability norm or number of particles (bosons) is conserved so that ${{\sum }^N}_{m=1}{\rho }_m=1$ where ${\rho }_m$ is a probability at level m. The system ground state has energy $E_g=E_1=1$. We assume that this system of linear modes m is coupled to an external thermostat or that an additional weak nonlinear interaction of modes (e.g. as in (3)) leads to dynamical thermalization.

Since we have classical modes or fields and two integrals of motion, being total energy ${ϵ}_0$, initially injected in the system, and the total norm being unity, the thermal distribution is the RJ one given by Eq.(2) with temperature T and chemical potential $\mu$ as it was shown in [26] for the NLIRM system. These two parameters are determined by two integrals of norm and total energy:

$$\sum _m{\rho }_m=T\sum 1/(E_m-\mu )=1;\sum _m{E}_m{\rho }_m=T\sum E_m/(E_m-\mu )={ϵ}_0;{ϵ}_0=NT+\mu ,$$

(4)

where the last equality directly follows from (2). If we have a case of a system coupled to a thermostat with temperature T then the chemical potential is determined by the first equality of norm being unity and then the total energy ${ϵ}_0$ is given by the last equality. In the case of dynamical thermalization we have initial system energy ${ϵ}_0$ being the total energy and since energy is the integral of motion then the first and second equalities in (4) determine the system temperature T and chemical potential $\mu$.

A case when energy levels grow approximately linearly with the level number is rather generic and standard. Due to that we call this simple model with equidistant levels as the RJ Standard (RJS) model.

As for the case of BEC [43,48] at low temperature we have chemical potential approaching to ground state energy $E_g=E_1$ from below ($\mu \le E_1$) leading to a concentration of probability at the ground state. Since $\mu$ is very close to $E_1$ we obtain from (4) the approximate condensation temperature $T_c$ at which the fraction of condensate at the round state $E_g$ is very close to unity:

$$T_c\approx ({ϵ}_0-E_g)/N.$$

(5)

Then the fraction at the excited states at $2\le m\le N$ is $W_{ex}$ and the condensate fraction is $W_c$:

$$W_{ex}\approx T{\sum ^N}_{m=2}1/(E_m-E_g)\approx (T/\Delta )lnN\approx ({ϵ}_{0}lnN)/N\Delta ;W_c\approx 1-W_{ex};\Delta =1.$$

(6)

Thus, in the RJS model for a fixed initial or total energy ${ϵ}_0$ the fraction of non condensate drops to zero as $W_{ex}\propto (lnN)/N$ with growing number of system levels N. When the condensate fraction is close to unity $W_c\to 1$ and $W_{ex}\ll 1$ we see that this situation takes place at very low temperatures $T\approx ({ϵ}_0-E_g)/N\ll E_m$ where $E_m\sim N\Delta /2$ is a typical system energy of excited levels. This shows that the RJ condensate is very different from the Fröhlich condensate [50,51] which exists at high temperatures $T\sim E_m\sim N\Delta /2$. We have $W_c\approx W_{ex}\approx 1/2$ at ${ϵ}_0\approx N\Delta /(2lnN)$. If an initially injected energy is ${ϵ}_0\sim B\Delta$ with a constant $B\sim 1$ (e.g. a few units of $\Delta$) then with increasing N almost all probability is located in the ground state while probability at $m\ge 2$ drops as $1/N$. Thus for such an initial energy ${ϵ}_0$ there is no ultraviolet catastrophe since the total probability on levels $m\ge 2$ drops as $(lnN)/N$.

4.2. Numerical Results

Here we present the numerical results that confirm the estimates of previous subsection. We assume that an initial energy ${ϵ}_0$ is injected in the above RJS model with eigenenergies $E_m=m\Delta$ of N modes. The two equalities for norm and total energy in (4) determine temperature T and chemical potential $\mu$ of the RJ thermal distribution (2). The solution always exists, at fixed ${ϵ}_0$ and $N\gg 1$ we have $\mu \to E_g$.

In Figure 5 we depict the curves of fixed condensate probability ${\rho }_1$ at the ground state $m=1$ at different values of initial energy ${ϵ}_0$ and total number of states N. The results show that very high values of ${\rho }_1=0.9;0.99$ can be obtained for high values of N and low values of ${ϵ}_0$ that corresponds to very low temperature $T\approx {ϵ}_0/N\ll {ϵ}_0$.

Indeed, the results of Figure 6 in the left panel show that high probability of condensate ${\rho }_1$ can be obtained at high N only at very low temperature $T<\Delta$ in agreement with above simple estimates (5),(6). The right panel shows the same ${\rho }_1$ but in the plane $(N,{ϵ}_0)$; it also shows that the condensate exists only at low temperature $T\approx {ϵ}_0/N<\Delta$.

The probability thermal distributions ${\rho }_m$ are shown in Figure 7 at various total number of modes N and fixed total energy ${ϵ}_0=2.5$. In agreement with estimates (5),(6), the probability at $m\ge 2$ drops with increase of N while condensate probability ${\rho }_1$ approaches unity. The dependence of total non condensate probability $W_{ex}=1-{\rho }_1$ at $m\ge 2$ is shown in Figure 8 for several values of ${ϵ}_0$. It decreases as $W_{ex}\propto (lnN)/N$ is agreement with estimate (6) and analytical result of next subsection.

We also should point that for the RJS model we consider the RJ condensate at low positive temperatures corresponding to the case when initial energy ${ϵ}_0<N\Delta /2$. However, there is a symmetric situation for negative temperature being close to zero from below $T<0$ when due to symmetry (see [26]) the RJ condensate appears on the top level $m=N$ (negative temperatures appear when ${ϵ}_0>N\Delta /2$). This RJ condensate at negative temperature appears in systems with bounded energy spectrum $E_m$ of linear eigenmodes. For our NSE billiard (3) the spectrum $E_m$ is unbounded and due to that we do not discuss here RJ condensation at negative temperatures (we remind that RJ distribution at negative temperatures are realized in fiber experiments [7]).

We note that a high fraction for RJ condensate, close to hundred’s percent, has been obtained in numerical studies of nonlinear oscillator models NLIRM and NLIRM plus extra diagonal in [26] (see there Fig.3 top left panel, Fig.S4 left panels and Fig.S14 top left, right panels). However, the dependence of this condensate fraction on initial energy ${ϵ}_0$ and total number of oscillators N was not analyzed there.

Above we considered the case of RJS model with mode energies $E_m=m\Delta$. But the same approach of RJ thermal distribution (2) can be used for other systems with other spectrum. Thus in Figure 9, top panels, we show the color plot of condensate probability in the planes $(N,T)$ and $(N,{ϵ}_0)$ for the spectrum $E_m$ of linear modes in the D-shape fiber of Figure 1 (a part of spectrum is shown in Figure 3). Here we use the modes of odd and even parity since they become both excited due exponential instability induced by nonlinearity in NSE. In the bottom panels of Figure 9 we show the condensate probability for the NSE in the Sinai oscillator studied in [70].

For comparison between BE and RJ distributions we show in Appendix Figure A1 the BE condensate probability ${\rho }_1$ in the RJS model of equidistant levels at various system parameters, this Figure A1 is analogous to the RJ case shown in Figure 6. For the BEC case ${\rho }_1$ is practically independent of number of levels due to an exponential decay of probability at high levels with $E_m>T$ that allows to have relatively high BEC probability ${\rho }_1\approx 0.5$ even at $T>\Delta$ in contrast to the RJ case in Figure 6. We note that relations between BE and RJ condensates are discussed in the recent work [61].

In the works [26,44,70] the dependence of entropy on system energy $S\left(E\right)$ was used to obtain a distinction between BE (1) and RJ (2) distributions. Thus in Figure 10 we present these two theoretical dependencies for the cases of RJS model of equidistant levels (left panels) and the system with energies of D-shape billiard. The theory shows that there is a significant difference between two dependencies that grows with the increase of total number N of levels in the system that boosts the condensate fraction (see Figure 8). At low energies E the RJ entropy is significantly smaller compared to BE case due to RJ condensation. The results obtained from the NSE dynamics (3) for the dependence $S\left(E\right)$ are discussed in the next Section 5.

4.3. Analytical Results

To describe analytically the population of eigenstates in the thermodynamic limit of a large number of modes N, we assume a microcanonical ensemble for a simplified model with evenly spaced eigen-energies ${ϵ}_n=n-1$. In this case the total energy is given by:

$$\begin{array}{c}\hfill ϵ=\sum _{n=1}^N(n-1){\rho }_n\end{array}$$

(7)

while the total population normalized to ${\sum }_{n=1}^N{\rho }_n=p$.

The probability distribution of the populations ${\rho }_n$ can be obtained from the partition functions:

$$\begin{array}{cc}\hfill s_N(ϵ,p)& =N!{\int }_{{\rho }_i>0}d{\rho }_1\dots d{\rho }_N\delta (ϵ-\sum _{n=1}^N(n-1){\rho }_n)\delta (p-\sum _{n=1}^N{\rho }_n)\hfill \end{array}$$

(8)

with the short notation $s_N\left(ϵ\right)=s_N(ϵ,p=1)$.

For example the probability distribution of the ground state population $p_1^{\left(N\right)}\left({\rho }_1\right)$ is given by:

$$\begin{array}{c}\hfill p_1^{\left(N\right)}\left({\rho }_1\right)=N{(1-{\rho }_1)}^{N-3}{\displaystyle \frac{s_{N-1}((ϵ+{\rho }_1-1)/(1-{\rho }_1))}{s_N\left(ϵ\right)}}\end{array}$$

(9)

The following analysis has close similarities with the approach from [71]. The Laplace transform of $s_N\left(ϵ\right)$ can be easily computed:

$$\begin{array}{cc}\hfill {\widehat{s}}_N(\beta ,\mu )& ={\int }_0^{\infty }dϵe^{-\beta ϵ}{\int }_0^{\infty }dp{e}^{-\mu p}{s}_N(ϵ,p)\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & =N!\prod _{n=0}^{N-1}{\displaystyle \frac{1}{\beta n+\mu }}\hfill \end{array}$$

(11)

we now need to inverse the Laplace transform to recover $s_N(ϵ,p)$.

The inverse Laplace transform over $\mu$, is known:

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mu }^{-1}\left(N!\prod _{n=0}^{N-1}{\displaystyle \frac{1}{\beta n+\mu }}\right)=N{\left({\displaystyle \frac{1-e^{-p\beta }}{\beta }}\right)}^{N-1}\end{array}$$

(12)

For a fixed N we can find the exact inverse Laplace transform over $\beta$ expanding the power in Eq. (12), this gives the explicit formula:

$$\begin{array}{c}\hfill s_N\left(ϵ\right)={\displaystyle \frac{N}{(N-2)!}}\sum _{n=0}^{N-1}{C}_{N-1}^n{(-1)}^n{(ϵ-n)}^{N-2}\theta (ϵ-n)\end{array}$$

(13)

but which is not convenient to study the $N\gg 1$ limit.

Instead, in this case we perform a saddle point approximation on the inverse Laplace transform integral:

$$\begin{array}{c}\hfill s_N(ϵ,p)={\displaystyle \frac{N}{2\pi i}}\int d\beta exp\left(ϵ\beta +(N-1)log{\displaystyle \frac{1-e^{-p\beta }}{\beta }}\right)\end{array}$$

(14)

We now fix $p=1$, in this case the saddle point approximation gives the following equation for $\beta$:

$$\begin{array}{c}\hfill {\displaystyle \frac{N-1}{e^{\beta }-1}}-{\displaystyle \frac{N-1}{\beta }}+ϵ=0\end{array}$$

(15)

Using the asymptotic expansion near the saddle point in the $ϵ\ll n$ limit, we can find an analytic expression for $p_1\left(x_1\right)$:

$$\begin{array}{ccc}\hfill p_1\left({\rho }_1\right)& ={\displaystyle \frac{N}{ϵ}}f\left(y\right),y={\displaystyle \frac{N}{ϵ}}\left(1-{\displaystyle \frac{ϵlogN}{N}}-{\rho }_1\right),f\left(y\right)\hfill & \hfill =e^{-y-e^{-y}}\end{array}$$

(16)

where $f(.)$ is the Gumbel distribution. The Gumbel distribution often appears in the description of rare events. It seems here it occurs because we are looking at the intersection of two microcanonical ensembles at fixed energy and at fixed population which is a rare event. This asymptotic agrees very well with both Monte Carlo simulations and the exact Eq. (13) (see Figure 11).

The mean of the Gumbel distribution is given by the Euler-Mascheroni constant $\gamma =0.57721$ This gives the following asymptotic expression for the population of the mean ground state population $\langle {\rho }_1\rangle$:

$$\begin{array}{c}\hfill \langle {\rho }_1\rangle =1-ϵ{\displaystyle \frac{logN+\gamma }{N}}\end{array}$$

(17)

which coincides from the behavior expected from the Rayleigh-Jeans distribution. Thus the condensate fraction $\langle {\rho }_1\rangle$ tends to 100% in the $N\gg 1$ limit. In the context of the time evolution of non-linear Schrödinger equation, this would correspond to a slow depletion eigen-modes at intermediate energies, in favor of a small and continuously decreasing population of modes of higher energy.

It is interesting to know what can limit the condensation in the large N limit in this statistical framework. Including a non-inform eigen-energy energy distribution as well as a weak nonlinearity in the expression of the total energy $ϵ$ as function of the eigen-level populations ${\rho }_n$ does not seem to stop condensation as long as there is a gap between the ground state and higher excited states. A possible constraint that can limit the condensation fraction is the requirement that the population ${\rho }_i$ do not enter the stability island around the functional ground-state of the nonlinear Schrödinger equation. Since we do not know the equation for the boundary of this stability island we cannot perform realistic simulations, but we checked that simple approximation for the shape of the stability island indeed limit the condensation. For example, if the stability island imposes ${\sum }_{i=2}^N{\rho }_i>{\rho }^{*}$, then the condensate fraction will converge to $1-{\rho }^{*}$, with ${\rho }_2\to {\rho }^{*}$ in the limit $N\gg 1$.

To summarize, for all systems discussed above we see that the RJ condensate fraction can be high only at low temperature $T\approx {ϵ}_0/N$ being comparable with a typical energy level spacing when the BE thermal distribution (1) cannot be approximated by the RJ distribution (2). This shows the difference of RJ condensate from the Fröhlich condensate [50,51] existing at high temperature in quantum systems with BE thermal distribution (1).

5. Dynamical Thermalization for Defocusing Fiber NSE

5.1. Thermalization Features

Here we present the numerical results for the time evolution in the defocusing NSE (3) for typical values of moderate and weak nonlinearity $\beta =10,1,0.5$. Due to chaos and exponential instability of motion the numerical integration corrections in NSE lead to symmetry parity breaking and for an initial odd state with $\psi (x.y)=-\psi (x,-y)$ the even component of positive parity $\psi (x,y)=\psi (x,-y)$ appears with time. Due to this we state with an initial state which contains from the beginning at $t=0$ both parities using usually $\psi (t=0)=({\psi }_m+{\psi }_{m+1})/\sqrt{2}$.

The time evolution of the projection probability to the linear ground state ${\rho }_g\left(t\right)={\rho }_1\left(t\right)=|<\psi (x,y,t)|{\psi }_1{(x,y)>|}^2$ is shown for various initial states in Figure 12 at $\beta =10$. The initial state with a minimal total energy $E=16.2$ is $\psi (t=0)=({\psi }_1+{\psi }_2)/\sqrt{2}$ so that ${\rho }_g(t=0)=0.5$. For all other initial states we have ${\rho }_g(t=0)=0$. The results of Figure 12 show that for all initial states the probability ${\rho }_g$ approaches to its limiting value being approximately in the range $0.75\le {\rho }_g\le 0.9$. Thus the probability fraction of RJ condensate at the linear ground state is very high but still it does not reach so high values as in the case of RJS model (see Figure 8 with the maximal value ${\rho }_g={\rho }_1\approx 0.98$). We discuss the reason of this limiting value of ${\rho }_g$ at the end of this Section. For the states with higher initial energy, e.g. $E=142.1;169.8$ the values ${\rho }_g$ are a bit smaller compared to those of lower energy, e.g. $E=20.6$. We attribute this to the fact that it takes more time to reach steady-state RJ thermal distribution (2) from high energy states with a transfer of probability from high energy to the ground state. In Figure 12 we also show the expected values of ${\rho }_g$ according to the BE ansatz (1) and to the RJ ansatz RJ (2) assuming the total energy is $E=20$ (and total number of states $N=500$ for the RJ case). We have for the BE ansatz ${\rho }_g\approx 0.38$ while for RJ ansatz ${\rho }_g\approx 0.98$ being very close to unity thus showing a drastic difference between BE and RJ thermalization.

In global the results of Figure 12 demonstrate a good agreement with the RJ thermal condensate theory.

In Figure 13 we show the probability distribution $<{\rho }_m\left(t\right)>$ over linear eigenmodes ${\psi }_m$ averaged over a long time interval $500\le t\le 1000$ for the same initial states as in Figure 12 at $\beta =10$. For all initial states we have a high fraction of RJ condensate at the linear ground state ${\psi }_g={\psi }_1$. The probabilities ${\rho }_m$ at higher modes with $m>1$ approximately follow the RJ thermal distribution (2) which is characterized by a decay ${\rho }_m\sim 1/(E_m-\mu )$ shown by a dashed black curve for initial energy $E=20$. The BE ansatz at the same energy $E=20$ is also shown in Figure 12 by the dashed red curve and it is obviously very different from the presented numerical results. Even if there are visible fluctuations of ${\rho }_m$ values, which we attribute to not very long times reached in heavy numerical simulations, on average we show that the dynamical thermalization and chaos in the defocucing fiber NSE (3) leads to the RJ thermal distribution (2) providing a good description of numerical results including the formation of RJ condensate.

The Boltzmann H-theorem of 1872 [8], see also e.g. [11,48], states that during a thermalization process a system entropy $S\left(t\right)$ cannot decrease with time t, thus $S\left(t\right)$ can only grow with time or stay constant at its maximal value in a thermal equilibrium. In our system (3) with the dynamical thermalization we have a strikingly different situation. The results presented in Figure 14 obtained at the same parameters as in Figure 12, Figure 13 show that the system entropy $S\left(t\right)=-{\sum }_m{\rho }_m\left(t\right)ln{\rho }_m\left(t\right)$ initially grows with time, reaches a maximal value and then decreases up to a certain smaller value (being up to a factor 3 smaller than at maximum) corresponding to the RJ thermal distribution (2) with RJ condensate. Our qualitative explanation of this surprising behavior of $S\left(t\right)$ is the following: initially $S(t\approx 0)\approx ln2$ since initially two linear eigenmodes are populated, then at higher times transitions to other linear eigenmodes are produced by nonlinear $\beta$ term and S is increased reaching its maximal value, but then a slow thermalizarion due to chaos leads to the RJ thermal distribution (2) where the fraction of RJ condensate is close to unity (see Figure 12), since RJ condensate fraction is so high the entropy S of steady-state has a relatively small value being significantly smaller than its maximal value. We argue that there is no contradiction with the Boltzmann H-theorem since it is stated for systems with only one energy integral of motion, while in our system we have two integrals of energy and norm.

In Figure 15 we compare the dependence of entropy S values, averaged over a certain time interval at low and high times, on initial energy of linear modes $E_m$. On short times $30\le t\le 40$ this dependence $S\left(E_m\right)$ is close to those of one given by the BE thermal distribution (1) while at longer times $500\le t\le 1000$ it is more close to those of RJ thermal distribution (2). The reason of this is related to the unusual entropy time dependence discussed for Figure 14 where S grows at initial times and later decays to a smaller value induced by the RJ thermal condensate. On the basis of obtained results we argue that the result of [44,70] that the dependence $S\left(E\right)$ is close to those of the BE distribution (1) instead of RJ one (2) is related to significantly shorter times reached there. Here the FREEFEM codes allow us to reach times being by a factor 30 longer where the effects of dynamical RJ thermalization and RJ condensate start to play the dominant role.

5.2. Stability Island Near The Ground State Linear Eigenmode

The obtained results exposed in Figure 12. Figure 13 show that the fraction of RJ condensate does not exceed 80%-90% even if for the RJS model (see Figure 8) this fraction goes to unity with the increase of total number of linear system modes N. Of course, such an effect can be related to a finite number of linear eigenmodes N effectively populated during the NSE (3) evolution. However, we argue that there is another reason behind the bounded fraction of RJ condensate. In fact this reason is related to existence of a stationary self-consistent solution of the NSE (3). To find this solution ${\psi }_s(x,y)$ we put the time derivative in (3) equal to zero. Then for the remained stationary Hamiltonian $H_{sf}(x,y)=-\Delta +\beta {\left|{\psi }_s(x,y)\right|}^2$ we find the self-consistent solution ${\psi }_s$. This is done by the recursive iterations ${\psi }_g\left(k\right)=\alpha {\psi }_g(k-1)+(1-\alpha ){\psi }_g(k-2)$ where ${\psi }_g(k-1)$ is the ground state of the Hamiltonian $H_{sf}(k-1)=-\Delta +\beta {\left|{\psi }_g(k-2)\right|}^2$. With the parameter $\alpha$ being close to unity, e.g. $\alpha =0.9$, this iterative process rapidly converges to the self-consistent stationary ground state ${\psi }_s(x,y)$ of the NSE (3). The same procedure allows to find other self-consistent eigenstates ${\psi }_{s2},{\psi }_{s3}\dots$ for higher modes.

The comparison of this stationary self-consistent state ${\psi }_s(x,y)$ at $\beta =10$ with the linear ground state ${\psi }_1(x,y)$ at $\beta =0$ is shown in Figure 16. The projected probability of ${\psi }_s$ on ${\psi }_1$ is 0.991 so that both states have similar shape but ${\psi }_s$ has also about 0.009 probability at other high energy linear eigenstates ${\psi }_m$.

Even if we find the stationary self-consistent eigenstates at finite $\beta =10$ the main question is if these states are stable or not. To determine their stability we follow the time evolution of NSE (3) for an initial state being a slightly perturbed self-consistent state. The time evolution of the deviation error, defined as $1-{|<\psi \left(t\right)|\psi (t=0)>|}^2$, is shown in Figure 17. The results show that the self-consistent ground state ${\psi }_s={\psi }_{s1}$ is stable in respect to perturbations while the excited self-consistent state ${\psi }_{s3}$ is exponentially unstable and in this case the deviation is growing exponentially with time. These results show that the self-consistent ground state ${\psi }_s$ is stable, e.g. at $\beta =10$, and hence there is a certain stability region around this state. In contrast, the higher energy self-consistent stationary states are unstable. The initial states started at high linear modes ${\psi }_m$ at $m>1$ are located in a chaotic component and thus they cannot populate the integrable stability region located around the stable self-consistent ground state ${\psi }_s$. We argue that this is the main reason due to which the fraction of RJ condensate remains bounded to 80%-90% (see Figure 12, Figure 13).

5.3. Quasi-Integrable KAM Regime

According to the KAM theorem [21,22,23,24,42] we expect that at small nonlinearity $\beta \ll {\beta }_{ch}$ the NSE dynamics is integrable in the major part of system phase space, even if certain small chaotic regions can remain, e.g. around thin separatrix layers with the Arnold diffusion. Indeed, the results of Figure 18 show that at small $\beta$ the entropy $S\left(t\right)$ is not growing with time showing only regular oscillations at $\beta =1$ or almost constant small value $S<1$ at $\beta =0.5$. Even if S values at $\beta =1;0.5$ are comparable with S value at long time $t=1000$ for $\beta =10$ the situation is qualitatively different: at $\beta =10$ the dynamics is chaotic leading to RJ thermalization with a high fraction of RJ condensate (that gives small S values at $t=1000$) while at $\beta \le 1$ the probabilities ${\rho }_m\left(t\right)$ are simply regularly oscillating near there initial values as it is shown in Appendix Figure A3.

It is not easy to determine the chaos border for the NSE (3) in a chaotic billiard. Similar to the arguments [44,72] we estimate that the developed chaos appears in the system when the nonlinear energy spread $\delta E\sim {\beta \left|\psi \right|}^2\sim \beta /A$ becomes larger than the level spacing $\Delta \approx 4\pi /A$ that leads to the chaos border:

$$\beta >{\beta }_{ch}\sim 4\pi .$$

(18)

This is an approximate estimate while the obtained results described above show that developed chaos and dynamical thermalization appear already at $\beta =10$. It is not so easy to determine the exact position of the chaos border due to complexity of various regimes present in the fiber NSE (3) and thus this question requires further studies.

6. Wave Collapse for Focusing Fiber NSE

It is known that for the focusing NSE on an unrestricted 2D plane the wave collapse takes place in a finite time if initial total energy of the system is negative [45,46]. Here we show that the collapse in a finite D-shape billiard can take place at moderate negative nonlinearity $\beta <0$ and certain positive energy. A typical example is shown in Figure 19 at $\beta =-12$ where the initial state at $t=0$ is taken as the ground state of linear system. Here the collapse happens at rather short time $t_{col}\approx 0.16$ even if the initial total energy at this $\beta =-12$ is positive being $E=2.076$. At $t_{col}$ the kinetic energy of the system $E_{kin}=<\psi \left(t\right)|-\Delta |\psi \left(t\right))>$ diverges growing to infinity. Indeed, approaching the collapse time $t_{col}$ the wave function is localized in smaller and smaller area $A_{col}$ in $(x,y)$ plane and due to norm conservation ${\left|\psi (x.y.t)\right|}^2\sim 1/A_{col}$ is growing so that negative nonlinear energy is also increasing as ${\beta \left|\psi (x.y.t)\right|}^2\sim \beta /A_{col}<0$. Thus due to the conservation of total energy we have divergence of kinetic energy $E_{kin}\sim \left|\beta \right|/A_{col}$. Special checks ensure that these results are independent of the numerical integration parameters $N_F$ and $\Delta t$. For the initial state taken at the next energy level with $\psi (t=0)={\psi }_2(x,y)$ at $m=2$ with the total energy $E=10.4$ there is no collapse and kinetic energy $E_{kin}$ simply oscillates around its initial value. A similar collapse is also present at $\beta =-10$.

We saw above that for $\beta >0$, e.g. $\beta =10$, due to dynamical thermalization with RJ thermal distribution (2), a significant fraction of total probability is accumulated in the RJ condensate at the ground state ${\psi }_1(x,y)$ with $m=1$. This probability ${\rho }_g=P_0$ can be around 80%-90% (see Figure 12, Figure 13). For negative values of $\beta$ we expect that the dynamical thermalization leads also to emergence of RJ condensate with a high condensate probability fraction. Thus it is important to see if the collapse can take place for an initial state with a high probability at the linear ground state and a relatively small probability at high energy state. To study such a situation we choose the initial state being $\psi (x,y,t=0)=({\psi }_1+C{\psi }_{181})/\sqrt{1+C^2}$ with $C\sim 0.18$. Thus at $t=0$ there are about 96.8% of probability in the ground state $m=1$ and 3.2% at the excited state $m=181$. For $C=0.18$ the results show that there are time oscillations of $E_{kin}$ with a significant increase of $E_{kin}$ but the collapse is avoided due to presence of high energy component of wave function (see top panel of Figure 20). However, for $C=0.15$ still the collapse takes place at $t_{col}\approx 0.18$ even if the initial total energy is rather high $E_{tot}=20.66$ (see Figure 19).

In the pre-collapse regime at $C=0.18$ there is an interesting behavior of wave function $\psi (x,y,t)$ that can be seen by its projection probability ${{\rho }^{*}}_n$ on the eigenbasis $1\le n<\infty$ of the Schrödinger equation (3) in which the nonlinear term is considered as the instantaneous effective attractive potential $U_{eff}=\beta {\left|\psi (x,y,t)\right|}^2\le 0$. The dependence of ${{\rho }^{*}}_n$ on time t is shown in the bottom panel of Figure 20. At small times $t<0.1$ there is about ${{\rho }^{*}}_1=0.95$ to $0.915$ probability located at the ground state $n=1$ of the instantaneous potential $U_{eff}(x,y,t)$. But $U_{eff}$ is changing with time that leads to a transition of almost all fraction of ${{\rho }^{*}}_1$ to other eigenstate at $n=2$ at higher moment of time t. Then with growth of t there are transitions to high and higher states n of instantaneous potential $U_{eff}$. The maximal n values (up to $n=16$) are reached at the minimal size of wave function corresponding to the maximal values of $E_{kin}$ in the top panel. Then, when the packet squeezing passed, its maximal value and the packet size starts to increase the population ${{\rho }^{*}}_n$ starts again go to minimal n values returning back to almost all probability at the ground state $n=1$ of $U_{eff}$ when the packet returns of its maximal site (minimal of $E_{kin}$ in the top panel). It is interesting to note that the transitions between states n happens in diabatic manner with almost all probability transferred from one n value to another one. These results show that in the pre-collapse regime there is rather nontrivial time evolution of wave function with its oscillating size. It would be interesting to study this behavior in more detail but this goes outside of scopes of this work.

Thus we show that the wave collapse in a focusing NSE inside a chaotic billiard can take place even at rather highly positive values of total system energy when a probability fraction on high energies is relatively small. Thus there should be a certain boundary in the space of total energy and probability on high levels ($C^2$ in Figure 19, Figure 20). In global this boundary determines the regime of collapse stability in respect of perturbations produced by high energy waves. Since the collapse remains stable only at rather small probability at high energies we expect that in the regime of dynamical thermalization and appearance of RJ condensate there will be no wave collapse at moderate negative $\beta$ values. The results of dynamical thermalization at negative $\beta$ are presented in the next Section.

7. Dynamical Thermalization in Focusing Fiber NSE

For negative nonlinearity $\beta <0$ in (3) the dynamical thermalization has the features similar to the case $\beta >0$ with certain specific aspects related to the wave collapse discussed above. Thus in Figure 21 we show the time dependence of projection probability ${\rho }_g$ on the linear ground state ${\psi }_1$ for a typical case with $\beta =-10$ and initial state $\psi (t=0)=({\psi }_m+{\psi }_{m+1})/\sqrt{2}$ at $m=11$. At long times t it approaches to ${\rho }_g\approx 0.8$ showing the formation of RJ condensate. This behavior is similar to those at $\beta =10$ in Figure 12 with the same initial state. In Figure 21 we add also the projection probability ${{\rho }^{*}}_g$ on the ground state of the instantaneous effective potential $U_{eff}=\beta {\left|\psi (x.y.t)\right|}^2<0$ as it was done in Figure 20 bottom panel for the case of wave collapse. At long times when the steady-state is established we obtain ${{\rho }^{*}}_g\approx 0.9$ being higher than ${\rho }_g\approx 0.8$. Thus a higher fraction of RJ condensate is captured by attractive instantaneous potential $U_{eff}$. However, these fractions of RJ condensate are still not sufficiently high and there is no wave collapse. We also find the similar ${\rho }_g\left(t\right)$ time dependence for other initial states $\psi (t=0)=({\psi }_m+{\psi }_{m+1})/\sqrt{2}$ at $m=5$ and $m=7$ (see Appendix Figure A4).

In Figure 22 we show the probabilities ${\rho }_m$ averaged in the time interval $400\le t\le 800$ with dashed curves showing the RJ ansatz (2) and BE ansatz (1) corresponding to energy E of the initial state (for RJ case we assume the total number of levels $N=500$). As for the case of Figure 13 at $\beta =10$ we find that the RJ ansatz provides a good description of the steady-state distribution even if the RJ fraction ${\rho }_g$ is higher comparing to the numerically obtained value ${\rho }_g\approx 0.8$ (due to that the RJ curve is located below the numerical ${\rho }_m$ values but the energy dependence ${\rho }_m$ on $E_m$ remains the same). Also the results show that the BE ansatz is not working. The similar results at $\beta =-10$ are shown in Appendix Figure A5 for other initial states $\psi (t=0)=({\psi }_m+{\psi }_{m+1})/\sqrt{2}$ at $m=5$ and $m=7$.

The time dependence of entropy $S\left(t\right)$ at $\beta =-10;-0.5$ and initial state $\psi (t=0)=({\psi }_m+{\psi }_{m+1})/\sqrt{2}$ at $m=11$ is shown in Figure 23. This dependence is similar to those in Figure 14 at $\beta =10$. At small $\beta =-0.5$ the entropy remains practically unchanged corresponding to the KAM regime. In this quasi-integrable KAM regime at $\beta =-0.5$ only regular probability oscillations in time (see Appendix Figure A6).

Thus we conclude that the process of dynamical thermalization is very similar for positive and negative values of nonlinearity $\beta$ in (3). For significantly negative $\beta$ value the collapse can take place even at positive total energies but then the initial state should have probability at the linear ground state being very close to unity.

8. Vortexes in Defocusing Fiber NSE

It is well known that the defocusing NSE in 2D unrestricted plane have vortexes at strong nonlinearity $\beta$ [65,66]. The properties of vortexes have been studied theoretically, numerically and experimentally with superfluidity in liquid helium [65,66] and with light fluid in nonlinear media of atomic vapor [62]. It is known that an effective size of vortex scales in dimensionless units as $1/\sqrt{\beta }$ [65,66]. A good agreement between numerical NSE modeling and experimental results has been demonstrated in [62].

Here we study the properties of vortex in a chaotic D-shape billiard in the frame of defocusing NSE (3). In the results presented above the initial state $\psi (x,y,t=0)$ was always chosen as a certain combination of linear modes. However, to study the vortex dynamics we need to start with initial state having initial vorticity. Thus we start with an initial state shown in Figure 24, it is close to a vortex field in an open 2D plane with $\psi (r,\theta )\propto rexp\left(i\theta \right)$ [65,66]. This state has a significant projection probability $P_1$ on the first excited linear mode ${\psi }_2(x,y)$ at $m=2$. Thus the high value of probability $P_1$ indicates that the vortex is well survived in time even if the classical dynamics of rays is chaotic in the billiard at $\beta =0$. In contrast, if the projection probability $P_0$ on th ground state linear mode ${\psi }_1(x,y)$ at $m=1$ becomes high and $P_1$ becomes low then this indicates that the vortex is destroyed by quantum chaos and nonlinearity.

The time evolution of projection probabilities $P_0,P_1$ for vortex time dynamics is shown in Figure 24 (left panel) for moderate nonlinearity $\beta =10$. For a small integration time step $\Delta t=dt={10}^{-4}$ and high number of FREEFEM net nodes $N_F=136241$ we have high numerical conservation of norm and total energy integral. For time $t\ge 80$ the results show that $P_0$ probability starts to grow significantly while $P_1$ probability starts to drop rapidly that tells about the complete destruction of vortex at times $t>80$. This is confirmed by color vortex images in the right panel of Figure 24: vortex structure is well visible at $t=10$ both for wave function amplitude $\left|\psi \right)x,y,t\left)\right|$ and its phase $arg\psi (x,y,t)$ showing clear phase rotation. In contrast at $t=100$ both $\left|\psi \right(x,y,t\left)\right|$ and $arg\psi (x,y,t)$ show disappearance of vortex. For times $t>80$ the vortex is destroyed and the dynamical RJ thermalization takes place for nonlinear field $\psi (x,y,t)$. Indeed, the results for $P_0$ show its significant growth so that its value becomes close to the probability fraction of RJ condensate and $P_1\approx 0.8-0.9$ (see Figure 12). At the same time the projection probability $P_1$ drops significantly similar to the cases shown in Figure 13. Thus the dynamical RJ thermalization works also for the initial vortex states at moderate nonlinearity $\beta =10$.

It is interesting to note in the left panel of Figure 24 that at smaller number of nodes $N_F=37906$ the evolution of vortex is more stable with higher projection probability $P_1$ at $t\approx 100$ comparing to the case with significantly higher integration parameters $dt={10}^{-4}$ and $N_F=136241$. We interpret this by assuming that the lower value $N_F=37906$ leads to a certain effective roughness of induced effective potential that reduces velocity of vortex center decreasing the quasi-particle excitation thus increasing vortex life-time.

The dynamics of vortexes at a moderate nonlinearity $\beta \sim {\beta }_{ch}$ tends to approach the dynamical thermalization even if the votrex life time $t_v$ can be relatively high. At high $\beta$ values the vortex size decreases as $1/\sqrt{\beta }$ and its interactions with billiard border becomes weaker with the growth of vortex life-time $t_v$ as it is show in Figure 25 where the projection probability $P_1$ remains high for longer and longer times. However, we find here a striking effect when this vortex life-time $t_v$ changes abruptly from $t_v\approx 70$ for $\beta =44.7$ to $t_v>600$ for $\beta =45$. This sharp transition in $t_v$ is also well seen in Figure 26 where the snapshots of nonlinear field $\psi (x,y,t)$ are shown at $t=200$: for $\beta =44.7$ the vortex is completely destroyed while for $\beta =45$ it is well preserved. The vortex time evolution is also shown in Supplementary Material (SupMat) video (file videovortex.mp4) for vortex time evolution at $\beta =44$ (left panel, vortex is destroyed at maximal shown $t\approx 30$) and at $\beta =45$ (right panel, vortex is well survived at $t>30$). This video shows that at $\beta =44$ the center of vortex moves more rapidly compared to those at $\beta =45$. Thus we make a superfluid transition conjecture: at $\beta <{\beta }_{sf}\approx 45$ the velocity of vortex center is higher than the critical superfluid velocity of light fluid and it starts to radiate quasi-particles breaking superfluidity, while for $\beta >{\beta }_{sf}\approx 45$ this velocity is below critical superfluid velocity and the vortex life-time is increased enormously or may become even infinite. The verification of this superfluid transition conjecture requires further studies which are beyond the scope of this work.

9. Nse Dynamics at Long Times

The question about asymptotic properties of NSE dynamics in the D-shape chaotic billiard at long time scales is rather complex. The NSE (3) can be rewritten in the basis of linear eigenmodes ${\psi }_m(x,y)$ of the billiard (at $\beta =0$) where it reads:

$$i\partial C_m/\partial t=E_m{C}_m+\beta \sum _{m_1{m}_2{m}_3}{V}_{m{m}_1{m}_2{m}_3}{C}_{m_1}{C}_{m_2}^{*}{C}_{m_3}.$$

(19)

Here $C_m\left(t\right)$ are wave function amplitudes in the eigenbasis ${\psi }_m$ so that $\psi (x,y,t)={\sum }_m{C}_m\left(t\right){\psi }_m(x,y)$ and interaction matrix elements induced by nonlinearity are $V_{m{\mathbf{m}}_{\mathbf{1}}{m}_2{m}_3}={\int }_{x,y}{{\psi }^{*}}_m{\psi }_{m_1}{{\psi }^{*}}_{m_2}{\psi }_{m_3}dxdy$.

This type of equation appears also in the models on nonlinear perturbation of quantum chaos [72], Anderson localization (DANSE model for Discrete Anderson Nonlinear Schrödinger Equation) [73,74,75,76] and Random Matrix Theory (RMT) [26] (NLIRM model). In these models the norm and total energy are exactly conserved as for NSE model (3). For DANSE type models [72,73,74,75,76] the linear eigenstates are exponentially localized due to quantum interference and Anderson localization so that the coupling matrix elements are $V_{m{m}_1{m}_2{m}_3}\sim 1/{\xi }^{3/2}$ for states ${\psi }_m,\dots {\psi }_{n_3}$ located on a distance of localization length $\xi$ from each other on a lattice, while outside of this range the values of V drop exponentially. Here $\xi$ is a number of lattice states in a localization length and for a d-dimensional lattice with disorder we have $\xi \sim {\ell }^d$ [74].

For such type DANSE models it was shown that for a nonlinearity above a certain chaos border $\beta >{\beta }_c$ there is a slow subdiffusive spreading of probability over the lattice sites n with the width growing as ${\left(\Delta n\right)}^2\sim t^{\nu }$ with the subdiffusive exponent $\nu \approx 0.3-0.4$ [72,73,74,75,76]. This growth continues during all computational times up to extremely long times $t\sim {10}^{10}$.

It is interesting to note that for the DANSE models considered in [72,73,74,75,76] the spectrum of linear modes is bounded in a certain finite energy band. However, the DANSE type model studied in [77] has a Stark lattice with on-site energies growing linearly with site n index $E_n\sim fn$ and still the unlimited spreading was present up to maximal times $t={10}^8$ with an exponent $\nu \approx 0.24-0.30$ (in a chaotic regime with $\beta >{\beta }_c,f<f_c$). Due to energy growth with site number the model [77] is similar to one studied here (3) since in a 2D billiard energies of modes are also growing with mode index $E_m\propto m$. Thus due to such a similarity between the model [77] and the present fiber model (3) one could expect to have unlimited norm spreading with time up to higher and higher modes m. However such an unbounded scenario has certain difficulties: indeed, if the total number of populated linear modes $N\sim \Delta m\left(t\right)\sim t^{\nu /2}$ grows with time then due to RJ condensation a main part of the total norm should be in the condensate phase being mainly at the linear ground state. From the view point of RJS model this norm should become closer and closer to unity with the growth of N (see Figure 8). Due to presence of stability island with self-consistent ground state we see that in the NSE fiber model (3) RJ condensate has not more than approximately 90% of the total norm (see Figure 12) so that only about 10% of norm can propagate to higher and higher N values. In the case of Stark ladder model [77] we have there about 50% of the norm propagating to the high energy modes, counting from the initial state location (and about 50% of norm moving to the low energy modes). Thus there can be a certain similarity between the Stark model [77] and our fiber NSE model (3). However, in [77] the coupling matrix elements V are approximately comparable for all modes (centered within a localization length $\xi$), independently of m values, while for the fiber case the dependence of V on m requires further analysis. Indeed, at high $E_m$ (or m) values the linear modes become more and more oscillating (see Figure 3) and thus we may expect that a number of effective components $\xi$ is growing with $E_m$ and thus couplings V are decreasing. A simple estimate based on the results for rough billiards [78] would give $\xi \sim \sqrt{E_m}\sim \sqrt{m}$ and thus $V\sim 1/{\xi }^{1/2}\sim 1/{E_m}^{1/4}$, but this expectation requires separate verification. If V decreases with $E_m$ then effective nonlinearity also decreases and at some high energies $E_m$ an effective nonlinearity $\beta V$ becomes below a chaos border that leads us to a bounded scenario with only finite number of effectively populated modes. At present we cannot present firm arguments in favor of bounded scenario of excitation to high energies or unbounded one.

Also the results presented in [26] for RMT with growing energy diagonal term (see SupMat Figure S14 in [26]) show that even in a system with $N=64$ modes the thermal RJ distribution is reached only at such high times as $t\approx {10}^8$ ($t\approx {10}^6$ is not enough). Such times cannot be reached in numerical simulations of fiber NSE (3). Such high times are also out of reach with fiber experiments.

For initial states with vortexes our results in Figure 24, Figure 25, Figure 26 show that for ${\beta }_{ch}<\beta <{\beta }_{sf}$ a vortex state is destroyed and we expect that the dynamical RJ thermalization takes place at long times which however are difficult to reach in numerical simulations due to high life time of vortex. For the superfluid phase at $\beta >{\beta }_{sf}\approx 45$ the question about vortex life time remains open.

10. Experimental Parameters for A Quantum Chaos Fiber

Here we discuss the requested experimental parameters for D-shape fiber to have sufficiently strong nonlinearity $\beta$ to reach regime of dynamical thermalization. We start from experimental conditions already realized in [1,2] with a circle fiber cross-section where the classical dynamics of linear system at $\beta =0$ in (3) is integrable. These fiber experiments were done at high laser beam intensity and the theoretical description of the complex field envelope $A(x,y,z,t)$ expressed in $[\sqrt{W}/m]$ is described by the NSE type equation in which we keep notations of [2]:

$${\displaystyle \frac{\partial A}{\partial z}}-(i/2k_0){{\nabla }^2}_{\perp }A=i\gamma {\left|A\right|}^{2}A.$$

(20)

Here z is distance along the fiber, corresponding to time in (3), $k_0=\omega n_{co}/c$ is the light wave vector $k_0=2\pi /\lambda$, $\gamma =\omega n_2/c$ is the fiber nonlinear coefficient with nonlinear Kerr index $n_2=3.2\times {10}^{-20}{m}^2/W$; $n_{co}\approx 1.5$ is the refractive index. Here we omit the term with time derivative $d^{2}A/d{t}^2$ which gives only linear in A contribution in (20) and assumed that the refractive index $n_{co}$ is constant inside the fiber of core radius $r_c=30\mu m$. Thus from (20) from the above typical parameters and a typical wave length $\lambda =1000nm$ we obtain the expression for coefficient $\beta$ in the NSE (3):

$$\beta =2n_{co}{n}_2{(2\pi /\lambda )}^2{r_c}^{2}I[W/c{m}^2],$$

(21)

that gives $\beta \approx 0.32$ for typical parameters of [2] with $I=10GW/c{m}^2$, $n_2=\times {10}^{-20}{m}^2/W$, $n_{co}=1.5$, $\lambda =1000nm$, $r_c=30\mu m$. In (20) the length unit is $z_0={r_c}^2/\lambda \approx 0.1cm$ for the above values $r_c$, $\lambda$. Thus for a fiber length of 1 meter we have the dimensionless length of $z/z_0\approx 1000$ that corresponds to the dimensional time interval t of NSE (3) being also 1000. In our numerical simulations of NSE (3) we are reaching similar time scales $t\approx 100$ that is comparable with the fiber length of $0.4m$ in [2].

We suppose that increasing by a factor 3 the parameters $r_c\approx 90\mu m$ and $I\approx 30GW/c{m}^2$ in (21) it would be possible to have nonlinear parameter $\beta \approx 10$ at which dynamical thermalization sets in the D-shape quantum chaos fiber described by the NSE (3).

11. Discussion and Conclusion

Our numerical and analytical studies of the NSE fiber light fluid dynamics in a chaotic D-shape billiard show that it is characterized by various unusual and rich phenomena. Above the chaos border at $\left|\beta \right|>{\beta }_{ch}\sim 10$ the dynamical thermalization leads to the RJ thermal distribution (2) with a formation of the RJ condensate in the ground state of linear system which fraction has about 80-90% of total probability. Certain similarities between RJ and Froöhlich condensates are discussed. We show that at $\beta >0$ near the ground state there is a self-consistent state with a certain island of stable motion around it. In contrast to the Boltzmann H-theorem [8,11,48] in the chaotic regime the entropy time evolution $S\left(t\right)$ is characterized by initial increase at short times and later relaxation from its maximal value to a significantly smaller value on long times in the steady-state regime that is explained by the formation of RJ condensate. We point that the presence of two integrals of motion, energy and norm, is at the origin of such a difference from the H-theorem formulated for systems with only one energy integral of motion. In absence of nonlinearity our system belongs to the systems of quantum chaos, thus having many features as in the Random Matrix Theory. On these grounds we argue that the dynamical thermalization induced by nonlinearity represents a generic features of such thermalization similar to those discussed for nonlinear perturbation of Random Matrix theory in [26]. We also note that even if the phenomenon of RJ condensation in nonlinear systems has been discussed for multi-mode fibers [58,59,60,61] its origin was not associated with the emergence of dynamical chaos above a certain chaos border below which the RJ thermalization and RJ condensate are absent and dynamics is quasi-integrable as in the KAM theory. We attribute the so-called self-cleaning phenomenon observed in multi-mode optical fibers [1,2] to the formation of RJ condensate induced by dynamical thermalization and chaos. A pictorial image of time development of dynamical thermalization and RJ condensate formation is shown in Figure 27.

In addition to dynamical thermalization at certain conditions the focusing fiber NSE (3) at $\beta <0$ can have the regime with wave collapse which in contrast to the case of open space [45,46] can happen even at relatively high positive energies (see Figure 19, Figure 20). In such regime there is a nontrivial interplay between chaotic component and wave collapse. However, the RJ condensate emerging due to dynamical thermalization does not lead to the collapse formation.

We also show that for the defocusing fiber NSE (3) there are long living excitations in the form of vortexes which finally are expected to relax to the RJ thermal state for nonlinearity ${\beta }_{ch}<\beta <{\beta }_{sf}$. However, for a strong nonlinearity $\beta >{\beta }_{sf}\approx 45$ the vortex life time is strikingly increased corresponding to the transition to the superfluid phase of light fluid.

We hope that these interesting and generic phenomena described here can be studied experimentally with multi-mode fibers similar to those of [1,2,3,4,5,6,7] (namely with a quantum chaos fiber cross-section, e.g. D-shape as in Figure 1). In fact certain studies of D-shape fiber have been reported recently in [37] but without discussion of dynamical thermalization. Also the quantum fluid experiments with laser beam propagating in a nonlinear medium formed by atomic vapor, like in [62], can be used to test interplay of nonlinearity, quantum chaos and dynamical thermalization.