One-Dimensional Fuzzy Dark Matter Models: Structure Growth and Asymptotic Dynamics
Tim Zimmermann, Nico Schwersenz, Massimo Pietroni, Sandro Wimberger
aa r X i v : . [ a s t r o - ph . C O ] F e b One-Dimensional Fuzzy Dark Matter Models:Structure Growth and Asymptotic Dynamics
Tim Zimmermann, ∗ Nico Schwersenz, † Massimo Pietroni,
2, 3, ‡ and Sandro Wimberger
2, 3, § Institut f¨ur Theoretische Physik, Philosophenweg 16, 69120 Heidelberg, Germany Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Universit´a di Parma,Campus Universitario, Parco Area delle Scienze n. 7/a, 43124 Parma, Italy INFN, Sezione di Milano Bicocca, Gruppo Collegato di Parma, 43124 Parma, Italy (Dated: March 1, 2021)This paper investigates the feasibility of simulating Fuzzy Dark Matter (FDM) with a reducednumber of spatial dimensions. Our aim is to set up a realistic, yet numerically inexpensive, toymodel in (1 + 1)-dimensional space time, that — under well controlled system conditions — iscapable of realizing important aspects of the full-fledged (3 + 1)-FDM phenomenology by means ofone-dimensional analogues. Based on the coupled, nonlinear and nonlocal (3+1)-Sch¨odinger-Poissonequation under periodic boundary conditions, we derive two distinct one-dimensional models thatdiffer in their transversal matter distribution and consequently in their nonlocal interaction alongthe single dimension of interest. We show that these discrepancies change the relaxation processof initial states as well as the asymptotic, i.e., thermalized and virialized, equilibrium state. Ourinvestigation includes the dynamical evolution of artificial initial conditions for non-expanding space,as well as cosmological initial conditions in expanding space. The findings of this work are relevant forthe interpretation of numerical simulation data modelling nonrelativistic fuzzy cold dark matter inreduced dimensions, in the quest for testing such models and for possible laboratory implementationsof them.
I. INTRODUCTION
Nonlinear Schr¨odinger equations are ubiquitous inphysics. Let us just think of interacting many-body prob-lems in nonrelativistic condensed-matter theory that arereduced to a mean-field approximation that usually endsup in a nonlinear effective Schr¨odinger equation, see e.g.[1–4]. A specific form of a such a nonlinear equation is theSchr¨odinger-Poisson, also known as Schr¨odinger-Newtonequation, describing a scalar massive quantum particlein its own gravitations field. It finds many applications,e.g., in nonlinear optics [5–8], decoherence theory [9], and,of course, in cosmology [10], whenever a nonrelativisticdescription of the particle suffices. Ref. [11] presents arecent review on the subject.Here, we are interested in the Schr¨odinger-Poisson (SP)equation as an alternative dark matter model to the es-tablished cold dark matter (CDM) paradigm. Hu et al. [12] dubbed the matter so described fuzzy cold dark mat-ter, or simply fuzzy dark matter (FDM), a notion thatwe will follow throughout the paper.The most compelling feature of FDM obeying SP infour dimensional (3+1) space-time is its distinct behavioron large and small spatial scales: Assuming a dark mat-ter particle mass of m ≈ − eV — canonical for FDM— cosmic structure growth under FDM is in accordancewith CDM on super-galactic scales, e.g., identical matterpower evolution or halo densities [13], sub-galactic scales ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] are influenced by quantum mechanical effects. In partic-ular, one expects the existence of a solitonic state, witha flat, high-density core region, that is roughly speak-ing obtained from a balance of quantum pressure andgravity, [12]. Structure smaller than the solitonic coreis suppressed, or smoothed out, by the uncertainty prin-ciple. Thus, FDM may provide a natural solution tothe small scale crisis of CDM [14], and in particular itscusp-core problem, [15–17], without the need of addingsophisticated baryonic feedback processes. From a fun-damental physics point of view, SP simulations couldpredict the mass scales of the bosonic particles possiblyconstituting dark matter, of course in comparison withand constrained by experimental observation data [18–21]. We note in passing that SP simulations may also beinterpreted as an alternative sampling approach of theCDM phase space evolution that may be controlled bythe phase space resolution ~ /m , see e.g. [22–24].Numerical considerations of the (3 + 1)-SP equation,see e.g. [13, 17, 25–28], identify the soliton state asa dynamical attractor in the time evolution of FDM.More precisely, overdense regions collapse under theirself-gravity and during this process radiate away excessmatter. The result of this process, sometimes dubbedgravitational cooling, [29], is a relaxed quantum matterdistribution in which a solitonic core is immersed in a ’seaof fluctuations’. The latter follows a power-law decay,with a radial profile, known as NFW profile, predictedin corresponding many-body simulations of CDM [30] toscale as ρ ∝ r − at large radii.This paper addresses the question whether the time-asymptotic behavior just described, or at least a similarscenario thereof, were recovered if only one spatial degreeof freedom is available. In other words, can we derivea one-dimensional, yet sufficiently realistic, toy modelthat realizes one-dimensional analogues of the (3 + 1)-FDM structure formation, and in particular the men-tioned solitonic core and its role as dynamical dynami-cal attractors? This includes, both the evolution to andthe possible reaching of the asymptotic state, hence therelaxation process itself as well as its final product. Arethe relaxation mechanisms in one-dimensional FDM thesame as for 3D FDM, i.e., is there a one-dimensional ana-logue of gravitational cooling? What are suitable astro-physical, quantum mechanical or statistical measures tojudge if the asymptotic state has been reached? Thesequestions have several motivations. Firstly, it is a pri-ori not clear whether a model with reduced dimensionwill lead to the same evolution as the full-fledged threedimensional one, simply because assumptions on the mat-ter configuration in the transversal dimensions must bemade and, in higher dimensions, interdimensional cou-pling and redistribution of mass is, in principle, possibledue to the non-separable nonlinearity [31]. Secondly, it isneedless to say that models with reduced dimensions lendthemselves to much more detailed investigations even onrelatively long spatial and temporal scales. Reliable high-precision numerical simulations are simply more efficientin one than in higher dimensions. Consequently, largerpatches of the system’s parameter space can be coveredwith manageable resources. The reduced numerical com-plexity also allows to generate large ensembles of sim-ulations once a, in some sense optimal, parameter setwas identified, thereby reducing statistical uncertaintieswithout the need to invoke assumptions like ergodicity.Obviously, both aspects, i.e., parameter space coverageand statistical reliability, are vital for comparison withobservational data.The way of the reduction of dimensionality will cer-tainly be important. Anticipating our results, we will seethat the relaxation process in the evolution as well asin the finally reached states do depend strongly on howmatter is organized in the transversal dimensions. Theusual approach of assuming a uniform matter distribu-tion in the transverse dimension, and hence in practicejust forgetting about the other degrees of freedom, doesnot yield a correspondence with the (3 + 1)-FDM predic-tions. This motivated our second approach in which weconfine the transverse degrees of freedom. Then, mostof these predictions are actually seen in the reduced one-dimensional model as well. Our findings are based on arelative simple but highly flexible, efficient and precise nu-merical method to integrate both dimensionality-reducedSP equations. This allows us access to long temporal evo-lutions in a constant as well as in an expanding universe.The paper is organized as follows: the next Sec. re-views the theory of FDM based on the SP equation. Wediscuss therein the importance of the reduction to (1 + 1)(space+time) dimensions, and we compare, in particular,the two ways of doing this: first, we assume a uniformextension in the transverse directions, then secondly acigar-shaped transversal confinement. Sec. III gives de- tails on our numerical methods and its properties andquality. Sec. IV reports our central results and com-pares the outcome of different numerical simulations ofthe two versions of reduction to (1 + 1). The last Sec. Vconcludes the paper with a short outlook on open ques-tions and future work. II. THEORETICAL BACKGROUND
Point of departure for our discussion is the dynamicsof Fuzzy Dark Matter (FDM) in (3 + 1) dimensions. Tothis end, we introduce its governing equation, commenton its origin and mention different interpretations of theFDM model. After a quick review of the hydrodynamicformalism in the linear evolution regime, the remainder ofthis section is devoted to the problem of dimensional re-duction. More precisely, two competing one dimensionalFDM representations are derived and their physical dis-crepancies and similarities are highlighted. This includesan in-depth discussion on the preparation and propertiesof the (1+1)-FDM ground state. We close this theory sec-tion by commenting on possible relaxation mechanismsand introduce suitable metrics for (1 + 1)-FDM to quan-tify if an equilibrated system state is reached.
A. Fuzzy Dark Matter in (3 + 1)
Dimensions
1. Governing Equations
For the sake of simplicity, envision a universe with neg-ligible contribution of baryonic matter and radiation to-wards the cosmic energy budget. Furthermore, let darkenergy be given by means of a time independent energydensity ρ Λ c . Under these assumptions only dark matterdynamics plays a nontrivial role.Starting from a massive, minimally coupled, scalarfield, one derives, e.g. [32], the (3 + 1)-Schr¨odinger-Poisson (SP) equation as governing equation of (3 + 1)-FDM applicable in the nonrelativistic limit c → ∞ : i ~ ∂ t Ψ = (cid:20) − ~ ma △ + m Φ (cid:21) Ψ , △ Φ = 4 πGa (cid:0) | Ψ | − ρ m (cid:1) , x ∈ Ω , (1)with Ω = Ω × Ω × Ω ⊂ R and Ω i = [0 , L i ). Here, Ψdenotes the nonrelativistic FDM field coupled to its owngravitational potential Φ. Fields are evaluated at comov-ing position x = ( x , x , x ) ⊺ and cosmic time t . Densi-ties are measured with respect to a comoving volume, sothat ρ m coincides with the present day total matter den-sity. m denotes the FDM particle mass. In accordancewith our initial assumptions on the composition of thecosmic energy budget, the scale factor a ( t ) obeys the flatspace, radiation free Friedmann equation: (cid:18) ˙ aa (cid:19) = 8 πG (cid:0) ρ m a − + ρ Λ (cid:1) .G is Newton’s constant.Analyzing the asymptotic behavior of the energy mo-mentum tensor associated with Ψ suggests the identifica-tion: | Ψ( x , t ) | = ρ m ( x , t ) ≡ ρ m + δρ ( x , t ) , (2)i.e. the scalar field encapsulated both the time-independent background density ρ m and deviations fromit, δρ ( x , t ). Since density deviations vanish upon averag-ing over Ω, Eq. (2) translates directly into a normaliza-tion condition on Ψ: ρ m = 1 L L L Z Ω d x | Ψ( x , t ) | = const . (3)We emphasize Eq. (3) is a physically relevant constraintas Eq. (1) is nonlinear . Therefore changing the normal-ization of Ψ will lead to a different time evolution.Eq. (1) still requires suitable boundary conditionsooup The natural choice for modelling an infinite systemis to impose the following periodic boundary conditions: ∂ mx i Ψ(0 , x j ) = ∂ mx i Ψ( L i , x j ) ,∂ mx i Φ(0 , x j ) = ∂ mx i Φ( L i , x j ) , m ∈ { , } ,i ∈ { , , } ,j ∈ { , , } \ { i } . (4)
2. Interpretations of (3 + 1) -Schr¨odinger-Poisson
It is instructive to shed some light onto the physi-cal character of the non-relativistic scalar Ψ. AlthoughEq. (1) has the mathematical structure of Schr¨odinger’sequation, there is a priori nothing quantum mechanicalabout the problem. Hence, the least arcane way to inter-pret Eq. (1) is in a literal sense, i.e. as the Euler-Lagrangeequation of a classical Lagrangian and ~ as a constantwith dimensions of an action but with a numerical valuenot constrained to Planck’s constant.That said, there is significant value in finding (formal)correspondences between (3 + 1)-SP and other, poten-tially non-cosmological, theories as it enlarges the num-ber of available tools with which FDM can be analyzed. a. Cosmic Bose-Einstein Condensate In fact,Eq. (1) can also be identified with the evolution equationof a self-gravitating Bose-Einstein condensate withnegligible local self-interaction and Ψ as the condensate wave function . The authors of [33] substantiates thisclaim by comparing the critical temperature of anultralight boson gas undergoing pair production withthe cosmic microwave background temperature.A quantum mechanical derivation for Eq. (1) couldthen depart from a second quantized many-body Hamilto-nian which is subsequently reduced to an effective Hamil-tonian for the order parameter, i.e. the condensate wave function Ψ, [4]. The result would be again Eq. (1) with ~ as Planck’s constant.We return to this quantum mechanical interpretationin Sec. III when the time integration is formulated asapproximation to time evolution operator or in Sec. II C 1when the quantum virial theorem is used to analyze thelong term FDM dynamics. b. Smoothed CDM Dynamics On the other hand, ifwe accept Eq. (1) as an abstract evolutionary problemand forget momentarily about its interpretation as al-ternative dark matter model, the dynamics of Ψ can beassociated with a smoothed version of the Vlasov-Poissonequation (VP) — the phase space description of CDM,[22–24].More precisely, if f V ( x , p ) denotes the solution to VPand f W ( x , p ) the Wigner phase space distribution, con-structed from Ψ via: f W ( x , p ) = Z d x ′ Ψ (cid:18) x − x ′ (cid:19) Ψ ∗ (cid:18) x + x ′ (cid:19) e i ~ p · x ′ , then the evolution of the smoothed, or convolved, distri-butions,¯ f V/W ( x , p ) = exp (cid:18) − x σ x − σ x ~ p (cid:19) ∗ f V/W , (5)obeys, [23]: ∂ t ( ¯ f W − ¯ f V ) = ~ ∂ x i ∂ x j ∇ x ¯ V ∂ p i ∂ p j ∇ p ¯ f W + O ( ~ , ~ σ x ) . Here, ~ is not Planck’s constant and ~ /m acts as a modelparameter that sets the maximum phase space resolution.Then, σ x is a artificial smoothing scale in comoving posi-tion space.In that sense, Eq. (1) can be understood as an alter-native sampling of the CDM distribution compared tothe N -body approach: Instead of following the evolu-tion of N test particles sampling f V , we coarse grainthe phase space distribution directly and use Ψ as a dy-namical proxy for its evolution. See Sec. II C 2 for anapplication of this interpretation.
3. Dynamics in the Linear Regime
The behavior of FDM in the linear growth regime iswell-established in the literature, see e.g. [32–34]. Weare therefore brief and only focus on aspects relevant forour interpretation of the numerical simulations later onin Sec. IV A.In short, Madelung’s ansatz, [35], of decomposingthe wave function into Ψ( x , t ) = p ρ ( x , t ) exp (cid:16) i S ( x ,t ) ~ (cid:17) ,turns Eq. (1) into the irrotational Euler-Poisson equationincluding an additional (quantum) pressure term: v ≡ ma ∇ S = ~ ma | Ψ | Im (Ψ ∗ ∇ Ψ) , (6a) ∂ t ρ + 1 a ∇ · ( ρ v ) = 0 , (6b) △ Φ = 4 πGa ( ρ − ρ m ) , (6c) ∂ t v + 1 a ( v · ∇ ) v + H v = − a ∇ Φ + ~ m a ∇ (cid:18) △√ ρ √ ρ (cid:19) , (6d)Linearizing Eq. (6b)-(6d) up to first order in v and ρ ( x , t ) = ρ m (1 + δ ( x , t )) yields a damped oscillator equa-tion for the density contrast δ in the reciprocal domain:0 = ∂ t ˆ δ + 2 H ( t ) ∂ t ˆ δ − πGρ m a " − (cid:18) kk J ( a ) (cid:19) ˆ δ (7)with k J ( a ) = (cid:18) πGρ m m a ~ (cid:19) , (8)and ˆ δ = ˆ δ ( k, t ). As for CDM, no mode coupling oc-curs and all perturbations evolve independently underFDM evolution. However, in contrast to CDM, largeand small scale modes behave differently. Most notably,Eq. (8) defines a time-dependent critical length scale, λ J ( a ) = 2 π/k J ( a ), — the Jeans scale— below which thequantum pressure counteracts gravity so that density per-turbations do not collapse under their self-gravity.The linear Jeans scale can also be understood as aconsequence of Heisenberg’s uncertainty principle, [12], mσ r σ v ≃ ~ . In hydrodynamic terms σ v may be inter-preted as a velocity dispersion and a simple way to esti-mate it in the linear regime is to follow a particle trappedinside a gravitational well of a matter distribution withdensity ρ m : σ v ≈ rt dyn ≈ ax p Gρ m a − , with t dyn as dynamical time scale estimate. Thus: σ x ≃ ~ mx √ Gρ m a . Setting x = σ x yields Eq. (8) up to a numerical constantof O (1).The interpretation then is that the source of the quan-tum pressure is Heisenberg’s uncertainty principle whichinduces an increasing velocity dispersion in the FDM con-densate once particles are confined to a space region thatis comparable to σ x .Although, the Jeans length of Eq. (8) is a purely lin-ear concept, the uncertainty principle is not. One should therefore expect a distinct signature in the matter powerspectrum, even deeply in the nonlinear regime. We re-turn to this in Sec. IV A 3.Once a perturbation mode ˆ δ ( k, t ) leaves the linearregime, Eq. (7) is not applicable anymore. One thenexpects mode coupling to take place and thus a redistri-bution of power across all scales that participate in thenonlinear evolution. The implications of this effect willbe analyzed in Sec. IV A 4. B. Fuzzy Dark Matter in (1 + 1)
Dimensions
Let us now turn the attention to one dimensional ap-proximations to (3 + 1)-SP and how such models maybe deduced from Eq. (1). A convenient starting point forthe dimension reduction procedure is to subsume Eq. (3)-(4) into a single dimensionless, nonlinear Schr¨odingerequation (NLSE). This simplifies the discussion and isachieved by (i) absorbing the solution of Poisson’s equa-tion by means of a convolution integral and (ii) adoptingdimensionless quantities. One arrives at: i ~ ∂ t Ψ = (cid:20) − △ + a ( t ) (cid:16) G ppp △ ∗ | Ψ | (cid:17)(cid:21) Ψ x ∈ Ω , (9)where we defined: x ′ ≡ (cid:16) m ~ (cid:17) (cid:20) H Ω m (cid:21) x , d t ′ ≡ a (cid:20) H Ω m (cid:21) d t , Ψ ′ ( x ′ , t ′ ) ≡ Ψ( x ′ , t ′ ) √ ρ m , V ( x ′ , t ′ ) ≡ a m ~ (cid:20) H Ω m (cid:21) − Φand dropped all primes subsequently. The nonlinear po-tential is then given by: G ppp △ ∗ | Ψ | = Z Ω d x ′ G ppp △ ( x − x ′ ) | Ψ( x ′ ) | , and G ppp △ denotes the Green’s function of the d = 3 Pois-son equation augmented with periodic boundary condi-tions in all three dimensions. Eq. (9) takes the form ofa non-autonomous, i.e. explicitly time-dependent, NLSEwith a long range (i.e. non-local) interaction kernel.We stress G ppp △ ( x , x ′ ) is not the canonical 1 /r -potential as it lacks the required periodicity andonly applies under free space boundary conditions,lim | x |→∞ | x | V ( x ) = − π , see [36]. Instead, we have: G ppp △ ( x , x ′ ) = 1 L L L X k n k > − k e i k · ( x − x ′ ) , (10)with n ∈ Z , k ∈ R and k i = πL i n i .The Newtonian 1 /r -potential may be recovered as free-space limit in two stages. By first taking L , → ∞ [37]recasts Eq. (10) into: G ppp △ ( x , x ′ ) L , →∞ −−−−−→ G ffp △ ( x , x ′ ) = 12 πL log | x ⊥ − x ′⊥ | − πL ∞ X m =1 cos ( k m ( x − x ′ )) K (cid:18) k m | x ⊥ − x ′⊥ | (cid:19) , (11)with K ( x ) denoting the 0 th -modified Bessel functionof the second kind and x ⊥ = ( x , x ) ⊺ . We returnto this mixed boundary condition Green’s function inSec. II B 2 b. Finally, take L → ∞ so that the Riemannsum in Eq. (11) approaches an analytically solvable inte-gral, [38]: G ffp △ ( x , x ′ ) L →∞ −−−−→ G fff △ ( x , x ′ ) = − π | x − x ′ | , yielding the expected result.The naive way of carrying out the dimension reductionof Eq. (9) is to simply drop all partial derivatives in x , x -direction. This appears to be the common approach inlow-dimensional studies on FDM [22, 24, 33, 39–41] andturns out to be true assuming we demand a uniform mat-ter distribution along the neglected dimensions. The ap-proach is equally applicable for d = 1 , d + 1)-SP equation.Maintaining Poisson’s equation as field equation hasimplications on how gravity acts in lower dimensionssince the periodic Green’s function in Eq. (10) dependson the dimensionality of the Laplace operator. Westress even if the aforementioned violation of the periodicboundary conditions would not exist for the Newtonianpotential, it would still be impossible to simply enforce a1 /r -interaction kernel in one dimension. Its singularity atthe origin remains too strong and consequently yields anill-defined convolution kernel. Hence, we ask whether areduction exists that approximately preserves the threedimensional interaction with only one spatial degree offreedom. This is realized by strongly confining matterorthogonal to the dimension in which the evolution isobserved.
1. General Reduction Procedure
We adapt the discussion under free-space conditionsoutlined in [42] to the periodic situation at hand. Tothis end, the NLSE in Eq. (9) is augmented by an artifi-cial, external potential V ext ( x ⊥ ; ǫ ) = ǫ V ext (cid:0) x ⊥ ǫ (cid:1) in theorthogonal x ⊥ -plane that is controlled by a confinementparameter ǫ . The extended Hamiltonian H = H x + H ǫ ⊥ is then decomposed into: H x = − ∂ x + a ( t ) (cid:16) G p △ ∗ | Ψ | (cid:17) ,H ǫ ⊥ = 1 ǫ (cid:20) − △ e ⊥ + V ext ( x e ⊥ ) (cid:21) = 1 ǫ H e ⊥ , with x e ⊥ = x ⊥ /ǫ . Notice how we keep the boundaryconditions of the interaction kernel G p △ deliberately un-specified in the x ⊥ -plane. These will be set in accordance with the external potential, i.e. if V ext ( x ⊥ ) is periodic onΩ × Ω , the long range interaction is G ppp △ , as in Eq. (10).If, however, the external potential assumes free-space con-ditions, the reduction departs from the mixed conditionGreen’s function G ffp △ in Eq. (11).Next, find the eigensystem { χ ǫk , λ ǫk } to the linear eigen-value problem H ǫ ⊥ χ ǫk ( x ⊥ ) = λ ǫk χ ǫk ( x ⊥ ) and enforce a fac-torization of the full fledged wave function Ψ( x , t ) accord-ing to Ψ( x , t ) = ψ ( x , t ) χ ǫk ( x ⊥ ) e − iλ ǫk t . Insert this ansatzinto the NLSE with the extended Hamiltonian, multiplyby ( χ ǫk ) ∗ and integrate over the x ⊥ -plane. The result isa one dimensional NLSE for ψ ( x , t ): i ~ ∂ t ψ = (cid:20) − ∂ x + a ( t ) (cid:0) U p ∗ | ψ | (cid:1)(cid:21) ψ x ∈ Ω , (12)alongside a one dimensional, long range interaction, U p : U p ( x , x ′ ) = 1 k χ ǫk k Z Ω × Ω d x ′⊥ Z Ω × Ω d x ⊥ | χ ǫk ( x ⊥ ) | | χ ǫk ( x ′⊥ ) | × G p △ ( x ⊥ , x , x ′⊥ , x ′ )(13)Since the spectrum of Eq. (10) is spherical symmetricaround k = 0 it is clear that, irrespective of the confiningpotential in the transversal plane, the one dimensionalinteraction kernel remains even and depends only on therelative distance: U p ( x , x ′ ) = U p ( | x − x ′ | ).
2. Uniform vs. Confined Transversal Density
Let us restrict the discussion to two distinct externalpotentials that either induce a complete delocalization ofmatter in the orthogonal plane or constrain all matteraround the x -direction. In what follows, both of thesemodels will be coined (1 + 1)-FDM collectively. If a dis-tinction has to be made, a more precise denotation willbe used. a. Uniform Matter Set V ext = 0. The correspond-ing eigenstates are plane waves and matter is thereforeassumed to be organized in homogeneous matter sheetsparallel to the x ⊥ -plane. Moreover, the external poten-tial is trivially L -periodic so that G p △ = G ppp △ . Eq. (13)then evaluates to, [38]: G p △ ( x , x ′ ) = 12 | x − x ′ |− (cid:20) ( x − x ′ ) L + L (cid:21) , (14)i.e. the Green’s function of the d=1 Poisson equationunder periodic boundary conditions.Eq. (12) together with the interaction of Eq. (14) isdenoted (1 + 1)-SP. b. Confined Matter We enforce an integrable har-monic confining potential V ext ( x ⊥ ) = x ⊥ and consideronly the ground state χ ( x ⊥ ) ∝ exp( − x ⊥ /
2) for thedynamics in the orthogonal plane. Our choice is moti-vated by cigar-shaped confinements often used for trap-ping Bose-Einstein condensates, see e.g. [3, 43–47]. Nei-ther the external potential nor its ground state are peri-odic. Hence, we set G p △ = G ffp △ in eq. (13) and find: U pconf ( x , x ′ ) = − πL X m U (cid:18) , , k m ǫ (cid:19) e ik m ( x − x ′ ) , (15)with m ∈ Z \ U ( a, b, x ) denoting the confluenthypergeometric function of the second kind, [48].Consider Fig. 1 for a graphical comparison of the longrange interaction induced by a test particle at x ′ = 0with and without confinement. All potentials are shownfor a periodic box of L = 100. Apart from the requiredfiniteness of both kernels at the origin, the behavior of U pconf and G p △ is quite disparate in the far-field region.A key difference lies in the effective interaction range R which may be defined as: ∂ r U p ( R ) ≡ δ max | ∂ r U p | with 0 < δ ≪ . (16)While the gravitational potential under confinement israther localized with R ∼ ǫ and quickly approaches thedesired Newtonian potential (black, dashed line) at largedistances | x − x ′ | , the effective interaction range for G p △ evaluates to R = (1 − δ ) L and is therefore comparableto the box size. As discussed in Sec. IV B 4, this hasimplications on which asymptotic states are accessible inthe (1 + 1)-SP case.
3. Symmetries and Conserved Quantities
Naturally, one is interested in conserved quantities ofEq. (12) for even kernels but irrespective of its exact form.By analyzing the action: S = Z d t Z Ω d x ( iψ ∗ ∂ t ψ − H [ ψ, ψ ∗ , ∂ x ψ, ∂ x ψ ∗ , t ]) , (17)with the Hamiltonian density: H = 12 (cid:0) | ∂ x ψ | + a ( t )( U p ∗ | ψ | ) | ψ | (cid:1) generating Eq. (12) upon variation, it is straightforwardto show that such a system obeys mass and momen-tum conservation. Moreover, Eq. (17) is invariant undergalilean boosts of the form: ψ ( x, t ) → e i ( vx − v t ) ψ ( x − vt, t ) . (18)If we consider a static space-time, i.e. a ( t ) = const . ,Eq. (17) is also time translation invariant and the total − − − π | x − x ′ | .
05 0 . x − − U p confined ǫ = 0 . ǫ = 0 . ǫ = 0 . FIG. 1. (Color online) Comparison of the one dimensional,long range interactions of Sec. II B 2 a and Sec. II B 2 b in-duced by a source particle located at x = 0 in a periodicbox of L = 100. The negative x-space illustrates the sit-uation for a uniform matter distribution in the transversalplane with G p △ , Eq. (14), as gravitational potential (blueline). The positive x-space depicts the gravitational potential U pconf , Eq. (15), obtained under harmonic confinement for var-ious confinement strengths ǫ (red lines). While both G p △ and U pconf are finite at the origin only U pconf approaches the New-tonian gravitational potential (black, dashed line) in the farfield. Note the different limits of the negative and positive x -axis: The negative half space covers the interval [ − L / , x = 0 .
15 to make theasymptotic behavior of U pconf better observable. energy, E [ ψ ] = Z Ω d x H [ ψ, ψ ∗ , ∂ x ψ, ∂ x ψ ∗ , t ] = h T i + a h V i , (19)is conserved as well. This is of course not true once space-time is allowed to expand.A symmetry unique to (1 + 1)-SP is the followingscaling transformation: If ψ ( x , t ) solves Eq. (12) with U p = G p △ then so does: e ψ ( x , t ) = λ ψ ( λx , λ t ) , λ ∈ R + . (20)Equivalent scaling symmetries exist for d = 2 ,
4. Properties of the (1 + 1) -FDM Ground State
We already alluded in our introductory remarks to theimportance of stationary states of (3 + 1)-SP, especiallythe role of its stable ground state that acts as dynamicalattractor and realizes a flat-density core inside relaxeddark matter structures. Let us extend this discussion tothe previously derived one dimensional FDM models byanalyzing properties of the ground state to Eq. (12) thatinfluence the asymptotic (1 + 1)-FDM dynamics. − − x | ψ | τ = 0 τ = 0 . τ = 0 . τ = 0 . τ = 1 (A) . . . . . τ − − − − E / E G S E h T ih V i h T i − h x∂ x V i (B) FIG. 2. (Color online) Exemplary gradient descent for U p = U pconf , i.e. under confinement, with ground state mass M = 50, a = 1 and ǫ = 10 − . Panel (A): Densities at variousstages of the minimization procedure. Panel (B): Evolutionof kinetic (red), potential (yellow) and total energy (black)normalized to the ground state energy E GS . Also shown ingreen is the free-space, quantum virial theorem, discussed inSec. II C 1. For practical purposes, ground states of mass N [ ψ ] = R d x | ψ | = M may be prepared by choosing an inter-action kernel and minimizing the grand canonical energy E grand [ ψ ] = E [ ψ ] + µ ( τ ) N [ ψ ] by means of a gradient de-scent with µ ( τ ) as chemical potential at descent parame-ter τ . The reader is referred to [49, 50] for numerical de-tails and [51] for a rigorous analysis on the existence anduniqueness of a minimizer to Eq. (19) for (1+1)-SP underfree-space conditions. We note in passing that the chosennumerical implementation of the gradient descent makesthe energy minimization approach equivalent to the well-known imaginary time propagation method [47, 52].An exemplary gradient descent under transversal, har-monic confinement is depicted in Fig. 2. As time pro-gresses the initial gaussian distribution focuses more andmore in position space until the descent converges to aspatially localized structure at τ = 1.Irrespective whether U p = U pconf or U p = G p △ , the(1 + 1)-FDM ground state shows the following properties: a. Solitary Wave It is easy to see that theequation governing the gradient descent, namely ∂ τ ψ = − δδψ E grand , reduces to a stationary form of Eq. (12), µ ( τ ∗ ) ψ GS = (cid:20) − ∂ x + a ( U p ∗ | ψ GS | ) (cid:21) ψ GS , (21)once the energetic minimum at τ = τ ∗ is reached.Therefore, the ground state may be written in itscanonical, linear quantum mechanics form, ψ GS ( x , t ) = ψ GS ( x ) e − iµ ( τ ∗ ) t , and we conclude ψ GS is a solitary wave ,i.e. a localized solution to a nonlinear equation with timeindependent envelope | ψ | . Obviously, its persistent formis also preserved for uniformly travelling configurationsobtained via Eq. (18).Two remarks are in order: Firstly, the above discus-sion derives ψ GS by minimizing the total energy func-tional. An alternative approach, e.g. [25], is to interpretthe nonlinear eigenvalue problem of Eq. (21) as bound-ary value problem and solve it via a shooting method.Secondly, recall the scaling symmetry of Eq. (20) uniqueto ( d + 1)-SP. Applying it to a single ground state givesaccess to an entire family of ground states parameterizedby their total mass M . Thus, for (1 + 1)-SP, it is suffi-cient to conduct the energy minimization only once foran arbitrary reference mass M ref . All other energy mini-mizers then follow from rescaling with Eq. (20). Lackingan equivalent symmetry, this is of course not true for theconfinement model. b. Inelastic Collisions So far, we demonstrated thesolitary character of the (1+1)-FDM ground state. Natu-rally, we are interested whether ψ GS can also be identifiedas a solitonic solution.Strictly speaking, the concept of a soliton calls for a rig-orous mathematical definition. To keep technical detailsto a minimum, we instead follow [53] and characterize asoliton as a solitary wave that is invariant under inter-actions with other solitons. Put differently, despite thenonlinear evolution, solitons obey a superposition prin-ciple and neither mass nor energy should be exchangedduring the interaction.Figure 3 investigates the behavior of an asymmetricconfiguration of two confined ground states with initialmasses M = 50 and M = 100 and ǫ = 10 − , boostedonto a collision course with: ψ ( x ) = ψ GS , ( x + x ) e ivx + ψ GS , ( x − x ) e − ivx . (22)Figure 3(A) depicts the matter density in a pre-collision time window, an instance during the groundstate interaction ( t ≈ .
4) as well as a post-collision timewindow around t = 32. It is evident that the superposi-tion principle is not satisfied. After the interaction tookplace both ground states propagate in a quasi-solitaryfashion in which linear dispersive and nonlinear focusingeffects are not exactly balanced anymore. Instead onefinds a periodic expansion and re-contraction of the mat-ter distribution once the dispersion or non-local, nonlin-earity dominates. Similar oscillatory behavior was foundfor the (3 + 1)-SP ground state, [25], once the density is − − − −
10 0 10 20 30 40 x | ψ | Light GS( M = 50) Heavy GS( M = 100) (A) t = 0 t = 32 . t ≈ . − . . . ∆ M ∆ M = 0 . (B) Heavy GSLight GS0 5 10 15 20 25 30 t − ∆ E ∆ E = 10 . (C) Heavy GSLight GS
FIG. 3. (Color online) Inelastic interaction of an asymmetric high-mass, low-mass ground state configuration under strong,transversal confinement ( ǫ = 10 − ) and a = 1. The evolution starts from Eq. (22). Panel (A): Density evolution. Initially, bothdensities travel as solitary waves (red), pass through each other (black, inset) and continue to propagate in a quasi-solitarymovement after the interaction (orange). By this we mean a state for which neither linear dispersive nor nonlinear focusingeffects induce a permanent deformation of the density. Instead, one observes an oscillation around a solitary wave. Similaroscillatory behavior was found for (3 + 1)-SP, [25], once the ground state density is perturbed. Panel (B): Time evolution ofthe mass deviation of both ground states relative to their initial masses cf. Eq. (23). Post-interaction, it is the high-massground state gaining additional matter at the expense of the lighter ground state. Panel (C): Time evolution of the total energydeviation of both ground states relative to their respective initial total energies. As for the mass, the interaction induces anenergy transfer from the low to high mass ground state. Notice, how the quasi-solitary behavior of the post-interaction densitiesis also seen in an oscillation of the energy deviation. Since the association of mass and energy contained in the positive andnegative box half to a particular solitary wave is ambiguous during the interaction, we deem data in the gray shaded intervalof panel (B) and (C) as not reliable. perturbed. Nevertheless, on average the post-interactionconfiguration is still comprised of two solitary, or station-ary, states.To investigate whether these post-interaction, solitarywaves are different from their initial composition, we an-alyze the deviation in mass (Fig. 3(B)) and total energy(Fig. 3(C)) from their initial values. For instance, themass difference for the M -ground is inferred as:∆ M ( t ) = L Θ( t − t coll ) Z − L Θ( t coll − t ) d x | ψ ( x , t ) | − M , (23)with t coll as collision time naively inferred from the uni- form velocity v at t = 0.One finds, a symmetric mass and energy gap after theinteraction: Both mass and energy were transferred fromthe low to high mass solitary wave. Clearly, such a mat-ter and energy transfer should not exist if the confinedground state were a true soliton. We note although thereported energy and mass differences are small they arerobust under variation of ∆ t , see Sec. III B for more in-formation. Qualitatively similar results were found for(1 + 1)-SP. Therefore, (1 + 1)-FDM ground states — atleast the ones considered here — are not solitons in thestrict sense of the word, but interact inelastically by ex-changing mass and energy during encounters, typicallyreshuffling them from the low-mass to the high-mass soli-tary wave.Of particular interest is the case of multiple successiveinteractions which, thanks to the periodicity of the box,is easily observed by increasing the integration time. Thereader is referred to Sec. IV B 2 for more details. c. Mass-Size Relation As we will see, understandingthe discrepancies between the attained asymptotic statesof (1 + 1)-SP and the confinement model, Sec. IV B 4,hinges on the ratio R ( L ) /σ ( M ), i.e. the interactionrange R given a periodic box of size L compared to thespatial extent σ of a mass M ground state.Deriving σ ( M ) is particularly simple in case of the un-confined FDM model as Madelung’s ansatz relates (1+1)-SP to a one dimensional version of the hydrodynamicdescription of Eq. (6a)-(6d). In the ground state’s restframe these reduce to the condition of hydrostatic equi-librium. Dimensional analysis then yields σ ∝ M − .The situation is more involved under harmonic confine-ment due to the missing PDE for the gravitational poten-tial. Thus, the spatial extent is deduced numerically bydefining: 0 . M ≡ Z σ − σ d x | ψ GS | , (24)and extracting σ for various, spatially centered groundstates of mass M . Figure 4 depicts the result for a staticspace-time with a = 1 for both (1 + 1)-FDM models.While (1 + 1)-SP shows satisfactory agreement with thedimensional analysis, σ ( M ) = 2 . · M − . , a stronglyconfined matter density at ǫ = 10 − results in a narrowerground state distribution at equal mass M . In this case σ ( M ) = 4 . · M − . . C. Relaxation Mechanisms and EquilibriumProperties
It is a priori not clear what dynamical mechanismsdrive (1 + 1)-FDM into its asymptotic equilibrium con-figuration let alone whether both reduction models obeythe same relaxation processes — recall the discrepanciesin the interactions of (1 + 1)-SP and the confinement sce-nario.Given the approximative CDM interpretation of FDMin Sec. II A 2 b, classical, non-collisional relaxation mech-anism may be a viable option, in particular a combina-tion of phase mixing and violent relaxation, see [54, 55].These processes induce a filamentation of the phase spacedynamics alongside a redistribution of energy inside theself-gravitating structure due to its fluctuating gravita-tional potential.On the other hand, (3 + 1)-FDM-typical mechanismslike gravitational cooling, [29], may be recovered even inone dimension, allowing collapsing matter structures torelax into an equilibrated state by radiating away excessenergy in form of small scale matter waves. M − σ uniformconfined ( ǫ = 0 . FIG. 4. (Color online) Spatial extent of the one dimensionalFDM ground states as a function of mass M at a = 1. Eachdata point corresponds to a solitary wave prepared by thegradient descent shown in Fig. 2. The spatial extent of thematter distribution is extracted according to Eq. (24). For U p = G p △ (blue) we find σ ( M ) = 2 . · M − . and conse-quently good agreement with dimensional analysis. Understrong harmonic confinement, i.e. U p = U pconf and ǫ = 10 − (red) one deduces σ ( M ) = 4 . · M − . over two orders ofmagnitude in M . Which relaxation channels are realized is discussed inSec. IV B 2. Here, we ask what properties the equili-brated system configuration should have and how theymay be measured such that the progress on the overallsystem evolution can be quantified.
1. Virial Equilibrium
Application of Ehrenfest’s Theorem for the virial opera-tor ˆ G = ˆ p ˆ x gives rise to a quantum analogue of the scalarvirial theorem. For bounded dynamics, i.e. h ˆ G i ( t ) < ∞ ,and periodic boundary conditions it reads:0 = 2( h T i ) ∞ − ( a h x ∂ x V i ) ∞ + L (cid:0) ( ψ ∗ (0 , t ) ∂ x ψ (0 , t ) − | ∂ x ψ (0 , t ) | (cid:1) ∞ (25)with ( A ) ∞ = lim t →∞ R t d t ′ A ( t ′ ). Relaxation into virialequilibrium, i.e. the regime where Eq. (25) is (approx-imately) satisfied, is then to be understood as a con-sequence of the evolution under Schr¨odinger’s equation.That said, any finite quantum system would virialize inthe limit t → ∞ .A couple of remarks are in order. Firstly, note Eq. (25)only holds in the limit t → ∞ . A notable exception arestationary states, like the (1 + 1)-FDM ground states ofSec. II B 4, which obey Eq. (25) without time averaging cf.[56]. If in addition fluctuations are present, we may assessvirialization of the total system by assuming Eq. (25)were approximately achieved after a finite thermalizationtime.0Secondly, we draw special attention to the boundaryterm, B ( t ) = ψ ∗ (0 , t ) ∂ x ψ (0 , t ) − | ∂ x ψ (0 , t ) | , (26)in Eq. (25). It emerges from the necessity to extend thedomain of the Hamiltonian onto states like xψ which arenot periodic but appear once Ehrenfest’s theorem is ap-plied to ˆ G , see [57]. Obviously, the boundary term is neg-ligible, if ψ decays rapidly towards the box boundaries,like in Fig. 2, or when artificial absorbing potentials areused to limit the physically relevant domain size, see e.g.[17, 27].
2. Maximum Entropy
From a statistical physics viewpoint, one generallyexpects the system to maximize its entropy. In fact,the idea of entropy maximization is close to the origi-nal approach of [54], showing that mixing processes un-der Vlasov-Poisson imply a quasi-stationary phase-spacedistribution that maximizes the system’s entropy on amacroscopic, i.e. coarse-grained level.To adopt this idea for FDM, we use Husimi’s distri-bution, i.e. ¯ f W in Eq. (5) and define a FDM entropyfunctional resembling Boltzmann’s entropy, [58]: S [ ¯ f W ] = − π Z d x d k ¯ f W log ¯ f W . (27)Thus, an equilibrated, thermalized system state isreached, once ∆ S = S ( t ) − S (0) saturates. We note thisapproach was also proposed by [24]. III. NUMERICAL METHOD
We briefly explain a simple, yet accurate spatial dis-cretization of Eq. (12) and sketch an approximation toits time evolution operator. The main properties of thepresented method are summarized. For more informa-tion, especially on the method’s behavior under expand-ing space-time conditions, the reader is referred to Ap-pendix A. Additional information on general NLSE nu-merics can be found in [59]. A recent survey of existingnumerical techniques on our subject is given by [60].
A. Spatial Discretization
Since Eq. (12) is defined on a periodic domain and in-volves only second derivatives in space, expansion of ψ ina truncated momentum eigenstate basis is a natural wayto discretize Eq. (12) in momentum space and diagonal-ize the kinetic part of the Hamiltonian. Discreteness inreal space is then achieved by evaluating the momentumstate expansion of ψ on sites { x j } j =0 ...N − with uniform spacing ∆ x = L/N . This translates Eq. (12) into thefinite dimensional, ordinary differential equation: i∂ t ψ ( t ) = " F † k k ⊺ F | {z } ˆ H K + a ( t ) V ( | ψ ( t ) | )] | {z } ˆ H V ( t, | ψ | ) ψ ( t ) , (28)with ψ j = ψ ( x j , t ), ( k ) n = πL n and F denoting thechange of basis matrix from the real space to the momen-tum basis. In practice the action of F and its inverse F † on ψ are implemented as discrete fast Fourier transform.The nonlinear, non-local potential is absorbed into thediagonal matrix V ( | ψ | ) and follows directly from theconvolution theorem,Diag (cid:2) V ( | ψ ( t ) | ) (cid:3) = F † d U π F| ψ ( t ) | , with the, in momentum space diagonal, kernel coefficientmatrix d U π :Diag hd U π i = n = 0 − k ) n n = 0 (1 + 1) SP − π U (1 , , ( k ) n ǫ ) n = 0 confined . B. Time Evolution Operator
Starting from Eq. (28) it remains to find an approx-imation to the time evolution operator ψ ( t + t ) =ˆ U K + V ( t , t + ∆ t ) ψ ( t ). For this, the idea of operatorsplitting is employed — a common choice for integratingNLSEs. Thus, we first find (approximate) evolution op-erators for the kinetic, ˆ H K , and potential Hamiltonian,ˆ H V , individually and combine them into an approxima-tion for ˆ U K + V afterwards.The solution to the kinetic problem is trivial and reads:ˆ U K (∆ t ) ≡ ˆ U K ( t , t + ∆ t ) = F † exp (cid:20) − i k k ⊺ t (cid:21) F . For the potential sub-problem we recall ˆ H V isdiagonal in real space, implying [ ˆ H V ( t ) , ˆ H V ( t ′ )] = 0.Thus, its time evolution operator can be written asˆ U V (∆ t ) = exp (cid:2) − i R t +∆ tt d t ˆ H V ( t, | ψ ( t ) | ) (cid:3) and it re-mains to approximate the time integral over ˆ H V .Fortunately, it is easily verified that evolutionunder the nonlinear Hamiltonian ˆ H V satisfies dd t | ψ | = 0. It is therefore sufficient to substituteˆ H V ( t, | ψ ( t ) | ) → ˆ H V ( t, | ψ ( t ) | ) which (i) reduces thetask of approximating the time integral over ˆ H V toapproximating the integral over the scale factor a ( t ) and(ii) allows for an explicit treatment of the nonlinearity,i.e. without the need to solve a nonlinear system ofequations.1Application of the midpoint rule yields the followingunitary approximation to ˆ U V (∆ t ):ˆ U V (∆ t ) = ˆ U V (∆ t, t ) + O (∆ t )= exp (cid:20) − ia (cid:18) t + ∆ t (cid:19) V ( | ψ ( t ) | )∆ t (cid:21) + O (∆ t ) . Finally, after composing both evolution operators in asecond order Strang scheme, we arrive at the approxima-tion to ˆ U K + V :ˆ U K + V = ˆ U K (cid:18) ∆ t (cid:19) ◦ ˆ U V (∆ t, t ) ◦ ˆ U K (cid:18) ∆ t (cid:19) + O (∆ t ) . C. Summary of Properties
The presented method is a simple extension to the kick-drift-kick scheme of [17]. In fact, for static space-timesboth methods are equivalent. Aside from its implementa-tional simplicity and resemblance of the symplectic Leap-frog method, it is unitary by design, explicit, second-order accurate in time, spectrally accurate in space (as-suming ψ is smooth) and provides a convenient unifiedapproach for both dimension-reduction models. The com-putational complexity per integration step is O ( N log N )due to the fast Fourier transformations, requiring O ( N )memory.Furthermore, for time-independent, nonlinear couplingconstants, the approximate evolution operator is time-symmetric, shows unconditionally stable numerical be-havior and conserves energy, Eq. (19), approximatelywith a bounded error, [61].For non-static background cosmologies, time-symmetry and energy conservation are broken bythe continuous problem. Concerning stability, our testsindicate an exponentially growing error at high redshifts.We expect this result to be to intrinsic to constant timestep integration methods under space-time expansion.Further informations on the convergence properties ofour numerical method are given in Appendix A. IV. RESULTS OF NUMERICAL SIMULATIONSA. Structure Growth under (1 + 1) -SP
We investigate the mean cosmic structure growth inan ensemble of FDM-only universes obeying (1 + 1)-SP.The purpose of this study is to explain characteristicproperties of the nonlinear FDM matter power spectrum P ( k ) = h| ˆ δ k | i .
1. Simulation Setup
To this end, we follow the evolution of N = 100 re-alisations of a gaussian random field δ ( x, a ) in a flat, radiation free FLRW background cosmology with Ω m =Ω DM + Ω baryon = 0 . H = 68 km s − Mpc − and powerspectrum: P ( k, a ) = D ( a ) T ( k ) T ( k ) P CDM ( k ) , (29)with k = | k | . Here, D ( a ) denotes the linear growthfactor, see [62], normalized to unity at z = 0, P CDM ( k )the linear CDM power spectrum produced by CAMB , see[63], at redshift z = 0, T (k) the CDM to FDM trans-fer function of [12] and T ( k ) = k π a transfer functionreducing the spectrum’s dimensionality to one spatial de-gree of freedom. The latter follows from demanding adimension independent real space variance.The initial phase function, S ( x, a ), is obtained fromEq. (6b) by solving: ∂ x S ( x, a ) = − ma H ( a ) δ ( x, a ) . Once { δ ( x, a ) , S ( x, a ) } are known, the initial wave func-tion follows from Madelung’s ansatz, see II A 3.Each realisation starts from z = 100 and is integrateduntil z = 0. For the FDM mass two fiducial values,namely m = 10 − eV and m = 10 − eV are analyzed.To guarantee sufficient resolution of P ( k ) at z = 0, thenumber of uniform spatial grid points is set to N = 2 ,implying for L = 100 Mpc a step size of ∆ x = 23 .
2. Overall Evolution of the Matter Power Spectrum
Figure 5 illustrates the evolution of the matterpower spectrum for m = 10 − eV, Fig. 5(A), and m = 10 − eV, Fig. 5(B), at various redshifts z . In bothpanels, the solid black line represents the linearly rescaledand dimensionally reduced reference spectrum of Eq. (29)from which the initial conditions of each realisation aredrawn. It is characterized by a flat, large scale regimequickly transitioning into a steep power law suppressionaround k J ( a ). Its exact functional behavior is encapsu-lated in the FDM transfer function T FDM ( k ).For both mass parameters the evolution of P ( k ) maybe summarized as follows: Early on, all modes behavelinearly, i.e. evolve according to Eq. (7). Recall that com-plex modes ˆ δ k with k > k J ( a ) are stabilized by quantumpressure and are therefore confined to a damped, oscil-latory motion with no increase magnitude | ˆ δ k | . Conse-quently, as long as all modes evolve linearly, one expectsthe power spectrum to stay close to its initial shape for k ' k J ( a ). The situation at z = 40 recovers this behav-ior for both mass parameters but is best seen in Fig. 5(B)in which P ( k ) slowly detaches itself from the power lawsuppression regime for k < k J . Recall linear FDM modesevolve independent but differently. It is therefore no sur-prise that the initial shape of P ( k ) is lost.As time progresses, k J proceeds to propagate outwarduntil all modes of interest are destabilized and collapseunder their own gravity. It is then that the linearizeddescription of Eq. (7) breaks down and nonlinear mode2 − − − − − P [ M p c ] k J ∝ k − z = 0 z = 1 z = 2 z = 5 z = 10 z = 40 (A) k [Mpc − ] − − − − − P [ M p c ] k J ∝ k − (B) FIG. 5. (Color online) Evolution of the matter power spec-trum P ( k ) = h| ˆ δ k | i inferred from the simulation ensemblespecified in Sec. IV A 1 with m = 10 − eV in (A) and m = 10 − eV in (B). The black solid lines show the initialpower spectrum, Eq. (29). After small scales pass the timedependent Jeans scale k J ( a ), their associated perturbationmodes ˆ δ k quickly grow in magnitude and couple to other non-linearly evolving modes. At late times, this leads to a distinctshape of P ( k ) with two characteristic regimes: At high k thecomoving uncertainty principle induces a sharp power sup-pression past k s ( a ), see Eq. (30) and dashed, vertical lines.By contrast, modes with intermediate values of k , induce ascale free power spectrum with P ( k ) ∝ /k cf. Sec. IV A 4. coupling is expected to set in — the independence of eachˆ δ k is lost. Driven by the focusing, nonlinear interaction,the result is a redistribution of matter power across allnonlinearly evolving perturbation modes. This manifestsitself in two observable effects in Fig. 5(A)/(B) for allredshifts z ≤
10: Firstly, an intermediate coupling regimeemerges that is well described by P ( k ) ∝ k − . Secondly,the power suppression regime steepens even further andcontinues to travel outward, leaving a distinct cutoff in P ( k ) at high k . Changing the mass parameter influencesthe transition scale between both regimes.
3. The Suppression Scale
We recall from Sec. II A 3 that the Jeans length onlyapplies in the linear regime and can therefore not explainthe observed transition from matter coupling to suppres-sion. On the other hand, the uncertainty principle re-mains applicable even under nonlinear evolution. Thus,we adapt the heuristic argument of Sec. II A 3 and againidentify the hydrodynamic velocity dispersion as a mea-sure for the velocity uncertainty σ v , but this time inferthe dispersion directly from the simulation. The suppres-sion scale k s then follows from: k s ( a ) = 2 πσ x ≈ πa m ~ p h v i − h v i (30)cf. Eq. (6a). Figure 5 illustrates its ensemble average h k s i ( a ) as vertical dashed lines for z <
5. The correspon-dence between h k s i ( a ) and the true suppression scale isconvincing. We conclude that the small scale evolutionremains well explained by the uncertainty principle evenin the nonlinear evolution regime.
4. The Coupling Regime
For scales k < k s ( a ) one expects FDM to quickly re-cover the evolution of cold dark matter (CDM). The ob-served power-law behavior of P ( k ) at scales k box ≪ k We now shift our attention to the subject of asymptoticdynamics and its equilibrated final (1 + 1)-FDM states.To restrict the complexity of the analysis two simplifica-tions are introduced.Firstly, artificial, i.e. non-cosmologically motivated,initial conditions are employed as these allow us to freelyset the degree of spatial localization. This is impor-tant since spatially delocalized initial conditions, such3as the gaussian random fields of Sec. IV A, usually comewith multiple overdense regions that collapse into multi-ple, high mass clusters which then undergo a subsequentmerger and collision phase. It is these violent, late timeevents that complicate the dynamics unnecessarily sincethey drive already relaxed clusters again out of equilib-rium, therefore increasing the required integration timeto re-relax into the asymptotic state. Moreover, startingfrom a localized configuration reduces the time to firstcollapse.Secondly, to assure swift relaxation times the scale fac-tor, which acts as a coupling constant to the nonlinearinteraction, is fixed to a = 1. Extensions of our staticspace time results to the expanding FLRW scenario aredeferred to section IV B 5. 1. Simulation Setup More precisely, all simulations reported in this sectiondepart from a gaussian initial density with zero initialvelocity: | ψ ( x , | ∝ e − x σ , Arg[ ψ ( x , . The standard deviation is chosen as σ = 6 k − J (1) cf.Eq. (8) assuring instability of the entire spectrum of ψ right from the beginning of the evolution. In or-der to comply with the periodic boundary conditionsup to floating point precision, the box size is chosen as L = 30 σ ≈ 127 and the numerical study is conductedfor both (1 + 1)-SP and the harmonically confined reduc-tion model. In both cases the number of grid points waschosen such that the entire wave function spectrum | ψ k | stays resolved throughout the integration. The integra-tion is stopped some time after Eq. (25) and/or Eq. (27)indicate the completion of the virialization and/or ther-malization process.All data is reported in dimensionless quantities. Toget some sense of scale, choosing the canonical FDMmass m = 10 − eV and adopting identical cosmologi-cal parameters as in Sec. IV A 1 implies a box size of L ≈ . 2. Relaxation a. (1 + 1) Schr¨odinger-Poisson Figure 6 illustratesthe relaxation process under (1 + 1)-SP in multiple ob-servables.For t / x = L / 2, (iii) decelerating until a turn-around radius ishit and finally (iv) recollapsing towards the origin. Thesecycles do not occur in a strictly sequential manner but areincreasingly superposed and thus induce a characteristicspiralization of the phase space distribution in Fig. 6(B).Notice that space regions exist in which multiple inwardand outward propagating matter streams coexists simul-taneously.The structure of the early phase space distributions,Fig. 6 (A)/(B) is qualitatively in good accordance withthe evolution of one dimensional collisionless N-body sys-tems, e.g. [55], and furthermore show the natural sig-nature of phase mixing and (less pronounced) violentrelaxation. While phase mixing manifests itself in theever tighter spiralization of Husimi’s distribution, vio-lent relaxation induces a small yet observable increasein the occupied phase space volume. The expansionin k -direction is best seen by comparing Fig. 6(A)/(B)whereas the white dotted lines in Fig. 6(D) show thespatial expansion. As for collisionless N -body systems,reason for this expansion is the time-dependency of thegravitational potential which dies out quickly after only1 − ψ when matterstreams intersect and are therefore maximally localizedin space. Once matter flows outward the system is lessbound thus increasing the expectation value of the poten-tial energy.Only 3 − t / 5. Thus, the thermaliza-tion and virialization time scale are essentially identicaland both metrics capture the convergence into the asymp-totic state equally well in this example. Furthermore, theboundary term B , Eq. (26), is negligible in this setting asthe long interaction range of G p △ does not allow ejectedmatter clumps to propagate till the domain boundaries.Hence ψ and its derivatives persist to be small at x = 0. b. Strong Harmonic Confinement Figure 7 depictsthe situation under strong, harmonic confinement with ǫ = 0 . 01. Again, a relaxation and quasi-stationary phasemay be identified. Their respective duration, however, isopposite to (1 + 1)-SP.During relaxation, i.e. when t / | ψ | in Fig. 7(D), around t ≈ 2. There, multiple, stable density excitations ofvarious masses depart from the position of first collapse,propagate outward, overcome the central gravitationalpotential at x = L / 55 65 75 x − − k (A) 55 65 75 x − − k (B) 55 65 75 x − − k f H f H, max (C) t x | ψ | (D) − t S (G) − − − E / | E G S | | (2 h T i − a h x∂ x V i ) t | | (2 h T i − a h x∂ x V i + B ) t | (F) − − E / | E G S | E h T i h V i (E) − − − t = 0 . t = 3 t = 1000 FIG. 6. (Color online) Relaxation process of an unstable gaussian density profile under (1 + 1)-SP visualized in multipleobservables. Panel (A)-(C): Husimi phase space distribution, Eq. (5), at characteristic stages of the evolution. For all panelsa constant smoothing scale of σ x = 1 / √ | ψ | . No qualitatively new features appear in the densitypast t ≈ x -space due to violent relaxation (see main text). Panel (E) Energy evolution. Theconservation of Eq. (19) is apparent. Panel (F): Assessment of the relaxation process in terms of the quantum virial theorem,Eq. (25). The absolute deviation from Eq. (25) decays quickly until t ≈ t ≈ 5. Comparing (F) and (G) reveals thatvirialization and thermalization occur on the same time scale. All energies were normalized to a M = L (1 + 1)-SP groundstate. 55 65 75 x − − k (A) 55 65 75 x − − k (B) 55 65 75 x − − k f H f H, max (C) 200 300 400 500 600 700 800 900 t | ψ | t S (G) − − E / | E G S | | (2 h T i − a h x∂ x V i ) t | | (2 h T i − a h x∂ x V i + B ) t | | ( B ) t | (F) − − E / | E G S | E h T i h V i (E) − − − x (D) t = 8 t = 352 t = 1000 FIG. 7. (Color online) Relaxation process of an unstable gaussian density profile under strong confinement, i.e. ǫ = 0 . | ψ | . One finds a stark contrast in the evolution of thestrongly confined reduction model compared to the unconfined scenario in Fig. 6. In particular, non-diffusive, stable excitationsare ejected from the collapse sight. Their bulk velocity is sufficient to leave the small interaction range of U pconf and thuspropagate freely through the entire domain without re-collapsing to the domain center. The subsequent evolution may thenbe summarized as series of inelastic solitary wave interactions involving different mass ratios and during which high massexcitations slowly consume small mass excitations until a final solitary wave persists at t = 1000. As a byproduct a completelydelocalized fluctuation background emerges that increases in magnitude up to O (1). Panel (A)-(C): Spatially re-centeredHusimi distributions, see Eq. (5). σ x chosen as in Fig. 6. Panel (A): Contrary to (1 + 1)-SP no phase space spiral develops.Instead circular, solitary excitations manage to leave the central gravitational potential. Panel (B): Example of a excitationmerger taking place in a spatially delocalized, fluctuating background. Panel (C): Quasi-stationary, final state after all solitarywaves merged into a single high mass stationary configuration. Panel (E): Energy evolution. Again total energy conservationis apparent. Variations in the potential and kinetic energy originate from solitary interactions and the associate expel ofexcess energy during their mergers. Panel (F): Deviation from the virial theorem in Eq. (25). Note without the boundaryterm the system would depart from virial equilibrium. All energies in (E)/(F) are normalized to the ground state shown inFig. 9. Panel (G): Entropy evolution cf. Eq. (27). Thermalization takes significantly longer compared to (1 + 1)-SP and is onlycompleted at t ≈ 900 when all solitary excitations have been consumed. The thermalization time is again comparable to thevirialization time. O (1) until allstationary states have been consumed by a single highmass solitary wave, see Fig. 26(C). At this point theasymptotic, relaxed system configuration is reached.Inspection of the virialization theorem, Fig. 7(F), andthe entropy evolution, Fig. 7(G), reveal that both ob-servables capture the relaxation process equally well andreport virialization or thermalization around t ≈ 950 re-spectively. Importantly, since the dynamics spans overthe entire domain, the boundary term in Eq. (26) cannotbe neglected. In fact, a naive application of the quantumvirial theorem omitting the boundary term would suggesta departure from the equilibrium state.Let us close this section by mentioning two imperfec-tions of the reported data. Firstly, | ψ | experiences anunphysical symmetry breaking around x = L / t = 100 being induced by small numericalerrors. The same symmetry breaking is also apparent inFig. 6 for (1 + 1)-SP. Assessing the situation in moredetail reveals an absolute momentum drift of h p i = 10 − .Given the long integration time and the high degree ofmobility seen in | ψ | , we find this momentum conserva-tion violation still to be acceptable.Secondly, the evolution of the virial theorem expe-riences unphysical jumps when significant amounts ofmatter travels across the domain boundary, e.g. at t ≈ ψ around these eventsmediates this problem. However, one cannot anticipatea priori when matter flows past the periodic boundary,making the non-stationary virial theorem cumbersometo work with in practice. 3. Final States a. (1 + 1) Schr¨odinger-Poisson Past t ≈ 5, one findsthe system in a quasi-stationary configuration. We cointhe attained state quasi-stationary since the entire matterdensity still undergoes significant time-dependent vari-ation yet does not produce qualitatively new featuresin the spatio-temporal evolution of | ψ | cf. Fig. 6(D). What remains is a circular phase space distribution com-prised of a high density core and halo of fluctuations sur-rounding it, see Fig. 6(B). One may think of this phasespace distribution as the result of smoothing out the fine-grained filament structure of an ever tighter spiralizeddistribution on the scale of Heisenberg’s uncertainty prin-ciple σ x σ k = .To answer whether the core coincides with a (1 + 1)-SPground state, we consider the mean spectral compositionof | ψ k | obtained by averaging 100 quasi-stationary wavefunctions past t > M = L . Theground states were generated independently, as discussedin Sec. II B 4. Evidently, no convincing agreement isachieved and the time asymptotic spectrum appears gen-erally too broad for any viable ground state with M < L .Changing the width and location of the averaging timeinterval does not yield any improvements.We conclude although (1 + 1)-SP realizes a cored den-sity profile, this core is not a ground state configurationof its Hamiltonian — a qualitative difference to (3+1)-SP.Nevertheless, there is still something to be learned aboutthe obtained long term density distribution. The comple-mentary, classical view point of Sec. II A 2 b suggests tocompare the results of (1+1)-SP with predictions for one-dimensional, collisionless N -body systems, in particulardensity profiles for dark matter haloes.For the situation at hand [55] observed how phase mix-ing and violent relaxation drives such system towardspower-law densities ρ ( x ) ∝ | x | − γ with γ ≃ . 5. In-spired by Einasto’s profile, [66], the authors of [67] ex-tended this halo model by an exponential suppressionfactor dominant past a cut-off radius r . Following thisargumentation, one expects: ρ ( x ) ∝ | x | − γ exp − (cid:18) | x | r (cid:19) − γ ! (31)for a d = 1 dimensional CDM halo located at x = 0.In accordance with [39], we consider the integrated andnormalized halo mass M ( | x | ) /L = L R | x | d x ρ ( x ) in-stead of ρ ( x ), thus sparing us to choose a particularvalue of the smoothing scale σ x — recall the necessityof smoothing to obtain a satisfactory Vlasov-Schr¨odingercorrespondence. Figure 8(B) depicts how this halo modelcompares to the simulated (1 + 1)-SP density at t = 1000.We observe a satisfactory correspondence with the fitmodel of Eq. (31) at γ ≈ . 62. Better results are achiev-able in the limit ~ → 0, see [39]. b. (1 + 1) Strong Harmonic Confinement Repeat-ing the spectral analysis in Fig. 9(A) yields convincingaccordance between the mean wave function spectrumand a strongly confined ground state of mass M = 48.The remaining spectral disturbances which were not com-pletely suppressed by the time-averaging are confined to | k | < 10. Comparison with the quasi-stationary phasespace distribution in Fig. 7(C) shows they originate from7 − − | x | / ( L / × − × − M ( | x | ) / L integrated densityprofile (B) − k − − − − | ˆ ψ k | M = M = L (1 + 1)-SP ground statesmean spectrum (A) (Eq. (31)) FIG. 8. (Color online) Analysis of the quasi-stationary stateof (1 + 1)-SP extracted from Fig. 6. Panel (A): Comparisonof the mean wave function spectrum obtained by averaging100 wave functions past t = 990 alongside two (1 + 1)-SPground state spectra. The latter were generated by meansof a gradient descent, Sec. II B 4. The poor correspondencebetween the spectra suggests that the minimal energy solu-tions to (1 + 1)-SP does not act as dynamical attractors inthe evolution. Panel (B): Integrated and normalized FDMdensity obtained from | ψ | at t = 1000 after radial averaging.One finds good correspondence with the halo density modelof Eq. (31) with γ ≈ . the delocalized background oscillations.We conclude that the inelastic excitation dynamics ex-perienced during relaxation does in fact drive the systemtowards a single high mass ground state. One may regardthe minimal energy solution as the fixed point in the longterm evolution under strong confinement — a result instark contrast the our observations for (1 + 1)-SP butqualitatively close to the (3 + 1)-FDM phenomenology.The halo matter undergoes spatial variation which di-minish after averaging multiple radial density profiles inthe same time window used for Fig. 9(A). The resultingmean density of Fig. 9(B) then indicates a power law haloprofile outside the core consistent with the decay behav-ior of the numerically obtained (1+1)-SP halo. This is tobe expected, since there is no reason to assume a stronglyconfined, one dimensional halo organizes into a canonical1 / | x | -NFW profile [30]. Intuitively, the behavior of adark matter halo should be influenced by (i) the interac- − | x | − − | ψ | ∝ | x | − . ground state ( M = 48)mean radial density (B) − k − − − − − | ˆ ψ k | ground state spectrum ( M = 48)mean spectrum (A) FIG. 9. (Color online) Analysis of the quasi-stationary stateunder strong harmonic confinement with ǫ = 0 . 01 extractedfrom Fig. 6. Panel (A): Comparison of the mean wave func-tion spectrum obtained by averaging 100 wave functions past t = 990 alongside the best matching ground state of mass M = 48 cf. Fig. 2. Also note that the remaining deviationsfrom the ground state spectrum are confined to | k | < 10— exactly the wave number regime in which the delocalizedbackground density is situated in Fig. 7. Panel (B): Radi-ally and temporally averaged density. Time averaging as in(A). Again a clear distinction between ground state core andadjacent halo density can be made. Moreover, the halo sur-rounding the core indicates the same power law scaling foundfor the (1 + 1)-SP halo of Fig. 8. tion in the far field and (ii) dimension dependent effectssuch as geometrical dilution.To substantiate this claim, we impose spherical symme-try to the full-fledged problem in Eq. (1) and adopt theregularization approach of [68]. The choice of sphericalsymmetry presents yet another reduction of dimension-ality to (1 + 1). Doing so allows us to implement thesame 1 /r -far field behavior as in our confinement modelwhile making the density dilution in radial direction man-ifest. Excess matter radiated away during the relaxationprocess is suppressed with a complex absorbing potentialsituated at the domain boundary, e.g. [25, 27]. We referto Fig. 10 for the relaxed, mean density profile obtainedfrom a ensemble of N = 20 gaussian initial conditionsof various masses. Application of the scaling symmetry,8 − r − − − | ψ | (3+1)-SP ground stateensemble avg. ρ ( r ) ∝ r − FIG. 10. (Color online) Ensemble averaged radial density de-duced from N = 20 realisations of different mass, initial gaus-sians. The integration was performed under the assumptionof spherical symmetry. Notice that in addition to the groundstate (’solitonic’) core, the remaining matter organizes into apower law halo resembling the behavior of a NFW profile atlarge radii. Eq. (20), allows us to rescale each realization to a com-mon peak density. Evidently, matter not included in theground state core now organizes into a NFW profile. 4. Control of Self-Organization Processes The foregoing discussion highlighted the superiority ofthe transversal confinement model in mapping the (3+1)-FDM phenomenology to one dimensional analogues —while both models indicate reasonable accordance withclassical predictions for the outer halo density, it is onlyunder confinement that the equilibrated state evolves to-wards a ground state core.In fact, the observed self-organization principle of(1 + 1)-FDM under strong confinement is not new. Theauthor of [69] showed how for a class of focusing, lo-cal nonlinearities of the NLSE perturbed uniform initialconditions have a single soliton as dynamical attractor.More precisely, the perturbed initial conditions develop anumber of small mass solitons which subsequently mergeinto a single high-mass soliton at late times. This phe-nomenon was coined soliton turbulence and it was arguedit is ”thermodynamically favorable” for the system to de-velop in this particular way. The authors of [70] later putthese findings on more theoretical grounds by developinga statistical theory around a mean-field approximation ofthe nonlinear Hamiltonian obeying a maximum entropyprinciple.The problem of non-local interactions was consideredin the context of nonlinear optics by [5]. Numericaland analytical arguments showed that the dynamics ismainly driven by the ratio between the interaction range R and the soliton size σ : If the interaction range istoo large, matter far away from a potential soliton, butstill within interaction range, contributes significantly to the convolution integral. Consequently, the delicate po-tential required to form a soliton gets averaged out bythe surrounding fluctuations. Hence, one expect soliton-turbulence-like behavior for R ≪ σ . In case of R ≥ σ ,the system organized into a ”spatially localized incoher-ent structure” coined incoherent soliton . Their resultsresemble our findings for the quasi-stationary state of(1 + 1)-SP.A limit not yet discussed, is the weak confinementregime, i.e. ǫ → ∞ . It is intuitively clear that in thiscase the interaction kernel U pconf should approach G p △ .A more careful analysis shows U pconf ( x , x ′ ) ∼ πǫ G p △ ( x , x ′ ) ( ǫ → ∞ ) . Increasing ǫ should therefore allow us to observe a tran-sition from the soliton to incoherent soliton turbulenceregime. To keep relaxation times comparable we alsosubstitute a → πǫ a so that the effective nonlinear cou-pling stays unity.Figure 11 compares the asymptotic state obtained un-der strong, weak, and no confinement alongside the re-spective interaction range R ( L ), Eq. (16), and solitonextent σ ( M ), Eq. (24).We find confinement parameters larger than unity toquickly approach quasi-stationary states comprised ofmany density maxima beating against each other aroundthe origin. These configurations are qualitatively iden-tical to the (1 + 1)-SP case. Comparing the maximal,boundary condition compatible soliton size σ max withthe interaction range R ( L ) at the chosen domain size L ≈ 127 shows σ/R < . 1, which according to [5] im-plies ”incoherent soliton” dynamics, and in particular nosolitonic attractor. On the other hand comparing bothlength scales for ǫ = 0 . 01, where a M = 48 soliton isformed, we have σ/R > 1, consistent with the solitonturbulence regime. 5. Space-time Expansion The foregoing results of Sec. IV B 4 allow us to ex-tend the discussion to non-static background cosmologies.We first note that the interaction range, as defined inEq. (16), is independent of the nonlinear coupling con-stant. The soliton size σ ( M ), on the other hand, is.This is intuitively clear: Decreasing the nonlinear cou-pling increases the importance of the diffusive characterof kinetic term in the Hamiltonian — we approach a freeSchr¨odinger equation. Hence, the radius at which thefocusing nature of the non-linearity balances the kineticterm is expected to increase as well.Under strong confinement cf. Fig. 11(B) an expand-ing background cosmology would therefore drive the sys-tem even further into the soliton turbulence regime σ ≫ R ( L ).Without confinement cf. Fig. 11(F) the typical in-crease of σ ( M ) experienced by starting from reasonable9 50 60 70 x | ψ | (A) 50 60 70 x (C) 50 60 70 x (E) M or L − l e n g t h (B) R ( L ) σ ( M ) 10 M or L (D) R ( L ) σ ( M ) 10 M or L (F) R ( L ) σ ( M ) ǫ = 0 . ǫ = 5 ǫ = ∞ (SP) FIG. 11. (Color online) Overview of asymptotic states as a function of the confinement parameter ǫ . From left to right, wehave (A)-(B) analyzing the strong confinement limit ( ǫ = 10 − ), (C)-(D) assessing the weak confinement regime ( ǫ = 5) and(E)-(F) evaluating the uniform reduction scenario in the limit ǫ → ∞ , i.e. (1 + 1)-SP. Upper panels: Snapshot of the attainedquasi-stationary states. Lower panels: Comparison of the interaction range R ( L ), Eq. (16), and the soliton size σ ( M ), Eq. (24).Crosses denote soliton sizes directly inferred from the gradient descent of various mass solitons cf. Fig. 4. Solitons with sizesinside the gray shaded area do not exist as they violate periodic boundary conditions. For ǫ > 1, we find a transition awayfrom turbulent soliton dynamics toward a ”incoherent soliton” configuration, i.e. a highly fluctuating state comprised of manydensity maxima beating against each other in real space. As argued in [5], this regime is entered once the interaction range R ( L ) is significantly larger than the soliton size σ ( M ). Even in the best case scenario for (D) and (F), i.e. when a solitonof maximal size could form, one still finds σ max /R ( L ≈ < . σ ( M = 48) /R ( L ≈ > initial redshifts, say z = 100, is insufficient to realize σ max ≈ R ( L ), around which a transition to the soli-ton turbulence regime should occur. Note that here σ max > R ( L ) is never achievable as it would violate theperiodic boundary conditions. That said, allowing for atime dependent coupling constant has, aside from numer-ical implications, also influence on the relaxation time.Preliminary analysis shows that although the strong con-finement scenario including a FLRW background doestrend towards a ground state (’solitonic’) spectrum, re-laxation is not completed at z = 0. Figure 12 illustratesthis result. A full-fledged analysis of the expanding modelis left to future work. V. CONCLUSION Purpose of this work was to conduct an extensive nu-merical study on the applicability of the Fuzzy Dark mat-ter (FDM) model in one spatial dimension. Particularemphasis was put on (i) properties of system’s long termevolution, (ii) the structure of the relaxed, asymptoticsystem state and how it compares to the core-halo struc- ture of (3 + 1)-FDM, as well as (iii) which model pa-rameters may be used to control its phenomenology. Tothis end, we derived two distinct one-dimensional FDMmodels by either allowing for a complete delocalizationof matter in the transversal plane, (1 + 1)-SP, or by con-fining all matter along one spatial direction. While bothmodels realize long range interactions free of singularities,it is only under strong confinement that the nonlocal in-teraction recovers the desired − /r interaction at largedistances.We proceeded to investigate the mean cosmic struc-ture growth in an ensemble of FDM-only, flat FLRW-universes obeying (1 + 1)-SP and starting from cosmo-logical initial conditions. By following the evolution ofthe matter power spectrum until present time, two dis-tinct spectral ranges were identified: Firstly, a suppres-sion range, in which the power spectrum is smoothedout by the uncertainty principle. Secondly, a couplingregime, where the redistribution of matter across nonlin-early evolving modes leads to a scale-free spectrum con-sistent with considerations for (1 + 1)-CDM. The tran-sition scale from coupling to suppression followed by aself-consistent application of the uncertainty principle.0 − k − − − − − | ˆ ψ k | ground state spectrum ( M = 48)spectrum at z=0 (A) z S (B) FIG. 12. (Color online) Strong confinement model with ǫ =0 . 01 undergoing collapse in a background cosmology as inSec. IV A 1. Panel (A): Final wave function spectrum at z =0 together with the same M = 48 ground state of Fig. 9.Panel (B): Entropy evolution. It is evident that relaxation isnot completed at z = 0. The analysis of the asymptotic system state was con-ducted under simplifying assumptions, i.e., a static back-ground cosmology and spatially localized initial condi-tions. Our study suggests that (1+1)-SP relies on violentrelaxation and phase mixing to approach its equilibrated,i.e. thermalized and virialized, system state. Althoughthe realized asymptotic configuration does form a core-halo structure, the central core cannot be identified withthe ground state solution of the (1 + 1)-SP Hamiltonian— a stark contrast to (3+1)-FDM for which the minimumenergy solution acts as dynamical attractor in the longterm evolution. The halo density, on the other hand, wasfound to be consistent with structure of one dimensionalCDM halos.By contrast, the evolution under strong confinementfavours the (3 + 1)-FDM typical relaxation mechanism ofgravitational cooling. Ultimately, the evolution then con-verges into a virialized and thermalized system state com-prised of a single high-mass ground state solution embed-ded in a delocalized fluctuation background that emergesfrom a series of inelastic ground state interactions. Theanalysis of the halo density suggested identical CDM-likebehavior as for (1 + 1)-SP. We conclude that under thechosen simulation conditions the strongly confined reduc-tion model is superior in mapping the three-dimensionalphenomenology to one spatial dimension.To understand the reason for the qualitative differencebetween the asymptotic behavior of (1 + 1)-SP and theconfinement model, we investigated the weak confine-ment limit of our reduction. In accordance with argu-ments from nonlinear optics, see e.g. [5], we found the system to converge towards a high-mass ground state ifthe effective interaction range is (considerably) smallerthan the spatial extent of the ground state. For (1 + 1)-SP and arbitrary but fixed coupling constant, no groundstate exists that satisfies this condition. It is for this rea-son, that we conjectured our static space-time analysisremains valid even for nontrivial background cosmologies.Although our work focused on an interaction kernel thatresembles Newtonian gravity, results from nonlinear op-tics suggest that the dependence of the asymptotic stateon the interaction range is a property also applicable forother long range potentials. For instance, the authors of[5] implemented a gaussian interaction kernel in Eq. (12),while [6] employed a screened Poisson equation as fieldequation. The latter approach implies an exponentiallydecaying Green’s function.Our work may be extended in multiple regards. Froma physical perspective, a full-fledged investigation of cos-mological initial conditions in various cosmological ex-pansion models is still pending and the phenomenologyunder strong confinement is presumably not exhaustedby our discussion. In this context, we mention the prop-erties of the delocalized fluctuation background as it maybe possible to understand it as an ensemble of small scaleplane waves obeying a dispersion relation akin to Bogoli-ubov’s excitation spectrum for Bose-Einstein condensates[4, 71].Moreover, additional conceptional optimizations of ourconfinement approach are worth exploring. For instance,our work only focused on a global, statically set con-finement parameter. However, incorporating the con-finement ansatz directly into SP’s generating action isexpected to yield additional information on the spatio-temporal evolution of the confinement strength itself,thereby allowing it to be set self-consistently and depen-dent on the wavefunction evolution [45].We also remark on our ongoing effort to optimize ournumerical approach by means of a more efficient basis-function method or splitting schemes with intrinsic errorestimates. With this we hope to (i) achieve a fully adap-tive spatio-temporal grid that is sensitive to nonlinearevolution and (ii) pave the way for a higher dimensionalanalysis. The latter should allow the investigation of ad-ditional relaxation channels unique to FDM, especiallythe emission of quantized vortices that may play an im-portant role for the asymptotic evolution on top of grav-itational cooling.Altogether, we hope that our investigation will leadto a better cross-fertilisation, see also [11], between cos-mology, statistical mechanics [72, 73], nonlinear dynam-ics [74, 75], nonlinear wave optics [5–8], and the quan-tum evolution of, e.g., Bose-Einstein condensates [76–82],with the scope to obtain further insight into the complexdynamics of the SP model and to search for possible labo-ratory experiments for the implementation and analoguesimulation of FDM models.1 ACKNOWLEDGMENTS The authors acknowledge support by the state ofBaden-W¨urttemberg, Germany through bwHPC. Ourspecial gratitude goes to Luca Amendola for his feedbackto this work, Javier Madro˜nero for various conversationson NLSE dynamics and to Jens Niemeyer for the insight-ful discussion on the subject in general. Appendix A: Convergence and Stability Let us now give a numerical justification for the accu-racy and stability claims mentioned in Sec. III C. To thisend, we conduct a convergence and stability study for thesimulation scenarios of Sec. IV A 1 (delocalized randomfield including space-time expansion) and IV B 1 (unsta-ble gaussian initial conditions in a static background cos-mology).The authors are not aware of a general analytical re-sult in any of these cases. Thus, we compute a referencesolution ψ ref on a fine spatio-temporal grid { N ref , ∆ t ref } and measure the error of ψ relative to ψ ref via ∆ ǫ = k ψ ref − ψ k / k ψ ref k . Note the error ∆ ǫ is a functionof N , ∆ t and integration time t . Since the results arequalitatively identical for (1 + 1)-SP and the confinementmodel we only report data for the former.We begin with the cosmological simulation scenarioof Sec. IV A, i.e. a random field evolving in a dynamicFLRW background. Figure 13 depicts the dependence ofthe numerical error as a function of the spatio-temporalgrid parameters { N, ∆ t } relative to the reference grid∆ t ref = 10 − and N ref = 2 (the white cells). The gridparameters of Sec. IV A are marked with a black cross.To assure comparability, all 56 { N, ∆ t } combinations areinitialized with every N/N ref point of the same referencegaussian random field δ ( x ).In the m = m case, we find a high degree of uniformityin N at fixed ∆ t throughout the integration. This is theresult of the spectral accuracy of the employed spatialdiscretization. Variations of the relative error ∆ ǫ in N at ∆ t = 1 − × − are common once we reach theconvergence plateau but may also be induced by a lack offidelity close to the reference solution. More pronounced is the loss of accuracy in ∆ t direction (at any considered N ) and we conclude the overall inaccuracy is dominatedby the temporal error. Solutions ∆ t ≤ × − can beconsidered as converged.The situation for the m = m scenario is qualitativelyidentical with an additional caveat at low red shifts. Here,spatial grids with N < proof to be insufficient toresolve the entire spectrum of | ˆ ψ k | . In this case spatio-temporal grids with N > and ∆ t ≤ × − aredeemed sufficient to achieve convergence.Concerning temporal accuracy and overall stability, werefer to Fig. 14 which evaluates the time dependence ofthe numerical error for various time increments ∆ t atfixed N = N ref = 2 . The reference solution is identicalto the one used in Fig. 13 and the spatio-temporal gridof Sec. IV A is depicted as red, dashed line.Evidently, allowing for an dynamic cosmological back-ground yields a numerical error that evolves roughly ex-ponentially. Nevertheless, errors are still of acceptablesize at present time which is why we deem our numeri-cal treatment of the non-autonomous Hamiltonian as ac-ceptable for the purposes of this work. Decreasing thetime step beyond ∆ t < − results in no significantgain in accuracy. Note that the non-converged time steps,i.e. ∆ t = 10 − , − approach the convergence plateauwith roughly quadratic speed. This can be inferred fromtheir relative offset which evaluates to about 2 orders ofmagnitude and is expected given that Strang splitting cf.Sec. III B is second order accurate in time.Let us contrast the error evolution of the fully cos-mological case with the static simulation conditions ofSec. IV B cf. Fig. 15. Here, the smaller dimension-less box size allows us to reduce the number of spatialgrid points required to fully resolve the spectrum of ψ .This in turn makes a smaller reference time incrementpossible. Figure 15 therefore uses a reference grid with∆ t ref = 10 − and N ref = 2 . Again, the grid parametersused in Sec. IV B correspond to the red, dashed line.As for the cosmological evolution, quadratic accuracyis achieved for non-converged time steps approaching theplateau. A notable difference, however, is the growthbehavior of ∆ ǫ ( t ) which only evolves linearly in time ifthe scale factor remains static. It is for this reason thatwe can extend the integration time up to t = 1000 cf.Sec. II C without losing reliability of our data. [1] A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems (DoverPublications, 2003).[2] G. D. Mahan, Many Particle Physics (Physics of Solidsand Liquids) (Kluwer Academic/Plenum Publishers, NewYork, 2000).[3] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari,Theory of Bose-Einstein condensation in trapped gases,Reviews of Modern Physics , 463 (1999).[4] L. P. Pitaevskii and S. Stringari, Bose-Einstein Conden- sation and Superfluidity (Oxford University Press, 2016).[5] A. Picozzi and J. Garnier, Incoherent Soliton Turbulencein Nonlocal Nonlinear Media, Physical Review Letters , 10.1103/physrevlett.107.233901 (2011).[6] R. Bekenstein, R. Schley, M. Mutzafi, C. Rotschild,and M. Segev, Optical simulations of gravita-tional effects in the Newton–Schr¨odinger system,Nature Physics , 872 (2015).[7] T. Roger, C. Maitland, K. Wilson, N. Westerberg,D. Vocke, E. M. Wright, and D. Faccio, Optical analogues N z = 90 z = 50 z = 10 z = 2 z = 1 z = 02 N z = 90 2 z = 50 2 ∆ t z = 10 2 z = 2 2 z = 1 2 z = 0 × − − − − − − − − ∆ ǫ FIG. 13. (Color online) Evolution of the relative error ∆ ǫ as a function of the spatio-temporal grid. The convergence studyadopts the initial conditions and parameters of Sec. IV A 1 using m = 10 − eV in the first row and m = 10 − eV in thesecond row. Space-time expansion is turned on. The white cells corresponds to the numerical reference solution ψ ref relativeto which ∆ ǫ is computed. The black cross represents the spatio-temporal grid used for Fig. 5.of the Newton–Schr¨odinger equation and boson star evo-lution, Nature Communications , 13492 (2016).[8] A. Navarrete, A. Paredes, J. R. Salgueiro, andH. Michinel, Spatial solitons in thermo-optical me-dia from the nonlinear Schr¨odinger-Poisson equationand dark-matter analogs, Physical Review A ,10.1103/physreva.95.013844 (2017).[9] L. Di´osi, Gravitation and quantum-mechanical localization of macro-objects,Physics Letters A , 199 (1984).[10] R. Ruffini and S. Bonazzola, Systems of Self-GravitatingParticles in General Relativity and the Concept of anEquation of State, Phys. Rev. , 1767 (1969).[11] A. Paredes Galan, D. Olivieri, and H. Michinel,From optics to dark matter: A reviewon nonlinear Schr¨odinger–Poisson systems,Physica D: Nonlinear Phenomena , 132301 (2019).[12] W. Hu, R. Barkana, and A. Gruzinov, Fuzzy Cold DarkMatter: The Wave Properties of Ultralight Particles,Physical Review Letters , 1158 (2000).[13] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Cosmicstructure as the quantum interference of a coherent darkwave, Nature Physics , 496 (2014).[14] J. S. Bullock and M. Boylan-Kolchin, Small-Scale Challenges to the ΛCDM Paradigm, An-nual Review of Astronomy and Astrophysics10.1146/annurev-astro-091916-055313 (2017).[15] B. Moore, Evidence against dissipation-lessdark matter from observations of galaxy haloes,Nature , 629 (1994).[16] W. J. G. de Blok, The Core-Cusp Problem, Advances inAstronomy , 789293 (2010). [17] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala,M. Boylan-Kolchin, A. Fialkov, and L. Hern-quist, Galaxy formation with BECDM – I.Turbulence and relaxation of idealized haloes,Monthly Notices of the Royal Astronomical Society , 4559 (2017).[18] N. K. Porayko, X. Zhu, Y. Levin, L. Hui, G. Hobbs,A. Grudskaya, K. Postnov, M. Bailes, N. D. R. Bhat,W. Coles, S. Dai, J. Dempsey, M. J. Keith, M. Kerr,M. Kramer, P. D. Lasky, R. N. Manchester, S. Os lowski,A. Parthasarathy, V. Ravi, D. J. Reardon, P. A.Rosado, C. J. Russell, R. M. Shannon, R. Spiewak,W. van Straten, L. Toomey, J. Wang, L. Wen, andX. You (PPTA Collaboration), Parkes Pulsar TimingArray constraints on ultralight scalar-field dark matter,Phys. Rev. D , 102002 (2018).[19] N. C. Amorisco and A. Loeb, First constraints on FuzzyDark Matter from the dynamics of stellar streams in theMilky Way (2018), arXiv:1808.00464 [astro-ph.GA].[20] A. Lidz and L. Hui, The Implications of a Pre-reionization 21 cm Absorption Signal for Fuzzy Dark Mat-ter, Phys. Rev. D , 023011 (2018).[21] J. C. Niemeyer, Small-scale structureof fuzzy and axion-like dark matter,Progress in Particle and Nuclear Physics , 103787 (2020).[22] L. M. Widrow and N. Kaiser, Using the SchroedingerEquation to Simulate Collisionless Matter,The Astrophysical Journal , L71 (1993).[23] C. Uhlemann, M. Kopp, and T. Haugg, Schr¨odingermethod as N-body double and UV completion ofdust, Physical Review D , 10.1103/physrevd.90.023517(2014).[24] M. Kopp, K. Vattis, and C. Skordis, Solving the t − − − − ∆ ǫ z ∆ t = · − ∆ t = · − ∆ t = 1 · − ∆ t = 2 · − (A) t − − − − ∆ ǫ ∆ t = · − ∆ t = · − ∆ t = 1 · − ∆ t = 2 · − (B) FIG. 14. (Color online) Evolution of the numerical error ∆ ǫ as a function of the integration time. The convergence studyadopts the initial conditions and parameters of Sec. IV A 1with m = 10 − eV in panel (A) and m = 10 − eV in panel(B). Space-time expansion a ( t ) is turned on. The space-timegrid of Sec. IV B corresponds to the black, dashed line. t − − − − ∆ ǫ ∆ t = · − ∆ t = · − ∆ t = 1 · − ∆ t = 2 · − FIG. 15. (Color online) Evolution of the numerical error ∆ ǫ as a function of the integration time. The convergence studyadopts the initial conditions and parameters of Sec. IV B 1.Space-time expansion is turned off, i.e., a = 1. The space-time grid of Sec. IV B corresponds to the black, dashed line.Vlasov equation in two spatial dimensions withthe Schr¨odinger method, Physical Review D ,10.1103/physrevd.96.123532 (2017).[25] F. S. Guzm´an and L. A. Ure˜na-L´opez, Evolution of theSchr¨odinger-Newton system for a self-gravitating scalarfield, Physical Review D , 10.1103/physrevd.69.124033(2004). [26] F. S. Guzman and L. A. Urena-Lopez, Gravita-tional Cooling of Self-gravitating Bose Condensates,The Astrophysical Journal , 814 (2006).[27] B. Schwabe, J. C. Niemeyer, and J. F. Engels, Sim-ulations of solitonic core mergers in ultralight ax-ion dark matter cosmologies, Physical Review D ,10.1103/physrevd.94.043513 (2016).[28] F. Edwards, E. Kendall, S. Hotchkiss, andR. Easther, PyUltraLight: A Pseudo-SpectralSolver for Ultralight Dark Matter Dynamics,Journal of Cosmology and Astroparticle Physics (10), 027.[29] E. Seidel and W.-M. Suen, Formation ofsolitonic stars through gravitational cooling,Physical Review Letters , 2516 (1994).[30] J. F. Navarro, C. S. Frenk, and S. D. M. White,The Structure of Cold Dark Matter Halos,The Astrophysical Journal , 563 (1996).[31] M. Modugno, C. Tozzo, and F. Dalfovo, Role of trans-verse excitations in the instability of Bose-Einstein con-densates moving in optical lattices, Physical Review A , 10.1103/physreva.70.043625 (2004).[32] P. H. Chavanis, Growth of perturbations in an expand-ing universe with Bose-Einstein condensate dark matter,Astronomy & Astrophysics , A127 (2012).[33] T.-P. Woo and T. Chiueh, High-Resolution Simulation onStructure Formation with Extremely Light Bosonic DarkMatter, The Astrophysical Journal , 850 (2009).[34] X. Li, L. Hui, and G. L. Bryan, Numerical and perturba-tive computations of the fuzzy dark matter model, Phys-ical Review D , 10.1103/physrevd.99.063509 (2019).[35] E. Madelung, Quantentheorie in hydrodynamischerForm, Zeitschrift f¨ur Physik , 322 (1927).[36] O. D. Kellogg, Foundations of Potential Theory (Springer Berlin Heidelberg, 1967).[37] S. L. Marshall, A periodic Green functionfor calculation of coloumbic lattice potentials,Journal of Physics: Condensed Matter , 4575 (2000).[38] Gradshteyn, I. S. and Ryzhik, I. M., Table of Integrals,Series, and Products (Elsevier LTD, Oxford, 2014).[39] M. Garny and T. Konstandin, Gravitationalcollapse in the Schr¨odinger-Poisson system,Journal of Cosmology and Astroparticle Physics (01), 009.[40] T. Zimmermann, M. Pietroni, J. Madro˜nero,L. Amendola, and S. Wimberger, A QuantumModel for the Dynamics of Cold Dark Matter,Condensed Matter , 89 (2019).[41] M. Garny, T. Konstandin, and H. Rubira, TheSchr¨odinger-Poisson method for Large-Scale Structure,Journal of Cosmology and Astroparticle Physics (04), 003.[42] W. Bao, H. Jian, N. J. Mauser, and Y. Zhang,Dimension Reduction of the Schr¨odinger Equationwith Coulomb and Anisotropic Confining Potentials,SIAM Journal on Applied Mathematics , 2100 (2013).[43] I. Bloch, J. Dalibard, and W. Zwerger,Many-Body Physics with Ultracold Gases,Rev. Mod. Phys. , 885 (2008).[44] M. Olshanii, Atomic Scattering in the Presence of an Ex-ternal Confinement and a Gas of Impenetrable Bosons,Phys. Rev. Lett. , 938 (1998).[45] L. Salasnich, A. Parola, and L. Reatto, Ef-fective wave equations for the dynamics ofcigar-shaped and disk-shaped Bose condensates,Phys. Rev. A , 043614 (2002).[46] S. Wimberger, R. Mannella, O. Morsch, and E. Arimondo, Resonant Nonlinear Quantum Trans-port for a Periodically Kicked Bose Condensate,Phys. Rev. Lett. , 130404 (2005).[47] S. Wimberger, R. Mannella, O. Morsch, E. Arimondo,A. R. Kolovsky, and A. Buchleitner, Nonlinearity-induced destruction of resonant tunneling in the Wannier-Stark problem, Phys. Rev. A , 063610 (2005).[48] D. B. Owen, M. Abramowitz, and I. A. Stegun, Hand-book of Mathematical Functions with Formulas, Graphs,and Mathematical Tables, Technometrics , 78 (1965).[49] W. Bao and Q. Du, Computing the GroundState Solution of Bose–Einstein Conden-sates by a Normalized Gradient Flow,SIAM Journal on Scientific Computing , 1674 (2004).[50] W. Bao, I.-L. Chern, and F. Y. Lim, Efficient and spec-trally accurate numerical methods for computing groundand first excited states in Bose–Einstein condensates,Journal of Computational Physics , 836 (2006).[51] P. Choquard, J. Stubbe, and M. Vuffray, Stationarysolutions of the Schr¨odinger-Newton model—an ODEapproach, Differential and integral equations , 665(2008).[52] T. Pang, An Introduction to Computational Physics (Cambridge University Press, 2006).[53] P. G. Drazin, Solitons: an introduction (Cambridge Uni-versity Press, Cambridge England New York, 1989).[54] D. Lynden-Bell, Statistical Mechanics ofViolent Relaxation in Stellar Systems,Monthly Notices of the Royal Astronomical Society , 101 (1967).[55] J. Binney, Discreteness effects incosmological N-body simulations,Monthly Notices of the Royal Astronomical Society , 939 (2004).[56] E. Weislinger and G. Olivier, The classi-cal and quantum mechanical virial theorem,International Journal of Quantum Chemistry , 389 (2009).[57] J. G. Esteve, F. Falceto, and P. R. Giri, Boundary con-tributions to the hypervirial theorem, Physical Review A , 10.1103/physreva.85.022104 (2012).[58] A. Wehrl, On the relation between clas-sical and quantum-mechanical entropy,Reports on Mathematical Physics , 353 (1979).[59] X. Antoine, W. Bao, and C. Besse, Computa-tional methods for the dynamics of the non-linear Schr¨odinger/Gross-Pitaevskii equations,Computer Physics Communications , 2621 (2013).[60] J. Zhang, H. Liu, and M.-C. Chu, Cosmological Simu-lation for Fuzzy Dark Matter Model, Frontiers in As-tronomy and Space Sciences , 10.3389/fspas.2018.00048(2019).[61] F. C. Sergio Blanes, A Concise Introduction to Geomet-ric Numerical Integration (Apple Academic Press Inc.,2016).[62] S. Dodelson, Modern Cosmology (Elsevier LTD, Oxford,2003).[63] A. Lewis, A. Challinor, and A. Lasenby, Ef-ficient computation of CMB anisotropies inclosed FRW models, ApJ , 473 (2000),arXiv:astro-ph/9911177 [astro-ph].[64] S.-F. Chen and M. Pietroni, Asymp-totic expansions for Large Scale Structure, Journal of Cosmology and Astroparticle Physics (06), 033.[65] M. Pietroni, Structure formation beyond shell-crossing:nonperturbative expansions and late-time attractors,Journal of Cosmology and Astroparticle Physics (06), 028.[66] J. Einasto, On the Construction of a Composite Modelfor the Galaxy and on the Determination of the System ofGalactic Parameters, Trudy Astrofizicheskogo InstitutaAlma-Ata , 87 (1965).[67] A. E. Schulz, W. Dehnen, G. Jungman, andS. Tremaine, Gravitational Collapse in One Dimension,Monthly Notices of the Royal Astronomical Society , 49 (2013).[68] X. Dong, A short note on simplified pseudospectral meth-ods for computing ground state and dynamics of spher-ically symmetric Schr¨odinger-Poisson-Slater system,Journal of Computational Physics , 7917 (2011).[69] V. Zakharov, A. Pushkarev, V. Shvets, and V. Yan’kov,Soliton turbulence, JETP Lett , 79 (1988).[70] R. Jordan and C. Josserand, Self-organization in nonlin-ear wave turbulence, Phys. Rev. E , 1527 (2000).[71] P. Maddaloni, M. Modugno, C. Fort, F. Mi-nardi, and M. Inguscio, Collective Oscillationsof Two Colliding Bose-Einstein Condensates,Physical Review Letters , 2413 (2000).[72] A. N. Kolmogorov, The local structureof turbulence in incompressible viscousfluid for very large Reynolds numbers,Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences , 9 (1991).[73] M. Kobayashi and M. Tsubota, Kolmogorov Spectrum ofSuperfluid Turbulence: Numerical Analysis of the Gross-Pitaevskii Equation with a Small-Scale Dissipation, Phys-ical Review Letters , 10.1103/physrevlett.94.065302(2005).[74] M. Tabor, Chaos and Integrability in Nonlinear Dynam-ics (John Wiley & Sons, New York, 1989).[75] M. Lakshmanan and S. Rajaseekar, Nonlinear Dynam-ics: Integrability, Chaos and Patterns (Springer Verlag,Heidelberg, 2003).[76] P. Jain, S. Weinfurtner, M. Visser, and C. W. Gardiner,Analogue model of a FRW universe in Bose-Einsteincondensates: Application of the classical field method,Phys. Rev. A , 033616 (2007).[77] F. Girelli, S. Liberati, and L. Sindoni, Gravi-tational dynamics in Bose Einstein condensates,Phys. Rev. D , 084013 (2008).[78] D. Proment, S. Nazarenko, and M. Onorato, Quan-tum turbulence cascades in the Gross-Pitaevskii model,Phys. Rev. A , 051603 (2009).[79] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, andT. Pfau, The physics of dipolar bosonic quantum gases,Reports on Progress in Physics , 126401 (2009).[80] R. Plestid, P. Mahon, and D. H. J. O’Dell, Violent re-laxation in quantum fluids with long-range interactions,Phys. Rev. E , 012112 (2018).[81] M. Combescot, R. Combescot, and F. Dubin,Bose–Einstein condensation and indirect excitons: a re-view, Reports on Progress in Physics , 066501 (2017).[82] L. Berezhiani and J. Khoury, Emergent long-range interactions in Bose-Einstein condensates,Phys. Rev. D99