PyUltraLight: A Pseudo-Spectral Solver for Ultralight Dark Matter Dynamics
Faber Edwards, Emily Kendall, Shaun Hotchkiss, Richard Easther
PPrepared for submission to JCAP
PyUltraLight: A Pseudo-SpectralSolver for Ultralight Dark MatterDynamics
Faber Edwards, Emily Kendall, Shaun Hotchkiss, and RichardEasther
Department of Physics, University of Auckland, Private Bag 92019, Auckland, New ZealandE-mail: [email protected], [email protected],[email protected], [email protected]
Abstract.
PyUltraLight simulates the dynamics of ultralight dark matter in a non-expanding background.
PyUltraLight can describe the evolution of several interactingultralight dark matter halos or one or more halos orbiting a central, fixed Newtonian po-tential, the latter scenario corresponding to dwarf galaxies orbiting a massive central galaxy.We verify
PyUltraLight by showing that it reproduces qualitative dynamical features ofpreviously published simulations and demonstrate that it has excellent energy-conservationproperties.
PyUltraLight is implemented in a Python-based Jupyter notebook, solving theSchrödinger-Poisson equation governing ultralight scalar field dark matter dynamics in thenon-relativistic regime using a symmetrised split-step pseudospectral algorithm. The note-book interface makes it simple to specify simulation parameters and visualise the resultingoutput but performance-critical routines are managed via calls to computationally efficientcompiled libraries.
PyUltraLight runs on standard desktop hardware with support forshared memory mutlithreading and is available on GitHub. a r X i v : . [ a s t r o - ph . C O ] O c t ontents PyUltraLight The realisation that there may be more to the universe than meets the eye is one of the mostprofound developments in 20 th Century astronomy and astrophysics. However, while thereare now multiple lines of evidence that dark matter outweighs baryonic matter by a ratioof approximately 5:1 [1], we have few clues regarding the physical nature of dark matter.Much theoretical and experimental effort has focused on WIMP models, motivated by theirconsistency with supersymmetric extensions to the Standard Model and their relatively sim-ple dynamics. However, advanced direct-detection experiments are putting increasingly tightconstraints on the WIMP parameter space [2, 3] and Λ CDM cosmology with simple, pres-sureless, noninteracting dark matter (a class including simple WIMP scenarios) is potentiallyat odds with observations at small astrophysical scales [4].The potential shortcomings of simple cold dark matter scenarios motivate investigationsof more novel dark matter scenarios. In particular, ultralight dark matter (ULDM), alsoknown as fuzzy dark matter (FDM), or BEC dark matter, is an increasingly well-studied pos-sibility; for a recent review of the potential advantages and characteristic attributes of thisscenario see Ref [5]. ULDM models are well motivated by fundamental theories possessingapproximate shift symmetries such as the theory of the QCD axion [6, 7]. Moreover, ULDMcan naturally resolve the small-scale problems of Λ CDM as the Heisenberg uncertainty princi-ple suppresses gravitational collapse on length scales shorter than the de Broglie wavelength– 1 –f the ULDM particle. In this regime the mass of the ULDM particle becomes correlatedwith astrophysical observables; if it is on the order of − eV, structure is suppressed atkiloparsec scales and below [8].Given the presence of a fundamental lengthscale, the behaviour of ULDM is more com-plex than that of simple dark matter scenarios whose cosmologically relevant interactionsare purely gravitational. Physically, the effective short-scale pressure and condensate-likeproperties of ULDM create new dynamical possibilities for ULDM scenarios, such as purelypressure supported soliton-like solutions [9] and superposition or interference during interac-tions between condensate-like halos [10]. Consequently, modelling dark matter dynamics inULDM scenarios is more challenging than in simple cold dark matter models, but is criticalto understanding the physical consequences of ULDM models.In the non-relativistic regime, the dynamics of ULDM can be reduced to the Schrödinger-Poisson system, where the complex variable ψ describes the local density of ULDM quantawhile the Poisson equation describes the local gravitational potential. Many approaches havebeen taken to this problem, including both modifications of existing cosmological simulationcodes and the development of new codes specifically designed for ULDM systems. One widelyused approach is the Madelung fluid formulation of the Schrödinger-Poisson system [11] whichhas a quantum pressure term that can be treated numerically in a variety of ways. In Ref. [12],the cosmological code gadget [13] is modified to treat the quantum pressure as an effectiveparticle-particle interaction and the resulting code, axion-gadget is publicly available [14].Ref. [15] modifies a non-public extension of gadget , p-gadget3 to treat the quantumpressure term via smoothed-particle hydrodynamics (SPH) routines. The SPH approach isalso used in Ref. [16], while a particle-mesh approach was implemented in [17]. N yx [18] wasmodified in [10] to study merging ULDM solitonic cores, G alacticus [19] was modified in[20] to study the effects of tidal stripping and dynamical friction on ULDM halos, arepo [21] was modified in [22] to study the core-mass relationship and turbulence characteristics ofULDM halos, and gamer [23, 24] was modified [25] to perform a detailed study of structureformation in ULDM cosmologies.While a large number of public codes can solve conventional dark matter scenarios, axion-gadget is the only currently available solver for ULDM dynamics. This paper in-troduces PyUltraLight , a stand-alone Python-based pseudospectral Schrödinger-Poissonsolver, and demonstrates that it reproduces many of the key findings of more complicatedcosmological simulation codes within a desktop computing environment. We anticipate thatas a publicly available resource,
PyUltraLight will serve as a valuable cross-check on morecomplex implementations, serve as a basis for further development of such codes within thecomputational cosmology community, and facilitate explorations of ULDM dynamics.
PyUltraLight is based on a symmetrised-split-step (leapfrog) solver for the time evo-lution, and uses a pseudospectral Fourier algorithm to solve the Poisson equation for the grav-itational potential at each step. This algorithm has nd order accurate time integration stepsand sub-percent level energy conservation, while the wavefunction normalisation is conservedto machine precision. As a pseudospectral code, linear differential operators are computed bydirect multiplication in the Fourier domain, while non-linear terms are evaluated in positionspace. Consequently, PyUltraLight is free from noise associated with spatial derivatives A similar methodology was described in Ref. [26]; at the time of writing this code has not been released.Spectral methods are often used to solve the Poisson equation in large scale structure simulations, while the
PSpectre code [27] provides a pseudospectral solver for the evolution of the inflaton and fields coupled to itduring parametric resonance and preheating after inflation [28–30]. – 2 –omputed via finite-differencing. There is a necessary computational cost associated with theFourier and inverse Fourier transforms but these transforms are optimised in
PyUltraLight through the use of the pyFFTW pythonic wrapper around the C-based FFTW subroutinelibrary [31, 32]. As the FFTW libraries offer full parallelisation,
PyUltraLight is currentlydesigned to take advantage of multiple cores on a user PC or shared-memory environment.Full MPI compatibility has not yet been implemented as we have not found a need to run sim-ulations in a distributed-memory cluster environment, however future releases may addressthis possibility.This paper is organised as follows. We first provide a short review of ULDM physics,including a derivation of the Schrödinger-Poisson equations from the underlying scalar-fieldLagrangian. We then describe their implementation in
PyUltraLight , before moving on todescribe testing and verification procedures applied to the code. We reproduce a selection ofresults from a variety of recent ULDM simulations and discuss the energy conservation andaccuracy as a function of spatial resolution.
The existence of an extremely light scalar field, minimally coupled to gravity, is the centralpremise on which ULDM models are predicated. Within the ULDM framework, it is proposedthat this scalar field exists as a Bose-Einstein condensate, described by a single wavefunctionwhich is governed by the Schrödinger-Poisson coupled differential equations. We begin byderiving this system of equations as a non-relativistic weak-field limit of a more generaltheory. We start from the action functional for a scalar field, φ , minimally coupled to gravityand in the absence of self-interactions, S = (cid:90) d x (cid:126) √− g (cid:26) g µν ∂ µ φ∂ ν φ − m (cid:126) φ (cid:27) , (2.1)where we have taken c = 1 but retain factors of (cid:126) at this stage. Applying the variationalprinciple to this action yields the Euler-Lagrange equations √− g ∂ µ (cid:2) √− gg µν ∂ ν φ (cid:3) − m (cid:126) φ = 0 . (2.2)We evaluate equation 2.2 using linear perturbation theory, adopting the perturbed FRWmetric in the Newtonian gauge: ds = − (cid:0) (cid:126)r, t ) (cid:1) dt + a ( t ) (cid:0) − (cid:126)r, t ) (cid:1) d(cid:126)r . (2.3)To linear order in Φ( (cid:126)r ) we obtain ¨ φ − (cid:0) (cid:1) a ( t ) ∇ φ + 3 H ˙ φ − φ + (cid:0) (cid:1) m (cid:126) φ = 0 , (2.4)where H = ˙ a ( t ) /a ( t ) . At late times in an expanding universe, H (cid:28) m/ (cid:126) and it is sufficient toset H = 0 and a ( t ) = 1 in equation 2.4. This is a good approximation even at relatively highredshifts, including the epochs of early structure formation. Alternatively, if we consider anon-expanding universe, these equalities are true by definition. In either case, we can remove– 3 –he third term in equation 2.4. The resulting equation can then be analysed using WKBmethods in the non-relativistic regime. This allows us to write an ansatz solution for thefield φ : φ = (cid:126) √ m (cid:0) ψe − imt/ (cid:126) + ψ ∗ e imt/ (cid:126) (cid:1) , (2.5)where ψ is assumed to be slowly varying in the sense that m | ψ | (cid:29) (cid:126) | ˙ ψ | , m | ˙ ψ | (cid:29) (cid:126) | ¨ ψ | , m | ψ | (cid:29) (cid:126) |∇ ψ | , and m |∇ ψ | (cid:29) (cid:126) |∇ ψ | . Since Φ is sourced by ψ , we also have that m | Φ | (cid:29) (cid:126) | ˙Φ | .Direct substitution of the ansatz solution into equation 2.4, discarding heavily suppressedterms, yields i (cid:126) ˙ ψ = − (cid:126) m ∇ ψ + m Φ ψ. (2.6)We have thus shown that ψ satisfies the Schrödinger equation in this limit, which is interpretedas the macroscopic wavefunction of a Bose-Einstein condensate. It follows that the particlenumber density of the condensate is given by | ψ | , so its mass density is simply m | ψ | . Thelocal gravitational potential thus satisfies the Poisson equation, ∇ Φ = 4 πGm | ψ | , (2.7)where G is Newton’s gravitational constant. The coupled equations 2.6 and 2.7 together formthe nonlinear Schrödinger-Poisson system which describes the dynamics of ULDM in thenon-relativistic regime. While Equations 2.6 and 2.7 are valid for open boundary conditions, PyUltraLight is designed to solve the Schrödinger-Poisson system under periodic boundaryconditions. In this case the correct form of equation 2.7 is ∇ Φ = 4 πGm (cid:0) | ψ | − (cid:104)| ψ | (cid:105) (cid:1) , (2.8)where we subtract the average density from the right hand side of the Poisson equation. Theform of Equation 2.8 is a consequence of Gauss’ law and the fact that the surface integral ofthe gradient of the field around the perimeter of the simulation grid is identically zero whenperiodic boundary conditions are imposed [34]. It is helpful to recast the Schrödinger-Poisson system (equations 2.6 and 2.7) in terms ofadimensional quantities. In keeping with Refs [26, 35] we introduce length, time, and massscales as follows: L = (cid:18) π (cid:126) m H Ω m (cid:19) ≈ (cid:18) − eV m (cid:19) kpc , (2.9) T = (cid:18) π H Ω m (cid:19) ≈ . , (2.10) M = 1 G (cid:18) π H Ω m (cid:19) − (cid:18) (cid:126) m (cid:19) ≈ × (cid:18) − eV m (cid:19) M (cid:12) , (2.11)where m is the mass of the ultralight scalar field, H is the present-day Hubble parameter, G is Newton’s gravitational constant and Ω m is the present-day matter fraction of the en-ergy density of the universe. We recast equations 2.6 and 2.7 in terms of the dimensionless For a detailed explication of the WKB approximation in the non-relativistic limit, see [33]. – 4 –uantities t (cid:48) = t T , (cid:126)x (cid:48) = (cid:126)x L , Φ (cid:48) = m T (cid:126) Φ , ψ (cid:48) = T √ mGψ. (2.12)Dropping the primes for notational convenience, we see that the coupled differential equationsof the Schrödinger-Poisson system under periodic boundary conditions reduce to i ˙ ψ ( (cid:126)x, t ) = − ∇ ψ ( (cid:126)x, t ) + Φ( (cid:126)x, t ) ψ ( (cid:126)x, t ) , (2.13) ∇ Φ( (cid:126)x, t ) = 4 π (cid:0) | ψ ( (cid:126)x, t ) | − (cid:104)| ψ ( (cid:126)x, t ) | (cid:105) (cid:1) , (2.14)where it is understood that all quantities involved are dimensionless. We can recover dimen-sionful quantities via the “dictionary” provided by equations 2.9 to 2.11. For example, theintegrated mass of the system, M tot , is given by M tot = M (cid:90) d x | ψ | . (2.15)Likewise, the mass density at any point is given by ρ = ML − | ψ | . (2.16)By dimensional analysis, we can easily restore dimensionful units to any of the quantitiescalculated by the code. In particular, in the following sections it is to be understood that E = ML T − E code , v = LT − v code , (2.17)where E and v represent energy and velocity, respectively. PyUltraLight works internallywith these dimensionless quantities but can receive initial conditions and generates output inphysical units. Henceforth, we will often refer to | ψ | as the density, where it is understoodthat this is in fact a dimensionless quantity related to the physical mass density via theconstant of proportionality given by equation 2.16. In this section we discuss the methodology used to calculate the dynamics of the adimensionalSchrödinger-Poisson system (equations 2.13 and 2.14) given user-specified initial conditions.We introduce the symmetrised split-step Fourier method, and schematically describe how thesystem is evolved at each timestep.
Dynamical evolution within
PyUltraLight progresses via a symmetrised split-step Fourierprocess on an N × N × N grid with periodic spatial boundary conditions. To understand thismethod, first consider the exact expression for the unitary time evolution of the wavefunctionaccording to equation 2.13, namely ψ ( (cid:126)x, t + h ) = T exp (cid:20) − i (cid:90) t + ht dt (cid:48) (cid:26) − ∇ + Φ( (cid:126)x, t (cid:48) ) (cid:27)(cid:21) ψ ( (cid:126)x, t ) , (3.1)where T is the time-ordering symbol. For a sufficiently small timestep h , the trapezoidal rulegives the approximation (cid:90) t + ht dt (cid:48) Φ( (cid:126)x, t (cid:48) ) ≈ h (cid:16) Φ( (cid:126)x, t + h ) + Φ( (cid:126)x, t ) (cid:17) . (3.2)– 5 –e can therefore write the approximate form of equation 3.1 as ψ ( (cid:126)x, t + h ) ≈ exp (cid:20) i h (cid:16) ∇ − Φ( (cid:126)x, t + h ) − Φ( (cid:126)x, t ) (cid:17)(cid:21) ψ ( (cid:126)x, t ) . (3.3)Note that the exponential in equation 3.3 omits the time-ordering symbol, and is only equiv-alent to its time-ordered counterpart to order h .The linear differential operator in equation 3.3 acts naturally in Fourier space, while thenonlinear potential term is simplest to evaluate in position space. By splitting the exponentialwe can evaluate each term in its natural domain. Such a splitting is valid when the timestepis small, and is represented as ψ ( (cid:126)x, t + h ) ≈ exp (cid:20) − ih (cid:126)x, t + h ) (cid:21) exp (cid:20) ih ∇ (cid:21) exp (cid:20) − ih (cid:126)x, t ) (cid:21) ψ ( (cid:126)x, t ) . (3.4)This splitting can be understood thusly: first, a half timestep is taken in which only thenonlinear potential operator acts, followed by a full timestep in the linear term. The potentialfield is then updated, and a final half timestep in the nonlinear term is performed. Using theBaker-Campbell-Hausdorff formula to express the product of exponentials in equation 3.4 asa single exponential and keeping only terms to order h we find: exp (cid:20) i h (cid:16) ∇ − Φ( (cid:126)x, t + h ) − Φ( (cid:126)x, t ) (cid:17) + h (cid:104) ∇ , Φ( (cid:126)x, t ) (cid:105) − h (cid:104) ∇ , Φ( (cid:126)x, t + h ) (cid:105)(cid:21) . (3.5)Making use of the fact that Φ( (cid:126)x, t + h ) ≈ Φ( (cid:126)x, t ) + h ˙Φ( (cid:126)x, t ) we see that the commutators inequation 3.5 cancel at O ( h ) and the expression matches 3.3, with the dominant error termappearing at O ( h ) .Evaluation of equation 3.4 within PyUltraLight thus proceeds as follows: Initially, thenonlinear term acts in position space for one half-timestep. The result is Fourier transformed,and a full timestep is taken with the differential operator applied in the Fourier domain. Thepotential field is then updated in accordance with equation 2.14. After an inverse Fouriertransform a final half timestep is taken with the updated nonlinear term acting in positionspace to give the new ψ field configuration. This procedure is known as the symmetrisedsplit-step Fourier method, and used widely in fields such as nonlinear fiber optics [36–38].The algorithm can be represented schematically as ψ ( (cid:126)x, t + h ) = exp (cid:20) − ih (cid:126)x, t + h ) (cid:21) F − exp (cid:20) − ih k (cid:21) F exp (cid:20) − ih (cid:126)x, t ) (cid:21) ψ ( (cid:126)x, t ) , (3.6)where the order of operations runs from right to left, F and F − denote the discrete Fouriertransform and its inverse, and k is the wavenumber in the Fourier domain. The potentialfield is updated following the inverse Fourier transform in equation 3.6, via Φ( (cid:126)x, t + h ) = F − (cid:18) − k (cid:19) F π | ψ ( (cid:126)x, t i ) | , (3.7)where ψ ( (cid:126)x, t i ) is the field configuration at this halfway point in the full timestep. We explicitlyset the k = 0 Fourier mode to zero prior to the final inverse Fourier transform; as a consequencethere is no need to subtract the global average density from the local value in Equation 3.7,in contrast to Equation 2.8. The final operation in equation 3.6 only changes the phase of– 6 – , so we could replace ψ ( (cid:126)x, t i ) with ψ ( (cid:126)x, t + h ) in equation 3.7 with no change in meaning. PyUltraLight makes an additional simplification to the symmetrised split-step Fouriermethod by combining the consecutive half-steps in the nonlinear term into a single full step.Consequently, only the first and last operations involve actual half steps. Schematically thisbecomes ψ ( t + nh ) = exp (cid:20) + ih (cid:21) (cid:32) n (cid:89) exp [ − ih Φ] exp (cid:20) ih ∇ (cid:21)(cid:33) exp (cid:20) − ih (cid:21) ψ ( t ) , (3.8)where Φ is updated at each step via equation 3.7; attention is drawn to the sign differencebetween the first and last operators.From a computational perspective, the numerical Fourier transforms are likely to bethe rate-limiting step in any pseudospectral code. In PyUltraLight the discrete Fouriertransform (DFT) and its inverse are implemented via pyFFTW, a pythonic wrapper for theC-based FFTW subroutine library which efficiently implements both real and complex DFTs[31, 32, 39]. This allows
PyUltraLight to combine the flexibility of a notebook basedmodelling tool with the efficiency of a carefully tuned, compiled numerical library. FFTWis fully parallelised and its support for multithreading is inherited by pyFFTW and accessedwithin
PyUltraLight ; the number of threads used by the pyfftw.FFTW class is determinedby the Python multiprocessing routines which are used to ascertain the number of availableCPU cores. In addition,
PyUltraLight uses the NumExpr package to parallelise operationson array objects within the simulation [40].
PyUltraLight specifies the initial dark matter configuration as a superposition of an ar-bitrary number of solitonic halos, with arbitrary (user-defined) velocities and phases. Thisis necessarily an idealisation, given that realistic dark matter halos will not map directly tothe solitonic solutions, but it provides an excellent “playground” in which to explore ULDMdynamics, and the initialisation routines within
PyUltraLight can be easily augmented toaccommodate a wider range of scenarios. The initial field configuration is built by loadinga NumPy array file encoding a solitonic solution to the Schrödinger-Poisson system and thecorresponding position mass, velocity, and phase parameters each specified by the user withinthe accompanying Jupyter notebook.In practice, only a finite range of halo masses can be supported within a given simulation– the radius of a solitonic halo is inversely proportional to its mass, so resolving a light halointeracting with a very massive halo would require an extremely fine spatial mesh. However,
PyUltraLight also allows the user to specify a fixed, external potential which does not takepart in the dynamics. At this point only a central /r potential is supported but this wouldbe easily generalised. It should be noted that because PyUltraLight enforces periodicboundary conditions, care must be taken in cases where solitons approach the boundaries ofthe simulation grid. If a soliton were to cross the boundary during a simulation in which aNewtonian central potential is implemented, the forces exerted during the crossing would beunphysical. For studies of orbital stability this is unlikely to cause any problems, as in thesecircumstances material collapses toward the centre of the simulation grid rather than crossingthe boundaries. However, the user should ensure that solitons are initialised sufficiently farfrom the boundary for the purposes of each simulation on a case-by-case basis. In situationswhere a significant portion of the total mass is expected to be ejected, such as the merger– 7 –f multiple solitons to form a larger halo, care should be taken to ensure that mass ejectedabove the escape velocity is not recaptured as it re-enters the grid from the other side. Forstudies of this kind, an absorbing sponge at the grid boundaries is perhaps more suitable thanperiodic boundary conditions, though this has not been implemented in
PyUltraLight atthis stage.The soliton profile used to generate the initial conditions in
PyUltraLight is found byfirst imposing spherical symmetry in the Schrödinger-Poisson equations and assuming timeindependence in the radial density profile [26]: ψ ( (cid:126)x, t ) → e iβt f ( r ) , Φ( (cid:126)x, t ) → ϕ ( r ) , (3.9)where r = | (cid:126)x | . Introducing ˜ ϕ ( r ) := ϕ ( r ) + β , equations 2.13 and 2.14 reduce to − f (cid:48)(cid:48) ( r ) − r f (cid:48) ( r ) + ˜ ϕ ( r ) f ( r ) (3.10) ϕ (cid:48)(cid:48) ( r ) + 2 r ˜ ϕ (cid:48) ( r ) − πf ( r ) (3.11)where primes denote derivatives with respect to r . Note that this system contains no arbitraryconstants, so the underlying profile is effectively universal and is loaded as a pre-computedarray by PyUltraLight , rather than computed from scratch with each code execution.The soliton profile numpy array file is included with
PyUltraLight , however, an auxiliaryprogram soliton _ solution.py is also supplied, from which this array can be generated; ituses a fourth-order Runge-Kutta algorithm to solve the coupled profile equations. We set f ( r ) | r =0 =1, while smoothness requires that first derivatives of f ( r ) and ˜ ϕ ( r ) vanish at theorigin. We then use the shooting method to search for solutions of f ( r ) and ϕ ( r ) satisfyingthe boundary conditions lim r →∞ ϕ ( r ) = 0 and lim r →∞ f ( r ) = 0 , varying ˜ ϕ ( r ) | r =0 until weobtain a solution of f ( r ) which approaches zero at the maximal specified radius, r m . Thevalue of β is then calculated by assuming that ϕ ( r ) goes as − c/r at large radii, where c is aconstant. Under this assumption, we can write ˜ ϕ ( r m ) = − cr m + β, c = r m ˜ ϕ (cid:48) ( r m ) . (3.12)We thus obtain the full solution ψ ( (cid:126)x, t ) = e iβt f ( r ) . Having initially chosen f ( r ) | r =0 = 1 ,we may then generalise to f ( r ) | r =0 = α , where α is an arbitrary positive real number. It iseasily verified that if e iβt f ( r ) is a solution to the spherically symmetric Schrödinger-Poissonsystem, then g ( r ) is also a solution, where g ( r ) = e iαβt αf ( √ α r ) . (3.13)We thus have a family of spherically symmetric soliton solutions to the dimensionless Schrö-dinger-Poisson system; the dimensionless soliton mass is proportional to √ α and the fullwidth at half maximum is proportional to / √ α . Since the size of the soliton scales inverselywith the mass, the most massive soliton in the solution puts a lower bound on the requiredspatial resolution.The Schrödinger equation is not trivially form invariant under Galilean boosts but wecan enforce Galilean covariance through the addition of a velocity-dependent phase factor, ψ ( (cid:126)x, t ) = αf (cid:0) √ α | (cid:126)x − (cid:126)vt | (cid:1) e i ( αβt + (cid:126)v · (cid:126)x − | (cid:126)v | t ) . (3.14)– 8 –o construct the initial field configuration, PyUltraLight , loads the NumPy array encodingthe radial profile f ( r ) for the f ( r ) | r =0 = 1 case. Equation 3.14 is then used to transform thethis solution into soliton(s) with user-specified values position, mass, and velocity specified,via the accompanying Jupyter notebook. The user may also add an additional constant phasefactor if desired. The Courant-Friedrichs-Lewy (CFL) condition is an upper bound on the timestep (as a func-tion of grid-spacing) that must be satisfied by many partial differential equation solversbased on finite-differencing [41] and is often cited in numerical analyses of ULDM via theSchrödinger-Poisson system, see e.g. Ref. [10]. However, the CFL condition expresses acausality constraint, and is generally only strictly applicable to hyperbolic PDEs, whereasthe Schrödinger-Poisson has only a first order time derivative, even though it is effectivelythe nonrelativistic limit of the Klein-Gordon equation. Moreover, because
PyUltraLight computes spatial derivatives via a pseudospectral method, technically it is unconditionallystable [42]. Our split-step algorithm is second order in the timestep, and its value will alwaysbe an empirical tradeoff between computational cost and convergence to the apparent limit inwhich step is arbitrarily small. Consequently, the user is encouraged to validate their choiceof timestep via case-by-case convergence testing.The default timestep in
PyUltraLight is fixed with reference to the fluid interpreta-tion of the Schrödinger-Poisson system [5]. The fluid interpretation is often used to recastthe Schrödinger-Poisson system in the form of the Madelung equations [43], which a hydro-dynamical representation of the system. The first step is to define ψ ≡ √ ρe iθ , (cid:126)v ≡ ∇ θ. (3.15)and to treat (cid:126)v as a fluid velocity. From this perspective, if the phase difference between twoadjacent grid points exceeds π the fluid will appear to move “backwards” across the grid. Wethus set the default timestep, ∆ t , so that fluid travelling at this maximum velocity traversesone grid space, ∆ x , per timestep, or ∆ t = (∆ x ) π . (3.16)This is a choice, rather than a strict constraint on ∆ t . However, if the “fluid” approachesvelocities where the phase appears to switch direction, the configuration is approaching thepoint where the simulation grid is too coarse to fully resolve the dynamics. Hence, a timestepmuch smaller than this value may offer little practical advantage. However, in some cases thebreakdown may occur in regions of the simulation volume that are of little physical interest,and the user is free to choose a larger timestep via the ’step_factor’ parameter in the Jupyternotebook.Alternatively, Ref. [22] fixes the timestep by ensuring that neither of the unitary oper-ators in Equation 3.6 lead to a phase change of more than π for a single grid point over onetimestep. However, because the pseudo-spectral algorithm does not compare the phase of asingle gridpoint at different points in time, this choice of timestep is not a requirement forstability. This method gives the following constraints: ∆ t < (cid:20) π Φ max , x ) π (cid:21) , (3.17)– 9 – − − − D i m e n s i o n l e ss d e n s i t y ( c o d e un i t s ) Soliton Collision Interference Pattern
Numerical profileTheoretical profile
Figure 1 . Comparison of theoretical and numerical density profiles at time of maximal interferencefor head-on collision of two solitons with mass ratio µ = 2 and no relative phase difference. Thesolitons have dimensionless masses 5 and 10, with an initial separation of 4 code units and relativevelocity of 20 code units. The simulation resolution is in a box of side-length 8. where the second of these constraints is generally the stricter of the two, and is equivalent toour default choice of timestep up to a multiplicative factor of O (1) . Our experience is thatspecifying the timestep via Equation 3.16 is suitable for the majority of simulation scenarios,and we explore convergence in more detail in Section 5.2. PyUltraLight
In this section we validate
PyUltraLight by reproducing results from previous studiesof ULDM dynamics, demonstrating interference effects and effective repulsive forces arisingfrom the wavelike nature of ULDM. In addition we study the evolution of the velocity fieldof a solitonic core orbiting within a Newtonian central potential, showing that the stableorbital configuration is an irrotational Riemann-S ellipsoid. Finally, we demonstrate that
PyUltraLight delivers sub-percent level energy conservation for a selection of dynamicalscenarios.
The outcomes of ULDM soliton collisions depend critically on whether the total energy of theisolated binary system is positive or negative. With a positive total energy the solitons passthrough each other, emerging largely undisturbed from their initial configurations and the– 10 –avefunctions describing the solitons are superposed during the collision, yielding distinctiveinterference patterns.Following [10], we consider the head-on collision of two solitions with mass ratio µ = 2 and high relative velocity. While we work in dimensionless code units, it should be noted thata dimensionful velocity can be restored from the code velocity by multiplying through by LT − , the scale parameters defined in Equations 2.9 and 2.10. This simple case of a head-onsoliton collision can be treated approximately. Starting from equation 3.14 we write the totalwavefunction of the binary system in terms of dimensionless quantities along the collision axisas ψ ( x, t ) = α f (cid:0) √ α | x − x − v t | (cid:1) e i ( α βt + v ( x − x ) − v t + δ )+ α f (cid:0) √ α | x − x − v t | (cid:1) e i ( α βt + v ( x − x ) − v t ) , (4.1)where x and x are the initial central positions of the solitions, v and v are the solitonvelocities, δ is a constant relative phase term and α = 4 α , parameterising the density profilesas discussed in Section 3.2. For convenience we set v = − v and x = − x . We expect thatthe interference effects will be maximised when two components of the wavefunction arecentred at the same location, such that x + v t = − x − v t = 0 . This corresponds to atime t o = | x /v | = | x /v | , where in this simplified model we do not account for distortionscaused by the accelerating or compactifying effects that the gravitational interaction has onthe soliton profiles as they approach one another. The dimensionless density is then given by | ψ ( x, t o ) | = α (cid:20) f ( √ α x ) + 16 f (2 √ α x ) +8 f ( √ α x ) f (2 √ α x ) cos (cid:18) − α β (cid:12)(cid:12)(cid:12)(cid:12) x v (cid:12)(cid:12)(cid:12)(cid:12) + 2 v x + δ (cid:19) (cid:21) . (4.2)Figure 1 shows the dimensionless density profile at the time of maximal interference for twosolitons with mass ratio 2 and phase difference δ = 0 . The numerical result obtained using PyUltraLight closely matches the theoretical prediction of equation 4.2. Small disparitiesbetween the numerical and theoretical profiles may be attributed to the effect of gravitationalcontraction not included in the theoretical prediction of equation 4.2 and to a small offset inthe true time of maximal interference due to the solitons accelerating as they fall together.We do not expect an exact match, but we have verified that
PyUltraLight qualitativelyreproduces the wave interference effects of the ULDM model. With the exception of [10], fewstudies of ULDM dynamics have investigated the interference patterns generated by collidingsolitons in this way. In some cases, this is because the algorithm employed to simulate thedynamics is not capable of reproducing such effects. An example of this is given in [17], whereit is demonstrated that the coarse-grained nature of the particle-mesh method renders thealgorithm incapable of reproducing detailed interference patterns such as those shown here.
As demonstrated in [26], the wavelike properties of ULDM give rise to effective forces whichcan dramatically affect the dynamics of core collisions. These effective forces arise as a resultof interference phenomena, rather than because of any local interactions the ULDM modelmight incorporate. Figure 2 shows the results of a head-on collision between two solitons,where in one instance the solitons have no initial phase difference, and in the other instancea phase difference of π is applied in the initial conditions. In this simulation, solitons of– 11 – . . . . . . . . . . . . Phase difference = 0 . . . . . . . . . . . . Phase difference = π T = 0 T = 0.04 T = 0.06 T = 0.08T = 0 T = 0.04 T = 0.06 T = 0.08
Figure 2 . Head-on collisions between solitons of mass 20, initial relative velocity 20, and initialseparation 1.2 in code units. Plots show contours of constant density; time progresses from left toright across each row and is indicated in code units for each frame. The upper panel shows solutonsinitially in phase; the effective repulsive force generated by a π phase shift can be seen in the lowerpanel. mass 20 are initialised with relative velocity 20 and initial separation 1.2 (code units). Thesolitons are allowed to collide, and contours of the density profile along the plane of symmetryare displayed. In one case (top) there is no phase offset between the initial solitons, whilein the second the phases differ by π . In the latter case, the phase shift creates an effectiverepulsive force between the two solitons. It can be seen in the second frame that as the solitonsapproach one another, the π phase shift results in a slowing of the approach accompaniedby a deformation of the density profile, acting so as to avoid contact between the solitons.Dissimilarly, in the case where there is no phase shift, the solitons readily collide and mergeto form a single contiguous density profile prior to re-separating. Further discussion of thisphenomenon and its possible observational consequences can be found in [26]. PyUltraLight allows the inclusion of a static potential equivalent to a point-mass at thecentre of the simulation region. There is no backreaction on this mass as a result of theULDM dynamics, and its “mirror images” within the periodic coordinate systems are notaccounted for within the overall gravitational potential. While a potential of this form doesnot necessarily accurately emulate that which we might expect from a realistic galaxy or darkmatter halo, it provides a starting point for a study of the stability of satellite dark matterhalos orbiting a much larger object. In particular, this includes the investigation of lifespansof dwarf satellite galaxies orbiting much larger objects (including the Milky Way) which area key to understanding whether ULDM models can resolve the so-called missing satellitesproblem [44].An extensive study of the tidal disruption of ULDM solitonic cores orbiting a central– 12 – igure 3 . Contours of constant density for a solitonic core after one revolution around a central po-tential. Contours are superimposed upon the internal velocity field (with the bulk motion subtracted).This velocity field is qualitatively that of an irrotational Riemann-S ellipsoid, and the deformation ofthe density profile of the soliton along the radial direction (red arrow) is visible. The host:satellitemass ratio is approximately 55; simulation resolution is . potential has recently been undertaken in [45] and we reproduce just one of their results here.To do this, we again adopt the definitions 3.15, namely ψ ≡ √ ρe iθ , (cid:126)v ≡ ∇ θ, (4.3)where we are working in dimensionless code units. Using these definitions, the Schrödinger-Poisson system can be recast in terms of hydrodynamical quantities in the so-called Madelungrepresentation. The Madelung equations resemble the continuity and Euler equations of clas-sical fluid dynamics, with the addition of a ‘quantum pressure’ term accounting for resistanceagainst gravitational collapse. The Madelung formalism is discussed in detail in [43, 46–48].Because this hydrodynamical formulation defines the fluid velocity as the gradient of thephase of the field ψ , problems arise when ψ = 0 , where the phase is not well defined. Becauseof this issue, the Madelung and Schrödinger representations are not strictly equivalent unlessa quantisation condition is imposed, as discussed in [49]. We do not consider the subtletiesof the Madelung representation here, as it is sufficient for our purposes to consider the fluidvelocity in the region of a solitonic core, where no field nodes are present. For a discussionof the possible remedies to the ‘nodes problem’, the reader is referred to Chapter 15.3 of It should be noted that, restoring dimensionful units, the fluid velocity (cid:126)v is related to the ususal quantummechanical probability current (cid:126)j through | ψ | (cid:126)v = (cid:126)j = (cid:126) / mi [ ψ ∗ ∇ ψ − ψ ∇ ψ ∗ ] – 13 –50]. Where the Madelung representation is well defined, i.e. where the phase is a smoothlyvarying function, the velocity field of the Schrödinger-Poisson system is strictly irrotational, ∇ × (cid:126)v = 0 . If a radially symmetric soliton is initialised in a circular orbit around a Newtonianpotential, there will be initial transient behaviour as the spherical profile becomes elongatedalong the radial direction of the central potential. Meanwhile, the velocity field correspondingto the overall orbital motion of the soliton will be superposed with the internal velocity field,combining so as to produce a net flow with vanishing curl.The family of Riemann-S ellipsoids describe non-axisymmetric uniformly rotating bodieswhose internal velocity fields have vanishing curl [51]. Therefore, it is the characteristicinternal velocity field of a Riemann-S ellipsoid which we expect to arise during our simulationof a soliton orbiting a central mass. It is found in [45] that an initially spherical solitoniccore without self-rotation will gradually spin up to form a tidally-locked ellipsoid with anirrotational internal velocity field when orbiting a host mass. We reproduce this result using PyUltraLight . Figure 3 shows the internal velocity field of a solitonic satellite after onecomplete revolution around a host mass. The soliton has become elongated along the radialline connecting it to the host, indicating that it is tidally locked, while the velocity field withinthe tidal radius is visibly irrotational and bears the qualitative trademarks of the Riemann-Sellipsoid as presented in Figure 2 of [52].It should be noted that the wider velocity field is not expected to be accurately predictedin a simulation of this kind, though the field within the tidal radius is well-modelled. This isbecause the initial soliton density profile is defined only out to a given cutoff radius, beyondwhich the ψ field value is set identically to zero. As mentioned previously, the Madelunghydrodynamical formulation of the Schrödinger-Poisson system is not valid where ψ = 0 .Because of this, we focus primarily on the internal velocity field within the high density regionof the solitonic core. As we have seen, in this region PyUltraLight is able to accuratelyreproduce the expected velocity field characteristics.
Physically, we expect that the overall energy in the system will be conserved. This provides atest on the numerical performance of
PyUltraLight , and we find that even at relatively lowspatial resolution we see sub-percent level energy conservation for all the dynamical scenariosconsidered here. In this Section we express the energy of the Schrödinger-Poisson system interms of the variables ψ and Φ and discuss its decomposition into individual constituentscalculated separately within the code. We then present results for a variety of configurations.We begin by defining a suitable action which yields the full Schrödinger-Poisson systemthrough its corresponding Euler-Lagrange equations. We find that variation of S = (cid:90) dt (cid:90) R d x − (cid:26) |∇ Φ | + Φ | ψ | + 12 |∇ ψ | + i ψ ˙ ψ ∗ − ˙ ψψ ∗ ) (cid:27) (5.1)with respect to Φ , ψ ∗ and ψ yields equations 2.14, 2.13, and the conjugate of equation 2.13,respectively. The integrand of equation 5.1 is the Lagrangian density, L , from which we canderive the conserved energy in the usual way: E tot = (cid:90) R d x (cid:26) ∂ L ∂ ˙ ψ ˙ ψ + ∂ L ∂ ˙ ψ ∗ ˙ ψ ∗ + ∂ L ∂ ˙Φ ˙Φ − L (cid:27) . (5.2)– 14 – D i s t a n ce T T T T T T T T T T Time − − E GP (self-interaction)Total energy E K + E Q E GP (central potential) Figure 4 . Left: Evolution of the density profile of a solitonic core at equally spaced times as itundergoes tidal disruption in a potential centred at the red cross. Right: Evolution of the componentsof the total energy of the system; times correspond to the labeled density profile snapshots. Allquantities are in dimensionless code units. . . . − − − − − E n e r g y ( c o d e un i t s ) E GP (self-interaction)Total energy E K + E Q Figure 5 . Evolution of the components of the total energy of a binary soliton system with eachsoliton in an elliptical orbit around the common centre of mass. – 15 – . . . − − ∆ E t o t a l / E t o t a l × − Figure 6 . Energy conservation as a function of simulation resolution for a binary solitons orbitingtheir common centre of mass. ∆ E total is the difference between the current total energy and theinitial total energy for the configuration. The ratio of this difference to the initial integrated energyis plotted on the y axis for each resolution. Evaluating this expression, we obtain: E tot = (cid:90) R d x (cid:26) |∇ Φ | + Φ | ψ | + 12 |∇ ψ | (cid:27) (5.3) = (cid:90) R d x (cid:26) ∇ (Φ ∇ Φ) −
12 Φ ∇ Φ + Φ | ψ | + 12 ∇ ( ψ ∗ ∇ ψ ) − ψ ∗ ∇ ψ (cid:27) (5.4) = (cid:90) R d x (cid:26)
12 Φ | ψ | − ψ ∗ ∇ ψ (cid:27) . (5.5)where in the last step we have used Stokes’ Theorem as well as the Poisson equation (2.14)to perform simplifications. Because we are working with the dimensionless quantities definedin equation 2.12, it is easy to see that this quantity is related to the physical energy throughmultiplication by a constant factor of L T − G − . It should be noted that equation 5.5 isnot equivalent to the expectation value of the Schrödinger Hamiltonian, which is itself not aconserved quantity of the Schrödinger-Poisson system and is given by (cid:104) ˆ H (cid:105) = (cid:90) R d x (cid:26) Φ | ψ | − ψ ∗ ∇ ψ (cid:27) . (5.6)The two terms in the integral 5.5 are calculated separately within the code. The first termis the gravitational potential energy of the Schrödinger-Poisson system, E GP . As discussed– 16 – T T T T Time − . − . − . . . . . . ∆ E t o t a l / E t o t a l × − Figure 7 . Energy conservation as a function of simulation resolution for a soliton undergoing signif-icant tidal disruption in a Newtonian central potential. ∆ E total is the difference between the currenttotal energy and the initial total energy for the configuration. The ratio of this difference to the initialintegrated energy is plotted on the y axis for each resolution. in [5], the second term may be decomposed into contributions which may be consideredseparately as kinetic and ‘quantum’ energies, E K and E Q . However, for our purposes itis sufficient to consider only their combined contribution. When PyUltraLight includesthe central potential of a point mass located at the centre of the simulation grid we haveadditional energy contributions and calculate the gravitational potential energy from self-interactions separately from the gravitational potential energy due to the central potential.Figures 4 and 5 demonstrate energy conservation for two scenarios. The first case showsthe evolution of the energy of a single soliton undergoing significant tidal disruption within aNewtonian central potential. For this simulation a soliton of mass 12 in code units was ini-tialised at a radial distance of 3 code units from the centre of a Newtonian central potentialgenerated by a central mass of 1000 code units. As the soliton is disrupted, the kinetic energyincreases, while the gravitational energy due to the central potential decreases, as expected.Meanwhile, the gravitational potential energy from self-interactions gradually increases to-ward zero as the disruption continues and the ULDM is halo spread over a greater area. Inthis case the sum of the individual energy components is conserved to − % at a resolutionof .Figure 5 demonstrates the evolution of the energy of a binary system of solitons in ellip-tical orbits around their common centre of mass over three orbital periods. In dimensionlesscode units, the soliton masses are 22, the initial separation 2, and the initial relative velocitywas 3.6. At points of closest approach the kinetic energy increases as the solitons speed up,– 17 – igure 8 . The configuration used to test the sensitivity of solutions to the spatial and temporalresolution. Two solitons of mass m = 20 are initialised at radial distances r = 2 from a centralmass with M = 1000 moving in opposite directions with initial speeds | v | = (cid:112) M/r , correspondingto clockwise orbits around the central mass. The box size is 10, while the total duration is 0.25 (allquantities in code units). Time runs from left to right. while the potential energy due to self-interaction decreases commensurately such that thetotal energy is conserved. In this scenario no central potential has been included. As thesolitons reach the first point of closest approach, they become slightly deformed, exciting os-cillatory modes which are manifest in the Figure as small scale oscillations superposed on theglobal behaviour. Figure 6 demonstrates the relationship between the total integrated energyand the grid resolution for the same binary system of solitons used to generate Figure 5.The vertical axis shows the ratio of the deviation in the total energy to the initial value ofthe energy, where the deviation is measured as the difference between the current and initialvalues. Energy is conserved at sub-percent level even at low resolutions ( ), and increasinggrid resolution greatly improves accuracy.Figure 7 demonstrates the improvement in energy conservation with increasing gridresolution for a single soliton tidally disrupted in a Newtonian central potential, with thesame set up as used in Figure 4. Namely, a single soliton of mass 12 code units is initialisedat a distance of 3 code units from a central mass, M = 1000 . The initial velocity of thesoliton is (cid:112) M/r where r is the radial distance of the soliton from the central mass. Theduration of the simulation is 0.5 code units so that the soliton undergoes significant tidaldisruption as demonstrated in Figure 4. While we see that energy is conserved at sub-percentlevel even for grid resolution, the qualitative behaviour of the mass density distributionin this case is not correct, so we conclude that this resolution is insufficient for convergencedespite good energy conservation. This highlights the importance of a multifaceted approachto convergence testing. At , energy is conserved to parts in − . We now examine the convergence of the ψ field configuration as a function of spatial resolutionand timestep in a typical simulation. We initialise PyUltraLight with two diametricallyopposed solitons orbiting a large Newtonian central potential, running until the solitons aretidally disrupted, as shown in Figure 8.To examine the sensitivity of the ψ field configuration to the spatial resolution, we firstrun at with the default timestep. We then re-run at resolutions from to withthe timestep fixed to the value and downsample the final outputs to . We sort theresulting values by the density at the corresponding spatial location, and plot differences inthe phase and the magnitude of ψ relative to the values of the run as shown in Figure– 18 – (bottom). The convergence is poor at , but improves with resolution, to the point thatthere is little difference between the and cases.To examine the sensitivity of the the ψ field configuration to the timestep, we take thesame default simulation at , and then compare this to runs with timesteps 0.1, 10, and50 times the default and down-sample the final output arrays to . We sort sort the arrayvalues in order of the ψ field magnitude in the run with the smallest timestep and in Figure10 we show the difference in the phase and magnitude of ψ as a function of the timestep.The difference between the results with the default timestep and a value 10 times smallerare negligible; and there is reasonable agreement between the default case and those with thetimestep boosted by a factor of 10. However, when the timestep is increased by a factor of the accuracy of both the phase and magnitude data are significantly reduced.Figure 11 shows profiles of the density through the simulation volume, as a functionof spatial resolution and timestep. Each plot represents the density profile down the axis ofsymmetry of the initial configuration (vertical axis in Figure 8) after approximately half arevolution around the central potential, or t=0.28 code units – slightly after the final framein Figure 8 - when the solitons have become distorted due to tidal forces, but are not yetcompletely disrupted. We see that as the timestep is varied from 0.1 to 50 times the defaultvalue, the results with the default and the shorter timestep are virtually indistinguishable, andresults are still reasonably accurate at 10 times the default timestep, with small deviationsat high densities. However, the results are significantly distorted at 50 times the defaulttimestep. We also see that as the spatial resolution is decreased from to , the lowestresolution performs poorly, but there is good convergence at resolutions of and above.– 19 – − − P h a s e d e v i a t i o n ( f r o m r e s u l t) D e n s i t y Figure 9 . Top: Deviation of the phase of ψ compared to the highest resolution result ( ). Fieldvalues are arranged in order of increasing magnitude from left to right. A slight improvement in phaseconvergence can be seen for higher density regions to the right. Bottom: Improving convergence of | ψ | with increased spatial resolution for the simulation shown in Figure 8 – 20 – − − P h a s e d e v i a t i o n ( f r o m . x d e f a u l t) D e n s i t y
50 10 1.0 0.1
Figure 10 . Top: Phase deviation of ψ , relative to solution with timestep 0.1 times default, sortedby the density. There is excellent agreement with the default timestep, and reasonable convergenceat steps up to times the default, with better accuracy in high density regions. Bottom: Differencein magnitude of ψ , relative to the solution with timestep 0.1 times default. Again we see goodconvergence with the default timestep, and tolerable agreement in high density regions when the stepis a factor of or less than default. – 21 – osition D e n s i t y Position D e n s i t y Figure 11 . Top: Effect of decreasing the spatial resolution on the density profile at half a revolution.Bottom: Effect of increasing the timestep on the density profile at half a revolution. – 22 –
Discussion and Outlook
PyUltraLight is an accurate, flexible and easy to use tool for studying the dynamicsof ultralight dark matter governed by the Schrödinger-Poisson system of equations. Thecode makes use of a pseudospectral symmetrised split-step Fourier methodology, in whichall spatial derivatives are treated via explicit multiplication in the Fourier domain, therebyavoiding difficulties associated with finite-differencing methods.Energy conservation within
PyUltraLight is excellent, at sub-percent level for sim-ulations run at , with even better performance as resolution is increased. The codecaptures complex phenomena resulting from the wave-like properties of ultralight dark mat-ter, including the interference patterns arising during high-velocity collisions of solitonic cores,and effective forces observed in cases where the colliding cores are out of phase. These phe-nomena can be clearly observed at relatively low spatial resolution, avoiding the need forhigh-performance computing infrastructure to study the fundamental behaviour of ULDMsystems in simple configurations. This makes PyUltraLight a useful tool for investigatingthe dynamics of ULDM systems.
PyUltraLight is Python-based, and as such is particularly simple to understand anduse. The accompanying Jupyter notebook allows for the efficient adjustment of simulationparameters, and offers a useful browser interface for quick visualisation of simulation results.While Python-based, the code makes use of low-level language resources, namely the FFTWlibraries through the use of the Pythonic pyFFTW wrapper and will operate at ∼ efficiency on a 16 core desktop workstation, suggesting that it is computationally efficient.The current implementation of PyUltraLight is already a useful tool for simulatingdynamical ULDM systems and exploring their dynamics However, there is much scope forimprovement. In particular, future releases may incorporate a variable timestep and more so-phisticated physics, including explicit self-interactions in the axion sector or additional mattercomponents. Augmented versions of the code may also include higher-order generalisations ofthe pseudo-spectral method, such as those used in [53].
PyUltraLight is publicly availableunder a BSD license.
Acknowledgments
We thank Xiaolong Du, Lam Hui, David Marsh, Nathan Musoke, Jens Niemeyer, and ChandaPrescod-Weinstein for valuable discussions on axion / ULDM cosmology, and thank Miro Erk-intalo for advice on Schrödinger-Poisson solvers in optical systems. We acknowledge supportfrom the Marsden Fund of the Royal Society of New Zealand, and the use of the New ZealandeScience Infrastructure (NeSI) high-performance computing facilities, which are funded jointlyby NeSI’s collaborator institutions and through the Ministry of Business, Innovation & Em-ployment’s Research Infrastructure programme . A Download and Licensing
PyUltraLight is publicly available under a BSD licence. The full repository, includingsupplementary files such as the code used to generate soliton profiles, is available on GitHubat https://github.com/auckland-cosmo/PyUltraLight . PyUltraLight makes use of thepyFFTW pythonic wrapper around the FFTW C-based fast Fourier transform libraries. BothpyFFTW and FFTW are freely-available and
PyUltraLight has been used successfully on– 23 –oth Mac OS and Linux systems, as well as a shared-memory cluster environment. Wewelcome advice and feedback from users.
References [1]
Planck
Collaboration, P. A. R. Ade et al.,
Planck 2015 results. XIII. Cosmological parameters , Astron. Astrophys. (2016) A13, [ arXiv:1502.01589 ].[2]
PandaX-II
Collaboration, A. Tan et al.,
Dark Matter Results from First 98.7 Days of Datafrom the PandaX-II Experiment , Phys. Rev. Lett. (2016), no. 12 121303,[ arXiv:1607.07400 ].[3]
LUX
Collaboration, D. S. Akerib et al.,
Results from a search for dark matter in the completeLUX exposure , Phys. Rev. Lett. (2017), no. 2 021303, [ arXiv:1608.07648 ].[4] P. Bull et al.,
Beyond Λ CDM: Problems, solutions, and the road ahead , Phys. Dark Univ. (2016) 56–99, [ arXiv:1512.05356 ].[5] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Ultralight scalars as cosmological darkmatter , Phys. Rev.
D95 (2017), no. 4 043541, [ arXiv:1610.08297 ].[6] J. E. Kim and G. Carosi,
Axions and the Strong CP Problem , Rev. Mod. Phys. (2010)557–602, [ arXiv:0807.3125 ].[7] D. J. E. Marsh, Axion Cosmology , Phys. Rept. (2016) 1–79, [ arXiv:1510.07633 ].[8] W. Hu, R. Barkana, and A. Gruzinov,
Cold and fuzzy dark matter , Phys. Rev. Lett. (2000)1158–1161, [ astro-ph/0003365 ].[9] D. J. E. Marsh and A.-R. Pop, Axion dark matter, solitons and the cusp-core problem , Mon.Not. Roy. Astron. Soc. (2015), no. 3 2479–2492, [ arXiv:1502.03456 ].[10] B. Schwabe, J. C. Niemeyer, and J. F. Engels,
Simulations of solitonic core mergers in ultralightaxion dark matter cosmologies , Phys. Rev.
D94 (2016), no. 4 043513, [ arXiv:1606.05151 ].[11] E. Madelung,
Eine anschauliche Deutung der Gleichung von Schrödinger , Naturwissenschaften (Nov, 1926) 1004–1004.[12] J. Zhang, Y.-L. S. Tsai, J.-L. Kuo, K. Cheung, and M.-C. Chu, Ultralight Axion Dark Matterand Its Impact on Dark Halo Structure in N -body Simulations , Astrophys. J. (2018), no. 151, [ arXiv:1611.00892 ].[13] V. Springel,
The Cosmological simulation code GADGET-2 , Mon. Not. Roy. Astron. Soc. (2005) 1105–1134, [ astro-ph/0505010 ].[14] “https://github.com/liambx/Axion-Gadget.”[15] M. Nori and M. Baldi,
AX-GADGET: a new code for cosmological simulations of Fuzzy DarkMatter and Axion models , Mon. Not. Roy. Astron. Soc. (2018) 3935, [ arXiv:1801.08144 ].[16] P. Mocz and S. Succi,
Numerical solution of the nonlinear Schrödinger equation usingsmoothed-particle hydrodynamics , Phys. Rev.
E91 (2015), no. 5 053304, [ arXiv:1503.03869 ].[17] J. Veltmaat and J. C. Niemeyer,
Cosmological particle-in-cell simulations with ultralight axiondark matter , Phys. Rev.
D94 (2016), no. 12 123523, [ arXiv:1608.00802 ].[18] A. Almgren, J. Bell, M. Lijewski, Z. Lukic, and E. Van Andel,
Nyx: A Massively Parallel AMRCode for Computational Cosmology , Astrophys. J. (2013) 39, [ arXiv:1301.4498 ].[19] A. J. Benson,
Galacticus: A Semi-Analytic Model of Galaxy Formation , New Astron. (2012)175–197, [ arXiv:1008.1786 ].[20] X. Du, C. Behrens, and J. C. Niemeyer, Substructure of fuzzy dark matter haloes , Mon. Not.Roy. Astron. Soc. (2017), no. 1 941–951, [ arXiv:1608.02575 ]. – 24 –
21] V. Springel,
E pur si muove: Galiliean-invariant cosmological hydrodynamical simulations on amoving mesh , Mon. Not. Roy. Astron. Soc. (2010) 791, [ arXiv:0901.4107 ].[22] P. Mocz, M. Vogelsberger, V. H. Robles, et al.,
Galaxy formation with BECDM - I. Turbulenceand relaxation of idealized haloes , Mon. Not. Roy. Astron. Soc. (2017), no. 4 4559–4570,[ arXiv:1705.05845 ].[23] H.-Y. Schive, Y.-C. Tsai, and T. Chiueh,
GAMER: a GPU-Accelerated Adaptive MeshRefinement Code for Astrophysics , Astrophys. J. Suppl. (2010) 457–484,[ arXiv:0907.3390 ].[24] “https://github.com/gamer-project/gamer/wiki.”[25] H.-Y. Schive, T. Chiueh, and T. Broadhurst,
Cosmic Structure as the Quantum Interference ofa Coherent Dark Wave , Nature Phys. (2014) 496–499, [ arXiv:1406.6586 ].[26] A. Paredes and H. Michinel, Interference of Dark Matter Solitons and Galactic Offsets , Phys.Dark Univ. (2016) 50–55, [ arXiv:1512.05121 ].[27] R. Easther, H. Finkel, and N. Roth, PSpectRe: A Pseudo-Spectral Code for (P)reheating , JCAP (2010) 025, [ arXiv:1005.1921 ].[28] M. A. Amin, R. Easther, and H. Finkel,
Inflaton Fragmentation and Oscillon Formation inThree Dimensions , JCAP (2010) 001, [ arXiv:1009.2505 ].[29] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg,
Oscillons After Inflation , Phys. Rev. Lett. (2012) 241302, [ arXiv:1106.3335 ].[30] S.-Y. Zhou, E. J. Copeland, R. Easther, et al.,
Gravitational Waves from Oscillon Preheating , JHEP (2013) 026, [ arXiv:1304.6094 Signature Changing Spacetimes and WKB Approximations in GeneralRelativity . Honors thesis, The University of Texas at Austin, 2015.[34] I. Dabo, B. Kozinsky, N. E. Singh-Miller, and N. Marzari,
Electrostatics in periodic boundaryconditions and real-space corrections , Phys. Rev. B (Mar, 2008) 115139.[35] H.-Y. Schive, M.-H. Liao, T.-P. Woo, et al., Understanding the Core-Halo Relation of QuantumWave Dark Matter from 3D Simulations , Phys. Rev. Lett. (2014), no. 26 261302,[ arXiv:1407.7762 ].[36] G. Agrawal,
Nonlinear Fiber Optics . Optics and Photonics Series. Academic Press, 2013.[37] Q. Zhang and M. I. Hayee,
Symmetrized Split-Step Fourier Scheme to Control GlobalSimulation Accuracy in Fiber-Optic Communication Systems , Journal of Lightwave Technology (Jan, 2008) 302–316.[38] O. V. Sinkin, R. Holzlohner, J. Zweck, and C. R. Menyuk, Optimization of the split-stepFourier method in modeling optical-fiber communications systems , Journal of LightwaveTechnology (Jan, 2003) 61–68.[39] M. Frigo and S. G. Johnson, The Design and Implementation of FFTW3 , Proceedings of theIEEE (2005), no. 2 216–231. Special issue on “Program Generation, Optimization, andPlatform Adaptation”.[40] “http://numexpr.readthedocs.io/.”[41] M. A. Ajaib, Numerical Methods and Causality in Physics , arXiv:1302.5601 .[42] T. R. Taha and M. J. Ablowitz, Analytical and Numerical Aspects of Certain NonlinearEvolution Equations. II. Numerical, Nonlinear Schrodinger Equation , J. Comput. Phys. (1984) 203–230. – 25 –
43] A. Suárez and T. Matos,
Structure Formation with Scalar Field Dark Matter: The FluidApproach , Mon. Not. Roy. Astron. Soc. (2011) 87, [ arXiv:1101.4039 ].[44] D. H. Weinberg, J. S. Bullock, F. Governato, R. Kuzio de Naray, and A. H. G. Peter,
Cold darkmatter: controversies on small scales , Proc. Nat. Acad. Sci. (2015) 12249–12255,[ arXiv:1306.0913 ].[45] X. Du, B. Schwabe, J. C. Niemeyer, and D. Bürger,
Tidal disruption of fuzzy dark mattersubhalo cores , Phys. Rev.
D97 (2018), no. 6 063507, [ arXiv:1801.04864 ].[46] A. Suárez and P.-H. Chavanis,
Hydrodynamic representation of the Klein-Gordon-Einsteinequations in the weak field limit , J. Phys. Conf. Ser. (2015), no. 1 012008,[ arXiv:1504.01164 ].[47] R. Johnston, A. N. Lasenby, and M. P. Hobson,
Cosmological fluid dynamics in theSchrödinger formalism , arXiv:0904.0611 .[48] M. Kopp, K. Vattis, and C. Skordis, Solving the Vlasov equation in two spatial dimensions withthe Schrödinger method , Phys. Rev.
D96 (2017), no. 12 123532, [ arXiv:1711.00140 ].[49] T. C. Wallstrom,
Inequivalence between the Schrodinger equation and the Madelunghydrodynamic equations , Phys. Rev.
A49 (1994) 1613–1617.[50] R. E. Wyatt,
Quantum Dynamics with Trajectories. Introduction to quantum hydrodynamics ,vol. 28. 2005.[51] S. Chandrasekhar,
The Equilibrum and the Stability of the Riemann Ellipsoids. I. , Astrophys.J. (Oct., 1965) 890.[52] T. Rindler-Daller and P. R. Shapiro,
Angular Momentum and Vortex Formation inBose-Einstein-Condensed Cold Dark Matter Haloes , Mon. Not. Roy. Astron. Soc. (2012)135–161, [ arXiv:1106.1256 ].[53] D. G. Levkov, A. G. Panin, and I. I. Tkachev,
Gravitational Bose-Einstein condensation in thekinetic regime , arXiv:1804.05857 ..