{'fraction_non_alphanumeric': 0.03504167953071936, 'fraction_numerical': 0.01589996912627354, 'mean_word_length': 4.493216506500848, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical thermalization of Bose-Einstein condensate in Bunimovich stadium

We study numerically the wavefunction evolution of a Bose-Einstein condensate in a Bunimovich stadium billiard being governed by the Gross-Pitaevskii equation. We show that for a moderate nonlinearity, above a certain threshold, there is emergence of dynamical thermalization which leads to the Bose-Einstein probability distribution over the linear eigenmodes of the stadium. This distribution is drastically different from the energy equipartition over oscillator degrees of freedom which would lead to the ultra-violet catastrophe. We argue that this interesting phenomenon can be studied in cold atom experiments.

Introduction. -Starting from the famous Fermi-Pasta-Ulam problem [1,2] the interest to understanding of thermalization process in dynamical systems with a finite number of degrees of freedom is continuously growing (see e.g. [3,4]). At present the experiments with cold Bose gas in atom traps and optical lattices open access to experimental investigations (see e.g. [5,6]) stimulating the theoretical and experimental interest to this phenomenon (see e.g. [7,8]).
The numerical analysis of dynamical thermalization in disordered nonlinear chains has been started recently showing that the quantum Gibbs distribution can appear at a moderate nonlinearity contrary to usually expected energy equipartition over linear modes [9,10]. Thus it is interesting to analyze the dynamics of a Bose-Einstein condensate (BEC), described by the Gross-Pitaevskii Equation (GPE), in a chaotic billiard where the quantum evolution corresponds to a regime of quantum chaos and energy levels statistics described by the Random Matrix Theory [11][12][13]. As an example we use a de-symmetrized Bunimovich stadium where the classical dynamics is fully chaotic (see e.g. [12,14]). We note that the chaotic optical billiards, created by a laser beam and containing cold atoms, have been already studied experimentally [15,16] and hence our model can be investigated experimentally.
Model description. -The model is described by the GPE for BEC in the de-symmetrized Bunimovich stadium billiard with Dirichlet boundary conditions: ih ∂ψ( r, t) ∂t = −h 2 2m ∆ψ( r, t) + β|ψ( r, t)| 2 ψ( r, t) where we considerh = 1, m = 0.5. The height of the stadium is taken as h = 1 and its maximal length is l = 2 (see Fig. 1). Thus the average level spacing is ∆ ≈ 4π/A ≈ 7.04 where A is the billiard area. At β = 0 the numerical methods of quantum chaos allows to determine efficiently about million of eigenenergies of linear modes and related eigenmodes [17]. For comparison, we also consider the case on a rectangular billiard with approximately the same area as for stadium and with the golden mean ratio l/h = (1 + √ 5)/2, h = 1. We note that for model (1) the spectrum of Bogoliubov excitations of BEC in a Bunimovich stadium had been studied in [18], but the question of dynamical thermalization has been not addressed there. We also point that the model (1) is described by the partial differential equation (continuous variables) thus being significantly more complex than the case of nonlinear chains studied in [9,10]. Indeed, even the prove of the existence of solution in (1) is a nontrivial mathematical problem which still remains open (e.g. the ultra-violet catastrophe would imply the absence of solution). In the following we restrict our analysis to the GPE p-1 arXiv:1505.05130v2 [cond-mat.quant-gas] 4 Sep 2015 case not going beyond this description. The time evolution of (1) is integrated by small time steps of Trotter decomposition of linear and nonlinear terms with a step size going down to ∆t = 4 × 10 −5 to suppress nonlinear growth of high modes. We use N c = 1076(1085) linear modes φ m of stadium (rectangular) for linear part of time propagation doing the nonlinear step with β term in the coordinate space with N p = 11207(12816) points inside the billiard. The lattice points are given by 79 × 79 = 6241 equidistant x-y coordinates for the square part of the stadium billiard (rectangular billiard), and rays with equidistance in angles for the circular part. The change of basis from coordinates to energies (and viceversa) is given by a unitary matrix in double precision. Similar to [19], a special aliasing procedure 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 largest value β = 20 and better than 1% at lower β values. At any step the wavefunction is expanded in the basis of linear modes φ m so that ψ(x, y, t) = m C m (t)φ m (x, y). The averaging of probabilities w m (t) = |C m (t)| 2 ( m w m = 1) over time gives the average probability distribution ρ m = |C m | 2 t . We note that the quantum evolution of GPE (1) has been studied in the frame of quantum turbulence for a rectangular billiard [20] and for a 3D-cube [21]. However, in these studies there is energy injection/absorption at low/high modes to generate the Kolmogorov energy flow in space of linear modes (see e.g. [22,23]). We also note that the time evolution of wave packet for the GPE in Bunimovich stadium had been simulated in [24] but only a spacial distribution had been considered there. In contrast to these studies we consider only unitary or Hamiltonian evolution given by (1) being interested in the distribution properties of probabilities ρ m over linear eigenmodes. In this respect our approach is different from other studies where the analysis had been concentrated on space fluctuations (see e.g. [25]). Also, as we will see below, there is a significant difference for the GPE evolution in chaotic and rectangular billiards.
Time evolution. -Examples of time evolution for two initial eigenmodes m = 5, 10 are shown in Fig. 1 at β = 10 (video is available at [26]). They show that, due to nonlinearity inside the stadium, there are complex irregular oscillations of wavefunction with time.
Another representation is obtained by considering the time evolution of probabilities w m (t) in the basis of linear modes shown in Fig. 2. At moderate value β = 2 the probability remains located only in a few modes without thermalization and spreading over many modes. For larger value of β = 10 the nonlinear spearing over modes goes in a more efficient way with many excited modes. Thus the dynamical thermalization is expected to be absent at small or moderate β < β c ∼ 1 while for large nonlinearities β ∼ 10 > β c we may expect the emergence of dynamical thermalization.
For the initial state m = 10 we have the approxi- 280 is a velocity of classical particle and a time interval between collisions is approximately τ col ≈ h/v ≈ 0.06. Thus during the time t = 40 we have approximately N col = t/τ col ≈ 670. Dynamical thermalization is reasonably achieved for time interval t ∈ [20, 40] as it is visible in middle and bottom panels of Fig.2 where w m (t) have thermalized like distribution with β = 10.
Bose-Einstein thermal distribution. -To characterize the dynamical thermalization in more detail we assume that a moderate nonlinearity acts as a certain thermalizer which drives the system to a thermal equilibrium over quantum levels of the stadium. At the same time we assume that the nonlinear term is not very strong so that it does not affect significantly the average linear eigenenergies. Indeed, on average we have β|ψ| 2 ≈ β/A ≈ ∆ for β ∼ 10, so that indeed, the nonlinear energy shift is moderate at such values of β.
Thus we expect that nonlinearity generates a dynamical thermalization over the quantum billiard energy levels. In such a case we should have the standard Bose-Einstein distribution ansatz over energy levels E m [27]: where E g = 13.25 is the energy of the ground state, T is the temperature of the system, µ(T ) is the chemical potential dependent on temperature. The values E m are the eigenenergies of the stadium at β = 0. The parameters T and µ are determined by the norm conservation ∞ m=1 ρ m = 1 (we have only one particle in the system) and the initial energy m E m ρ m = E. The entropy S of the system is determined by the usual relation [27]: S = − m ρ m ln ρ m . The relation (2) with normalization condition determines the implicit dependencies on temperature E(T ), S(T ), µ(T ).
The advantage of energy E and entropy S is that both are extensive variables, thus they are self-averaging and due to that they have reduced fluctuations. Due to this feature S and E are especially convenient for verification of the thermalization ansatz. To check this ansatz we start from an initial linear mode m which corresponds to the system energy E ≈ E m and follow the GPE time evolution of probabilities w m (t) determining the value of entropy S from obtained average probabilities ρ m . Considering the initial states with 1 ≤ m ≤ 50 we obtain the numerical dependence S(E) shown by symbols in Fig. 3. This dependence is compared with the analytic curve following from the Bose-Einstein ansatz (2) which gives the dependencies E(T ) and S(T ) and hence provides the analytic dependence S(E) shown by the red curve in Fig. 3.
The data of Fig. 3 show that even at large β = 20 there is no thermalization for the rectangular billiard. We attribute this to the fact that the ray dynamics is integrable in this billiard and thus it is much more difficult to reach onset of chaos for the GPE in this billiard at moderate nonlinearity studied here. The situation is different for the stadium: at β = 2 only a few modes are populated, at β = 5 the number of modes is increased but still the numerical data for S are very far from the thermalization red curve (at least on the time scale reached in our numerical simulations). However, for β = 10, 20 we find that the numerical data at large times t > 15 follow the theoretical curve S(E) given by the Bose-Einstein thermalization. A small visible deviation from theory is still visible since the numerical points are systematically slightly below the theory curve. We attribute this to a finite computation time which apparently is not long enough to visit all regions of multi-configurational space with sufficiently large statistics. On the basis of obtained data we can conclude that the dynamical thermalization in the stadium sets in for β > β c ≈ 7 ≈ ∆. We also checked that the initial states, which represent a linear combination of a few eigenmodes, also follow the theoretical red curve in Fig. 3 at β > β c (e.g. two modes m = 10, 15 at β = 10).
Another way to check the thermalization predictions is to determine T from the numerical values of E, S which, according to the Bose-Einstein ansatz, give independent values T 1 (E) and T 2 (S). The average value T = (T 1 + T 2 )/2 is shown in Fig. 4 as a function of numerical values of E and S. We see that the numerical points are in a good agreement with the analytic curves (apart of small systematic displacement of numerical points discussed in the paragraph above). In a similar way we make a comparison between the theory and numerical data for dependencies µ(E) and µ(T ) shown in bottom panels of Fig. 3. Again we find a good agreement between the numerical data and the Bose-Einstein thermalization distribution (2).
The validity of the Bose-Einstein thermalization ansatz (2) leads to a striking paradox pointed already for nonlinear chains in [10]: formally the GPE in stadium gives a system of equations for nonlinear coupled oscillators (we have nonlinear coupling between linear oscillator modes of the stadium) with a moderate nonlinearity. The usual expectations of the statistical mechanics predict the energy equipartition between these modes [10,27]. If the number of modes is infinite then we should have ultra-violet catastrophe with probabilities ρ m ∼ 1/E m at high modes m and the global temperature approaching zero as T ∼ 1/m max for initial excitation with m ∼ 10 (here m max is the maximal mode number). This classical thermalization ansatz for m max = 50 is shown by a dashed curve in Fig. 3 and the data clearly show that it is very different from the numerical data which are close to the Bose-Einstein ansatz (we note that numerically the quantum Gibbs distribution [10] over quantum levels of stadium gives the results being rather close to those of (2) since at low temperatures and large m values both distributions are rather similar). Thus our data clearly show the emergence of the dynamical thermalization described by the Bose-Einstein distribution in a chaotic billiard for moderate nonlinearity β > β c ∼ ∆.
A more detailed check of the Bose-Einstein distribution requires a direct comparison of numerically obtained probabilities ρ m (m ) with the theoretical expression (2) for each initially excited mode m . We show such data in Fig. 5 for β = 2, 5, 10. It is clear that there is no thermalization at β = 2, 5 since a large fraction of probability remains at the initially populated state m = m . For β = 10 we see that the probability at initial state m drops significantly indicating emergence of dynamical thermalization. However, still the numerical probabilities at large m have larger values compared to those of the theory (2) shown in the bottom right panel of Fig. 5. We attribute this to the fact that our total computation time is not large enough to have good statistical data for average values of ρ m (m ) which require good averaging and long computation times. Such a problem had been visible in the numerical simulations with nonlinear chains [9,10] where the time of simulations have been by a few orders of magnitude larger than here. At the same time the extensive property of energy E and entropy S makes them self-averaging and more stable in respect to fluctuations thus allowing to compare them with the theory (2) at significantly shorter time scales.
Unfortunately the large scale simulations of the GPE for stadium are rather heavy and time consuming since they require transformations from coordinate to linear mode space and small integration time step with aliasing procedure to suppress numerical instabilities of high modes. It is possible that the numerical codes can be improved allowing to reach larger time scales but this requires further studies going beyond the scope of this work.
Finally we discuss a preparation of one or a few initial eigenstates considered above for the time evolution and dynamical thermalization. It is clear that the ground state of the billiard is relatively easy to prepare since it is the final state in a process of relaxation and also since it is compact in space being close to a coherent state of a harmonic trap. An excited state can be produced from the ground state applying a monochromatic driving (oscillation) of the billiard that creates an effective ac-potential V ac = f x cos(ωt) if one goes to the oscillating frame (see e.g. [28]). In a chaotic billiard dipole matrix elements have transitions between all energy eigenstates [28] and thus a resonant transition will populate one or a few states being close to the resonance E n ≈ E 0 +hω. We note that such a method demonstrated already its efficiency for excitation of high energy states for chaotic Rydberg atoms (see e.g. [29,30]).
Discussion. -Our studies of the GPE in the Bunimovich stadium billiard show that for a moderate nonlineariy parameter above a certain threshold β > β c ≈ ∆ the nonlinear Hamiltonian dynamics leads to emergence of dynamical thermalization over the linear billiard modes which is well described by the Bose-Einstein distribution. This distribution is strikingly different from the usually expected energy equipartition over modes [10,27] which would lead to a violet catastrophe with a significant probability transfer to higher and higher modes of the chaotic billiard. The established validity of the Bose-Einstein distribution, together with the previous studies of dynamical thermalization in nonlinear chains [9,10], leads to an unexpected conclusion about emergence of quantum distributions over linear energy modes in systems of coupled nonlinear oscillators at moderate nonlinearity. This result is drastically different from the standard energy equipartition picture expected for nonlinear dynamics of oscillator systems [1][2][3][4]27].
The described picture of "quantum" dynamical thermalization for the GPE in chaotic billiards requires a better understanding of nonlinear dynamics in systems with many degrees of freedom. It is known that slow chaos, like the Arnold diffusion [31][32][33], and an anomalous diffusion in disordered nonlinear chains (see e.g. [19,[34][35][36]) generate a number of features which still wait their deep under-standing. We hope that our results will stimulate further research in this field of fundamental aspects of nonlinear dynamics and thermalization onset in systems with large but finite number of degrees of freedom.
The modern progress in the cold atom experiments allows to investigate a dynamics of Bose gas and Bose-Einstein condensates [5,6] while the chaotic billiard for such atoms can be created by optical beams [15,16]. Thus we think that the model (1) can be realized with cold atom experiments. and middle panels show zoom marked in top panels (adding 15 trajectories in the integrable region to time t ≤ 10 4 ), one invariant curve of each panel is highlighted with red (gray) color inside stability island at middle panels. Bottom panels show dynamics in (x, y) plane with stable orbits from middle panels (same red color) and chaotic orbits (orange color also shown in top panels) up to times t = 10 2 ; dashed horizontal lines mark y = 0.5 used for the Poincaré sections.
Another promising possibility can be an experimental realization of a harmonic Sinai billiard, or Sinai oscillator. An example of such a billiard is described by a classical Hamiltonian H = (p x 2 + x 2 )/2 + (p y 2 + 2y 2 )/2 with a rigid disk of radius r d = 1 located at x = y = −1 (thus the ratio of frequencies in x and y is irrational). The harmonic potential can be realized by optical traps while the repulsive rigid disk can be created inside by an additional laser beam with such a frequency detuning that it acts as a repulsive potential for cold atoms. Examples of typical trajectories in such a billiard are shown in Fig. 6. In the same figure we also show the Poincaré sections constructed in (x, p x ) plane at y = 0.5, p y > 0 with energies E = 2 and E = 18, when the size of oscillations of atom is larger than r d . In this regime almost all phase space is chaotic (the domains of integrable dynamics are very small). Thus we think that the GPE in such a harmonic Sinai billiard will show all the effects of dynamical thermalization discussed above for a more convential case of the Bunimovich stadium. We expect that such a system can be more simple for experimental investigations. Also in such a billiard a coherent state of the harmonic potential can be created experimentally and can be used as an initial state with energy being close to the energies of linear eigenmodes of such a billiard. We expect that experimental investigations of the GPE in a harmonic Sinai billiard (or in the Bunimovich stadium) will allow to understand the fundamental aspects of dynamical thermalization.
We thank D.Guéry-Odelin for useful discussions of cold atom physics.