Nonthermal fixed points and solitons in a one-dimensional Bose gas

Single-particle momentum spectra for a dynamically evolving one-dimensional Bose gas are analysed in the semi-classical wave limit. Representing one of the simplest correlation functions these give information about possible universal scaling behaviour. Motivated by the previously discovered connection between (quasi-)topological field configurations, strong wave turbulence, and nonthermal fixed points of quantum field dynamics, soliton formation is studied with respect to the appearance of transient power-law spectra. A random-soliton model is developed to describe the spectra analytically, and the analogies and difference between the appearing power laws and those found in a field theory approach to strong wave turbulence are discussed. The results open a view on solitary wave dynamics from the point of view of critical phenomena far from thermal equilibrium and on a possibility to study this dynamics in experiment without the necessity of detecting solitons in situ.


I. INTRODUCTION
Many-body systems far away from thermal equilibrium can show a much wider range of characteristics than equilibrium systems. However, they buy this abundance for the price of transiency. Among the wealth of possible non-equilibrium many-body configurations most interesting candidates for theoretical and experimental study are those at which generic time-evolutions get stuck for an extraordinarily long time. This is in particular possible at critical points where universal properties like powerlaw behaviour of correlation functions appear, leading to slowing-down phenomena as infrared modes become dominant. Out of equilibrium, fluid turbulence is among the earliest described as a critical phenomenon [1][2][3][4].
In the realm of quantum physics, far-from-equilibrium many-body dynamics, from the formation of Bose-Einstein condensates in ultracold gases to quark-gluon plasmas produced in heavy-ion collisions and reheating after early-universe inflation, exhibits many interesting phenomena. In this context, increasing attention [5][6][7][8][9][10][11][12][13][14][15][16][17][18] is being given to wave turbulence phenomena [19][20][21]. Critical points far from equilibrium, so-called nonthermal fixed points, were proposed [6] and related to strong wave turbulence [7][8][9][10] and the formation of quasi-topological defects [15][16][17]. Such defects play an important role in superfluid turbulence, also referred to as quantum turbulence (QT), which has been the subject of extensive studies in the context of helium [22,23] and dilute Bose gases [24][25][26][27]. In contrast to eddies in classical fluids, vorticity in a superfluid is quantised [28,29], and the creation and annihilation processes of quantised vortices are distinctly different [22,23]. * Electronic address: t.gasenzer@uni-heidelberg.de Considering a single vortex or vortex line, geometry imposes a power-law shape of the angle-averaged flow velocity away from the vortex core. This spatial powerlaw dependence implies the power-law momentum spectrum predicted to occur at the nonthermal fixed point [15,17]. In this article we focus on solitary waves in onedimensional quantum gases and show that these exhibit scaling behaviour in analogy to the universal properties of vortex ensembles in two-and three-dimensional systems. Self similarity here is due to the absence of a scale in a soliton which has zero width and is randomly positioned in space. See Refs. [30,31] for related studies in relativistic field theory.
Solitons are non-dispersive wave solutions which can arise in many non-linear systems spanning a wide range from the earth's atmosphere [32] or water surface waves [33] to optics [34]. The characteristics of single or few solitons in ultracold Bose gases have been studied during the last decade regarding their movement and interaction in traps [35][36][37][38], their formation and creation [39][40][41][42] and their decay [43,44], see Ref. [45] for a review. Here, we propose to detect and characterise soliton excitations in ultracold gases by measuring the single-particle momentum spectra for such systems. We emphasise that these can give a strong indication for the presence of solitons even in cases where they cannot be observed in situ.
Superfluid turbulence plays an important role in the context of the kinetics of condensation and the development of long-range order in a dilute Bose gas. This, as well as turbulence in its acoustic excitations has been discussed in Refs. [46][47][48][49][50], and more recently in Refs. [15,17,18,51]. A possible observation of QT in ultracold atomic gases presently poses an exciting task for experiments [52][53][54]. We stress that the experimental study of superfluid turbulence and, more generally, of nonthermal fixed points in ultracold Bose gases is central for the understanding of the processes important during arXiv:1203.3651v1 [cond-mat.quant-gas] 16   the build-up of coherence and degeneracy in particular also in one-dimensional systems [55][56][57].
In the following we study the formation of soliton excitations in trapped one-dimensional Bose gases by means of simulations in the classical-wave limit of the underlying quantum field theory. We analyse the momentum spectra of the time-evolving system and discuss them with respect to nonthermal fixed points discussed in field theory. In Section II, we develop a model of independent, randomly positioned grey solitons, locally being solutions of the Gross-Pitaevskii equation describing a degenerate, one-dimensional dilute Bose gas. Single-particle momentum spectra are derived for such systems, both in a homogeneous system and within a trapping potential. A protocol for the formation of such soliton ensembles within a trapped gas is described in Section III, and the resulting states are analysed with respect to the predicted momentum spectra and power-law signatures of a nonthermal fixed points. Our conclusions are drawn in Section IV.

II. MOMENTUM SPECTRA OF SOLITON ENSEMBLES
In the classical-wave limit, a dilute ultracold Bose gas is well described by a positive definite Wigner phase-space distribution function W [φ, φ * ] for the complex field φ at each point in position or momentum space [58,59]. The dynamics of the gas is determined by the time evolution of this Wigner function according to the classical field equation (1) which is identical in form to the Gross-Pitaevskii equation (GPE) for the field expectation value. Here, m is the mass of the bosons, V an external trapping potential, and g 1D the coupling constant in one spatial dimension (1D).
The 1D Gross-Pitaevskii equation (1) possesses quasitopological soliton solutions which may travel with a fixed velocity but are non-dispersive, i.e., stationary in shape [45]. For positive coupling constant g 1D > 0 solitons are characterised by an exponentially localised density depression with respect to the surrounding bulk matter and a corresponding shift in the phase angle ϕ of the complex field φ = |φ| exp{iϕ}. Depending on the depth of this depression, the soliton is either called grey or, for maximum depression, black. On the background of a homogeneous bulk density n it is described by where x s (t) = x 0 + νt is the position of the soliton at time t. Here, ξ = [2mng 1D ] −1/2 is the healing length, and we express time in units of the inverse speed of sound c s = [ng 1D /m] 1/2 , i.e., t =t/c s , and drop the overbars. γ = 1/ √ 1 − ν 2 is the 'Lorentz factor' corresponding to the velocity v of the grey soliton in units of the speed of sound, ν = v/c s = |φ ν (vt, t)|/ √ n. Being related to the density minimum, ν is also termed the 'greyness' of the soliton, ranging between 0 (black soliton, |φ ν (vt, t)| = ν √ n = 0) and 1 (no soliton, |φ ν (vt, t)| = √ n).
In Sect. III we will study the formation and evolution of soliton excitations in trapped one-dimensional Bose gases by means of the GPE (1). An example of one run is shown in Fig. 1: A sudden initial quench of the coupling g 1D creates strong oscillations of the bulk gas in the trap and lets solitons form out of the short-wave-length collective oscillations. These solitons are seen as black lines surviving within the oscillating gas for a long time.
In this article we are mainly interested in the characterisation of the ensemble of solitons emerging in our simulations in terms of single-particle spectra in momentum space. Hence we refer to Sect. III for more details on the simulation protocol and results and start here with a discussion of the possible spectra.
The single-particle momentum spectrum is determined by evaluating ensemble averages different times during the period where the initial breathing oscillations are still present. Note the double logarithmic scale. At low momenta, the spectrum shows a plateau while at high momenta it falls of exponentially. At the point of time corresponding to the red curve, the cloud is further expanded such that solitons are more separated from each other which causes an intermediate power law dependence n(k) ∼ k −2 to appear, as indicated by the straight line. The plateau, the power-law and the exponential decay form the characteristic signature for solitons which we discuss in more detail in the following.
A. Random-soliton model: uniform gas To obtain an analytical understanding of the possible spectra in the context of nonthermal fixed points we discuss, in the following, the case of a dilute ensemble of well-separated solitons with random velocities and positions. The wave function of a single grey soliton in a homogeneous bulk background condensate is given in Eq. (2). Using this we write down an expression for a set of N s uncorrelated solitons with density minima at {x i , i = 1, ..., N s } and (dimensionless) velocities ν i , on a background of constant bulk density n: Note that due to the neglection of correlations this field does in general not represent a solution of the GPE in which the solitons remain non-dispersive.
We make use of the assumption that the ensemble is dilute, i.e., that the distance between each pair of neighbouring solitons is much larger than the healing length. This assumption, for grey solitons, is not valid as soon as two oppositely moving solitons encounter each other but for simplicity we will assume that these collisional configurations can be neglected in view of a majority of well-separated solitons. Since, for any i, 1 we can rewrite the spatial derivative of the field (4) as Here x i (t) = x i (0) − ν i t, β j = 2 arccos ν j , denotes the convolution over the spatial dependence on x, θ(x) the Heaviside function, and we have neglected an irrelevant overall phase. Note that the sign of β i indicates the direction of the propagation of the ith soliton. The term in square brackets in the second line of Eq. (5) is proportional to the spatial derivative of the field describing an ensemble of N s infinitely thin solitons (ξ → 0), at the positions {x i }, where n s is the number of solitons per unit length and the prefactor containing γ i takes into account that the phase jump by gives a sum of terms, each being proportional to a delta distribution, only one of these remains when evaluated at x i , which gives the term in square brackets in Eq. (5). We take the Fourier transform Here, P is the probability for finding a soliton with greyness ν = cos(β/2), and averaging over the random positions of all solitons other than those at x i and x j has been done in order to obtain the exponential decay of the coherence function. Combining Eq. (7) with the Fourier one derives the single-particle momentum distribution for a set of N s solitons defined by greyness and position, Here, the inverse volume δ(0) = L −1 appears as we first choose the φ * νi and φ νj in Eq. (7) different and take the identity limit only at the end. α is the average over all Assuming the dependence of γ i / sinh(πγ j kξ/ √ 2) on ν i to be negligible we obtain an approximate stationary distribution with a yet to be determined parameter γ. For black solitons (ν i ≡ 0) one obtains the exact expression For an ensemble of grey solitons of identical |ν i | ≡ ν, traveling with probabilities P into the positive x-direction and Q = 1 − P into the negative direction one finds relative distance on the average is much larger than their widths we chose the phase-space distribution to allow for a maximum greyness, |ν i | < |ν max | < 1, such that the diluteness criterium requires an approximate minimum box length of L = 4(1 − ν 2 max ) −1/2 N s . Fig. 3a shows the single-particle momentum spectrum n(k) on a double-logarithmic scale for an ensemble of 5 × 10 3 configurations with N s = 20 solitons each distributed according to a flat distribution across phasespace defined by the positions in the box and the maximum greyness |ν max | = 0.99. Solid (black) squares represent the results of the numerical ensemble average while the solid line corresponds to the analytical formula (10), with fitted parameters α = 0.7, γ = 1.05. Compare this to the analytical average α = 2/3. For comparison, we give the results for the same number of purely black solitons (red squares and line) as well as for a fixed greyness |ν| = 0.707 (blue squares and line), choosing an equal number of right-and left-movers. The comparison validates the approximate expressions (10)- (12) which exhibit scaling behaviour in the regime of momenta k ns γ −2 k k ξ γ −1 . In Fig. 4 we show the single-particle momentum spectrum for an ensemble of 5×10 3 configurations with N s = 20 solitons each distributed according to a phase-space distribution with an unequal weight for solitons with positive ν (right-movers) and negative ν (left-movers). We specifically restricted the greyness to the interval 0 = ν min ≤ ν ≤ ν max = 0.99 and besides that chose the same parameters as for Universal power-law behaviour in a many-body system far from equilibrium points to the appearance of turbulence phenomena. In the following, we summarise a few basic ideas of wave-turbulence theory for quantum gases and discuss the power-law spectra derived for soliton ensembles in this context.
A dilute, degenerate Bose gas is compressible such that collective sound-wave excitations can occur. This allows for so-called weak wave-turbulence which can occur in a regime where kinetic theory applies: The wave-kinetic equation for a Bose gas, which describes the evolution of the system under the collisional interactions between different wave modes, is known to have non-trivial stationary solutions. These solutions are nonthermal and exhibit a power-law dependence of mode occupations on momentum [19][20][21]. As in fluid turbulence, such solutions imply that energy flows from, e.g., small to large momentum scales. In between these momentum scales of the source and the sink, the energy passes through the so-called inertial interval, where the distribution over momenta is stationary and follows a power law with a universal exponent predicted by weak-wave-turbulence theory [19][20][21]. The stationary solution is metastable, i.e., requires a constant and equal flow in and out of the inertial interval. We note that weak wave turbulence does not bear, in general, vortical excitations. It is defined by scaling and stationary transport between different scales and hence can also occur in one spatial dimension where vortices are absent.
Applying the above to a degenerate Bose gas one faces, however, the problem that the description in terms of wave-kinetic equations breaks down in the infrared (IR) regime of long wavelengths where amplitudes, i.e., talking about a Bose gas, single-particle occupation numbers grow large and the description in terms of, e.g., elastic two-to-two collisions becomes unreliable [19]. As a consequence of this, so-called strong wave turbulence is expected to occur in the IR regime. Recent developments presented in Refs. [6][7][8][9][10]13] allow one to set up a unifying description of scaling, both in the ultraviolet (UV) regime, where the wave-kinetic and Quantum Boltzmann equations apply, and in the IR limit. Recall that at a critical point the IR modes dominate the system's behaviour. In this IR regime, new scaling laws were found by analysing non-perturbative Kadanoff-Baym dynamic equations, as in the UV, with respect to nonthermal stationary power-law solutions [6][7][8]10]. These solutions were termed nonthermal fixed points. Analogous predictions for dilute Bose gases were given in [9], proposing what was termed strong matter-wave turbulence in the regime of long-range excitations.
In [15,17], it was then shown that these nonthermal fixed points can also be understood, in two and three dimensions, in terms of vortex excitations of the superfluid: In the infrared limit of large wave numbers the incompressible, superfluid component of the gas dominates, and the predicted IR power laws appear due to the algebraic radial decay of the flow-velocity around the vortex cores. As a consequence, within a window in momentum space, which is limited by the inverse mean distance between different vortices and the inverse core size, the single-particle occupation number spectrum shows the power law predicted in [9].
Before we discuss this in more detail let us first turn back to the soliton spectra derived in the previous section. Assuming an equal number of solitons traveling with positive and negative velocities, P = Q = 1/2, i.e., assuming Im α = 0, the single-particle spectrum (12) is characterised by a maximum of two scales. Consider the case n s γ/( √ 2πξ). For momenta greater than the reduced soliton density but smaller than the reduced inverse healing length, k ns γ −2 k k ξ γ −1 , with k ns = 2n s , and k ξ = √ 2/(πξ), the momentum distribution exhibits a power-law behaviour, n(k) ∼ k −2 . This reflects, first, the random position of the kink-like phase jump across the center of the soliton, and second, that these momenta cannot resolve the spatial width of the kink. In other words, looking within a spatial window of size between γ 2 k −1 ξ and γk −1 ns , the appearance of a single sharp solitonic phase jump inside the window is observed in a random manner. Hence, for any of these window sizes, the system looks identical, it appears selfsimilar. This self-similarity is at the base of the scaling momentum distribution. Already for a single soliton, the momentum distribution does not know anything about the position of the kink (this information appears as a phase in the momentum-space Bose field which is irrelevant for n(k)), and thus the single-soliton distribution is self-similar, too. For black solitons, the self-similar scaling region is limited by the scales k ns and k ξ . Below k ns , the distribution is constant because too low wave numbers cannot resolve the kink-structure. This corresponds to the first-order coherence function decaying exponentially in space, with the decay scale set by the soliton distance 1/n s , Above k ξ , the momentum spectrum resolves the finite width of the soliton density dip which results in an exponential suppression of the mode occupations. We recall that also in equilibrium, at sufficiently high temperatures where quasiparticle mode occupations are large, n (qp) (k) = k B T /(c s k) 1/2, the first-order coherence function decays exponentially, g (1) where the scale is set by the coherence length r 0 = 2n/(mk B T ). Hence, the corresponding momentum spectrum has the same shape n eq (k) ∼ 2r −1 0 /(r −2 0 +k 2 ) as for a set of random thin solitons. We emphasise, however, that the transition scale k ns above which scaling ∼ k −2 sets in (see, e.g. (11)) can be made larger than for a thermal ensemble by increasing the soliton density n s above the inverse thermal coherence length 1/r 0 . This allows to identify non-equilibrium soliton vs. thermal scaling in experiment.
We finally compare the universal and non-universal aspects of the soliton momentum spectra found here with the corresponding spectra in d = 2 and 3 dimensions. As discussed in detail in Refs. [15,17], the universal exponent ζ = d + 2, d = 2, 3, found for the particle spectra n(k) ∼ k −ζ during the dynamical relaxation of an initially strongly quenched gas reflected the appearance of vortices in two and vortex lines in three dimensions. In two dimensions, this can be seen in an easy way looking at the flow velocity field v ∼ ∇ϕ = e θ /r at the distance r from the core of a singly quantised vortex, where e θ is the local tangential unit vector and ϕ the phase angle of the complex Bose field. The r-dependence implies a k −1 scaling of |v(k)| and thus a k −2 scaling of the kinetic energy E(k) ∼ k 2 n(k) ∼ v(k) 2 , i.e., n(k) ∼ k −4 [17,60]. Similar arguments lead to n(k) ∼ k −5 for the radial momentum distribution in the presence of a vortex line in three dimensions [17]. Extending these arguments, the scaling n(k) ∼ k −d−2 was shown to appear for ensembles of randomly positioned vortices/vortex lines in a range of momenta l −1 v k ξ −1 between the inverse of the intervortex distance l v and the inverse of the healing length ξ which is a measure for the core width.
As pointed out above, the power-law spectrum n(k) ∼ k −d−2 , in turn, had been predicted by use of nonperturbative field-theory methods in [6,9] where it resulted for a strong-wave-turbulence cascade in the IR, characterising the scaling behaviour at a nonthermal fixed point. This cascade was shown in [17] to be caused by particles being transported towards the IR where they build up high mode occupations and thus coherence in the sample, see also [46][47][48][49][50][51]. Note that in the picture of the evolving Bose field this momentum-space transport corresponds to the mutual annihilation of vortices and anti-vortices in the system which results in an increase of the intervortex distance and thus of the range over which phase coherence is established.
Having recalled all this, we note that there is a discrepancy between the predicted scaling ∼ k −d−2 , which was found consistent with the vortex picture for d = 2 and 3, and the scaling ∼ k −2 obtained here for the solitons in d = 1. To gain more insight into this issue we consider the spatial decay of the phase coherence for a system in two dimensions over distances considerably larger than the intervortex spacing. In analogy to the soliton ensembles discussed in Sect. II C, we consider randomly positioned and well-separated vortices. One finds that the decay of the phase coherence follows the same exponential law, where n v,1 is the one-dimensional uniform vortex density along the straight line through x and y. Fourier transforming (14) with respect to x−y results in a momentum spectrum n c (k) ∼ (4n 2 v,1 + k 2 ) −3/2 which scales as k −3 for momenta considerably larger than the inverse vortex distance n v,1 = 1/l v . Analogously, one finds n c (k) ∼ k −4 for vortex line tangles in three dimensions. We use the subscript 'c' to distinguish these spectra from the n(k) discussed above.
The apparent contradiction between the respective scalings of n(k) and n c (k) whose exponents differ by 1, is resolved by observing the following: While the exponential decay (14) is valid over large distances |x − y| l v , it is the algebraic dependence of the flow field v as a function of distance from the vortex core which matters on length scales below the inter-vortex distance, causing the steeper power law n(k) ∼ k −d−2 at momenta k 1/l v . Hence, we find an important qualitative difference between the physical properties underlying the universal scaling n(k) ∼ k −2 in one dimension and the scalings n(k) ∼ k −d−2 in d = 2 and 3 dimensions: While black solitons in one dimensions are at rest and no particle flow can occur, in higher dimensions, transverse flow circling around the vortex cores gives rise to an additional contribution to the kinetic energy. As far as this transverse flow dominates over a possible additional longitudinal flow component for which v · ∇v = 0, see [17], this comparison also holds when allowing for grey solitons, which are moving opposite to the (longitudinal) particle flow across the soliton dip.
In summary, there is a principal difference between the scalings of the momentum spectra in one-and in higher dimensions, giving rise to a deviation of the scaling ∼ k −2 in d = 1 from the field-theory prediction k −d−2 which, in turn, is valid in d = 2 and 3 dimensions. To recover this discrepancy within a field-theory approach to strong wave turbulence is beyond the scope of this article.

C. Random-soliton model: trapped gas
In our dynamical simulations we will consider soliton formation in a harmonic trapping potential rather than in a homogeneous system, see Sect. III. We therefore need to take into account, in the random-soliton model, the inhomogeneous bulk distribution of the gas.
Assuming a sufficiently shallow harmonic potential, l ho /ξ = (mω ho ξ 2 ) −1/2 = (2g 1D nω ho ) −1/2 n −1 s we can describe the Bose field in local-density approximation with respect to a bulk density distribution given in Thomas-Fermi approximation, n TF (x) n 0 [1 − (x/R) 2 ], with R = √ 2c s /(ω ho ξ) = 2g 1D n/ω ho being the Thomas-Fermi radius in units of ξ. We take the maximum density n 0 large enough to ensure k ξ k ns for solitons not too close to the edge of the cloud. In such a bulk density distribution, single solitons oscillate harmonically between classical turning points where the solitons "touch ground", i.e., momentarily turn black [44]. Their oscillation frequency is by a factor of ∼ 1/ √ 2 smaller than the trap frequency, ω s ω ho / √ 2. In leading order in ∼ẍ s /ẋ s , the field of a single soliton can locally be written in the simple form given in Eq. (2), with ν, and thus γ replaced by local quantities, [44]. ν s,max is the maximum greyness the soliton acquired in the center of the trap, and x s,max = ν s,max R is the distance of the soliton's turning point from the trap center. Only solitons whose velocity does not exceed the Landau critical velocity, i.e., for which ν s,max ≤ 1, can oscillate in the trap for more than a quarter of the period T s = 2π/ω s . This limits the maximum greyness at a distance x from the trap center to a range between 0 and ν max (x) = 1 − (x/R) 2 . At a given time t, we assume a particular set {x i (t)} of N s well-separated solitons across the trapped gas. The single-particle momentum spectrum corresponding to an ensemble of such sets depends on the distribution of the solitons over the greyness for each position in the trap. This distribution is best visualised in phase space which is parametrised by the (x, v), or, equivalently and in dimensionless form, by (x/R, ν), with both, x/R and ν ranging between −1 and 1. In this space, the trajectory of a single soliton is a circle with radius ν s,max which is traced out with constant angular velocity ω s . Hence, a stationary distribution of the N s solitons is given by a circularly symmetric distribution in phase space, i.e., a distribution over the different possible maximum greynesses ν s,max or turning points x s,max . The simplest assumption would be that of a uniform distribution of the solitons in phase space, which amounts to a uniform distribution over the different possible ν at each distance x from the trap center and an integrated soliton density distribution n s (x) ∝ 1 − (x/R) 2 .
Following the above considerations we can obtain approximate expressions for the momentum spectrum. For instance, for a uniform density P [x, ν] ≡ n s,0 /R s of solitons within a radiusR s = R s /R in phase space, i.e., for all (x, ν) with √x 2 + ν 2 ≤R s ≤ 1, with n s,0 = N s /(R s π), the first-order coherence function for thin solitons becomes = n TF (x)n TF (ȳ) exp{−2n s,0 x y dz α(z)} (16) with a local average dephasing of Using this, the integral over α to linear order inx reads . This approximation is best forR s → 1, in which limit we obtain Analogously, one calculates the exponential dephasing factor for more complicated soliton phase-space distributions. Taking the above results together one derives the momentum distribution in local-density-approximation as the convolution of the spectrum for a homogeneous distribution of thin solitons with the Fourier transform of the bulk density, multiplied with the momentum spectrum of a single soliton: where the parameter γ is to be determined and denotes the convolution with respect to k.
In the case that single soliton distributions contributing to the ensemble are not uniform throughout phase space, the ensemble-averaged n(k) would rather be a sum of (k → −k)-asymmetric distributions such that on the average the momentum distribution can have local maxima at finite |k| > 0.
To study the quality of the above analytic expressions we construct ensembles of phase-space distributions of spatially well-separated solitons inside a harmonic trap and compute the ensemble average (3). For this we multiply N s single-soliton solutions (2) with positions x i and greyness ν i chosen according to a given phase-space probability distribution and ensure that their relative distance on the average is much larger than their widths.  In the following we study the formation of soliton ensembles by use of semiclassical simulations, with Gaussian noise for the initial field modes. Moments of the phase-space probability distribution at a later time are determined by sampling the initial distribution, propagating each realisation according to the GPE, and averaging over many such trajectories [58,59]. At the initial time, we take the gas to be noninteracting and thermalised and impose an interaction quench. To allow the emerging collective excitations to form solitons at a desired density we furthermore apply evaporative cooling by opening the trapping potential at the edges in a controlled fashion. During the first, cooling period, t ≤ t c , the potential is given by the inverted } with its maximum being ramped down by sweeping U (t) = U 0 + (U c − U 0 )t/t c linearly in time from U 0 to U c . At the same time, highly energetic particles near the edge of the potential are removed by adding a loss term iΓ(x, t)/2 = iΓ ∞ [V (x, t)/U (t)] r /2 to the trapping potential. Thereafter, during the interval t c ≤ t ≤ t max the loss is switched off and the potential is ramped up again to harmonic shape accross the extension of the gas, . We choose r = 10, U 0 = 2.75, Γ ∞ = 0.1 U 0 , U c = U 0 /3, and U max = 10 U 0 . The times t c and t max vary and are given in the following. This protocol corresponds to the one used in Ref. [61]. Different cooling schemes have been used in experiments, see, e.g., [55,57], but as we are primarily interested in the one-dimensional dynamics, we here restrict ourselves to purely 1D calculations.
For the simulations, we map the system onto a grid of N = 1024 lattice sites with a lattice constant a s . If not stated otherwise, quantities are given in grid units based on a s , and the parameters chosen are a dimensionless coupling constant g 1D = 2ma s g 1D = 7.3 × 10 −3 , a cooling time t c = t c /(2ma 2 ) = 9.1 × 10 3 , and a harmonic oscillator length l ho = l ho /a = 8.5. Lattice momenta are k n = 2 sin(πn/N ), n ∈ {−N/2, ..., N/2 − 1}. We drop overbars in the following.
Three stages of the induced dynamical evolution can be observed, see also The smallest time scale is the oscillation period in the trap T ho ≈ 230, which leads to an initial collective oscillation with period T oscillation ≈ 300 (cf. Fig. 1). The collective breathing motion dies out after τ ≈ 2000. The oscillation period of a soliton in the Thomas Fermi bulk is T s ≈ √ 2 T ho ≈ 320. The longest time scale in our setup is the cooling time t c = 9100. Comparing these time scales to the total time of the simulation, at the end of which solitons are still present, we see that the solitons are quasi-stationary in the system. They emerge soon after the initial quench and remain throughout the whole evolution while thermalisation of the high-momentum modes is proceeding as we will see in the following.

B. Number of solitons after cooling ends
In order to study the statistics of the solitons emerging during the evolution we have set up an efficient tracking algorithm which identifies the trajectories of the solitons oscillating in the gas. The algorithm scans the wave function for density minima coinciding with a phase jump around them. Fig. 7 shows the evolution of the mean number of solitons, for an ensemble of 200 runs. The three stages described above can be identified. The strong initial oscillations give an oscillating number of solitons until t 10 3 . The number of solitons decreases while the gas is evaporatively cooled. After the end of the cooling, at t = t c the decay is considerably slowed down and the number of solitons remains largely stable. Three different cooling times and two ramp speeds are shown, with t c = 3.9 × 10 3 (red), t c = 4.55 × 10 3 (blue), and t c = 9.1 × 10 3 (black), where the same speed is chosen to obtain the blue and black data.

Kibble and
Zurek have predicted that the number of defects created in the near-adiabatic crossing of a phase transition scales with the crossing rate according to a power law which depends on the universal properties of the transition [62,63]. This was studied numerically in [61] using the cooling protocol described above. While the interacting gas was chosen to be in thermal equilibrium initially, with a temperature well above the critical point, we start our simulations, motivated by earlier work on vortex dynamics [15,17], with an interaction quench driving the system strongly out of equilibrium. To compare the dynamics induced in this way with the results of [61] we show, in Fig. 8, the dependence of the number of solitons created on the cooling ramp time t c . We find that, within the error bars which indicate the variance over 200 runs, the data is rather fitted by an exponential dependence N s (t c ) = f 0 exp{−γt c } than by a power law N s (t c ) = g 0 /t c as predicted in Ref. [39]. We emphasise however that in our system, solitons mainly form during the initial stage following the interaction quench. . The function f0 exp{−γtc}, with f0 = 9.7 and γ = 9.77 · 10 −5 , was fitted with χ 2 = 0.006 (black line). The function g0t −1 c , with g0 = 1.55 · 10 5 , was fitted with χ 2 = 1.28 (red dashed line). Hence, other than for the Kibble-Zurek scheme of Ref. [61], cooling after the initial quench results in an exponential dependence of Ns(tc) on tc.
C. Time evolution of single-particle spectra We finally discuss the relaxation dynamics with respect to the evolution of the respective single-particle momentum spectra (3). The initial state chosen in the simulations is given by a thermal canonical ensemble of distributions over the single-particle eigenstates of the trap. In Fig. 9 we show the momentum spectrum at time t = 5·10 3 . Solitons have formed at high density such that the scales k ns 0.15 and k ξ 0.2 are close together, as indicated in the graph. The solid line represents a fit of the analytical model spectrum (12). A Rayleigh-Jeans tail is absent as the cooling is still on. Due to the proximity of k ns and k ξ no k −2 power law is seen in between the low-energy plateau and the high-energy exponential fall-off.
In Sect. II, Fig. 2 we showed the spectrum for a wider trap with l ho = 17 which allows the solitons to be diluted more across the trap and results in k ns 0.15 and k ξ 0.65. This allows for a k −2 scaling to appear clearly, indicating a self-similar random distribution of the solitons. For t > t c , the gas enters the final stage: The number of particles and the energy are now conserved and the healing length is 'frozen out'. Fig. 10 shows the development of n(k, t) from t t c to late times. Once the cooling and thus the removal of particles with high energy is terminated, a transport process from low to high momenta starts and thermalisation takes place. The influence of solitons is still present with two effects: First, there are still solitons in the gas for the times displayed in Fig. 10 which contribute with their spectral profile to the total spectrum and broaden the plateau at low momenta up to k ξ 0.1. Second, the momentum distribution starts to oscillate between situations with a stronger weight on the positive and on the negative side. Fig. 11 shows how k ns σ −1 k ξ Occupation number n(k) FIG. 9: Momentum spectrum n(k, t) at time t = 5 · 10 3 with statistical errors (ensemble average over 100 runs). k and t are in grid units as defined in Sect. III A where also all other simulation parameters are given. Solitons have formed at high density such that the scales kn s and k ξ are close together, as indicated in the graph. Solid line: fit of the analytical model spectrum (19), with ns,0 = 0.076, γ = 1, σ − 1 = 0.036 where Gaussian of width σ was used to describe the bulk distribution in position space. A Rayleigh-Jeans tail is absent as the cooling is still on. Due to the proximity of kn s and k ξ no k −2 power law is seen in between the low-energy plateau and the high-energy exponential fall-off. π/R k ξ Occupation number n(k) t = 9 × 10 3 t = 10 5 t = 10 6 k −2 FIG. 10: The spectrum n(k) is shown for three different times from t ≈ tc = 9 × 10 3 to t = 10 6 , for an ensemble of 5 × 10 3 (10 3 for t = 10 6 ) realisations. k and t are in grid units as defined in Sect. III A where also all other simulation parameters are given. Momentum scales defined by the inverse Thomas-Fermi radius R as well as by the healing length ξ are marked by dashed lines. A power-law dependence ∼ k −2 of the thermal tail is indicated by the black solid line.
the remaining solitons, which oscillate in the trap, influence the spectrum with their own momentum appearing as the oscillating maxima in the spectrum. The oscillatory pattern appears despite a large number of runs contributing to the statistical ensemble. A small number of runs with a few oscillating solitons dominate the averaged spectrum.

IV. CONCLUSIONS
We have studied the formation of dark solitary waves in one-dimensional Bose-Einstein condensates as well as their relaxation dynamics towards equilibrium. The corresponding single-particle momentum spectra were predicted in the framework of an instantaneous model of well-separated grey solitons, whose width is considerably smaller than their mutual distances in the bulk. For comparison with these predictions, semiclassical simulations of the relaxation dynamics of one-dimensional Bose gases after an initial interaction quench and a cooling period were used to determine the respective spectra numerically. The so found spectra compared well with the analytical predictions, giving insight into the many-body dynamics from the point of view of universal properties and critical physics far from equilibrium. We emphasise that the particular protocol used to produce the solitons is irrelevant in that the properties characterising the fixed point do not depend on how it is reached.
We have discussed the power-law behaviour of the momentum spectra, which appears in a range of momenta between the inverse of the inter-soliton distance and the inverse healing length, with regard to the uni-versal scaling laws predicted in non-perturbative fieldtheory approaches to strong wave turbulence. In the onedimensional case studied here, the derived power-law exponent ζ = 2 differs by 1 from the exponent ζ = d+2 predicted for a system in d spatial dimensions and previously recovered in the scaling due to vortex excitations in a d = 2 and 3-dimensional Bose gas. We trace this discrepancy back to the different flow patterns possible in d > 1 versus d = 1 dimensions: While in the one-dimensional gas, particle flow cannot choose its orientation except for a sign, vortical excitations in two and three dimensions are characterised by flow circling around the vortex cores. This flow is transverse, i.e., it changes its strength perpendicular to its direction. It dominates, at sufficiently low energies and momenta, over any additional longitudinal flow caused by compressible sound-wave excitations. Moreover, for geometric reasons, the transverse-flow velocity field changes algebraically in space and causes the single-particle momentum spectra to scale as predicted in strong-wave-turbulence theory.
We point out that the single-particle momentum spectra discussed here could be used in experiment to study solitary-wave dynamics in one-dimensional Bose gases without the necessity to detect solitons in situ. Studying in this way universal properties during the relaxation dynamics from a non-equilibrium initial state or under a constant driving force opens a new access to strong wave turbulence and nonthermal fixed points.