AxioNyx: Simulating Mixed Fuzzy and Cold Dark Matter
Bodo Schwabe, Mateja Gosenca, Christoph Behrens, Jens C. Niemeyer, Richard Easther
AAxioNyx : Simulating Mixed Fuzzy and Cold Dark Matter
Bodo Schwabe, ∗ Mateja Gosenca, † Christoph Behrens, ‡ Jens C. Niemeyer,
1, 2, § and Richard Easther ¶ Institut f¨ur Astrophysik, Universit¨at G¨ottingen, Germany Department of Physics, University of Auckland, New Zealand (Dated: July 17, 2020)The distinctive effects of fuzzy dark matter are most visible at non-linear galactic scales. Wepresent the first simulations of mixed fuzzy and cold dark matter, obtained with an extended versionof the
Nyx code. Fuzzy (or ultralight, or axion-like) dark matter dynamics are governed by thecomoving Schr¨odinger-Poisson equation. This is evolved with a pseudospectral algorithm on theroot grid, and with finite differencing at up to six levels of adaptive refinement. Cold dark matteris evolved with the existing N-body implementation in
Nyx . We present the first investigationsof spherical collapse in mixed dark matter models, focusing on radial density profiles, velocityspectra and soliton formation in collapsed halos. We find that the effective granule masses decreasein proportion to the fraction of fuzzy dark matter which quadratically suppresses soliton growth,and that a central soliton only forms if the fuzzy dark matter fraction is greater than 10%. The
Nyx framework supports baryonic physics and key astrophysical processes such as star formation.Consequently,
AxioNyx will enable increasingly realistic studies of fuzzy dark matter astrophysics.
I. INTRODUCTION
The physical nature of dark matter is a major openquestion for both astrophysics and particle physics.Axion-like particles (ALPs) are viable candidates [1–12];see Refs [13, 14] for recent reviews. Among them, ultra-light axions predict the strongest differences from colddark matter (CDM) on galactic scales. String theory mo-tivates the existence of large numbers of scalar fields aris-ing from the compactification of extra dimensions withinherently low masses and weak interactions [15, 16]. Sce-narios with more than a single species of dark matter, ormixed dark matter (MDM) models, arise naturally fromthis perspective.Independent of their fundamental origin, extremelylight scalar particles with negligible interactions pro-duced in a coherent non-thermal state via the misalign-ment mechanism are typically referred to as Fuzzy DarkMatter (FDM) [4] or Ultralight Dark Matter (ULDM) [5].Unlike QCD axions, their masses can be low enough toexhibit wave-like behaviour on astrophysical scales, upto distances on the order of a kiloparsec given currentbounds on the axion mass. Chaotic interference struc-tures with strong, short-lived density variations and ha-los with gravitationally-bound central solitons are amongthe potentially observable wave-like phenomena [17].These structures have no direct analogues in CDM. How-ever, on scales greater than the de Broglie wavelength, λ dB = h/mv , the properties of collapsed structures areindistinguishable from those of pressureless CDM, via theSchr¨odinger-Vlasov correspondence [18, 19]. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] Small-scale differences open the possibility of observa-tional comparisons of FDM and CDM. Moreover, the ob-served small scale properties of galaxies may be in tensionwith simple CDM scenarios, and FDM appears to amelio-rate some of these issues [17, 20–24]. On the other hand,FDM solitons alone may not explain the observed radialdensity profiles of galactic cores [25–28]. However, thesescales are also heavily influenced by baryonic physics, starformation, supernovae, environmental effects and galaxy-galaxy interactions, all of which operate in the collapsed,nonlinear regime. Consequently, obtaining accurate pre-dictions for FDM dynamics in realistic astrophysical en-vironments is a challenging task even for state-of-the-artsimulations on modern supercomputers.FDM is composed of non-relativistic scalar matter de-scribed by a wavefunction ψ , which both sources and in-teracts with the Newtonian gravitational potential. Con-sequently, FDM is governed by the coupled Schr¨odinger-Poisson equation. The local FDM velocity is representedby the gradient of the phase of the complex-valued wave-function. Extracting this from simulations requires thattheir spatial resolution is fine enough to resolve λ dB , notonly in small high-density regions, but also in exten-sive low density regions if high-speed flows are present.This makes full cosmological simulations with FDM sig-nificantly more challenging than their CDM analogues,which can be efficiently evolved in phase space with N-body algorithms.Several different approaches for solving theSchr¨odinger-Poisson equations numerically have beenemployed in the context of cosmology, see [14] for anoverview. In very large volumes that fail to resolve λ dB , N-body simulations with FDM initial conditionsadequately reproduce the suppression of small-scaleclustering [29–33]. Modified hydro solvers including a“quantum pressure” term motivated by the Madelungtransformation of the Schr¨odinger equation have beenused in the weakly nonlinear regime [33–37].However,only methods that solve the Schr¨odinger equation a r X i v : . [ a s t r o - ph . C O ] J u l directly, capture the fully nonlinear wave-like dynamicsin collapsed FDM structures, either in the entire com-putational domain [17, 38–40] or in small subvolumes ofa hybrid N-body-Schr¨odinger scheme [41]. Recently, thefirst hydrodynamical simulations with FDM includingbaryonic feedback on galactic scales were presented inRefs. [42, 43].This paper introduces AxioNyx , a Schr¨odinger-Poisson solver on adaptively refined regular meshes tofacilitate detailed explorations of the astrophysical con-sequences of FDM. Built within the existing cosmologycode Nyx [44], it inherits support for CDM and baryonicmatter by means of a particle-mesh N-body scheme and ahigher-order unsplit Godunov method for the gas dynam-ical equations, respectively.
Nyx itself is built upon AM-ReX, a powerful framework for block-structured adaptivemesh refinement (AMR) applications on massively par-allel supercomputers [45], boosting resolution in regionsof interest and allowing the modeling of structure for-mation over many length scales.
AxioNyx implementsboth pseudospectral and finite-difference methods for theSchr¨odinger-Poisson equation on the root grid, and usesa finite-difference method in refined regions.We examine the collapse of a spherical overdensity ina universe containing a mixture of cold and fuzzy darkmatter as a test problem for
AxioNyx . This scenariois motivated by the many naturally light and weakly-interacting scalar particles predicted by string theory[15]. In our setup, they are represented by two dark mat-ter components, one assumed to be sufficiently massiveto be well-described by the standard N-body method andthe other exhibiting wave-like behavior governed by theSchr¨odinger-Poisson equation.
AxioNyx is a cosmological solver but the Schr¨odinger-Poisson equation arises in several contexts, including bo-son stars [39, 40, 46] and QCD axion miniclusters [47].
AxioNyx can easily be adapted to address these cases,and boson star condensation [48] is one of the testsused to verify the code (see Appendix A 3). Further, inmany early-universe scenarios the post-inflationary uni-verse can host the formation of transient, gravitationallycollapsed overdensities, whose evolution is governed bythe Schr¨odinger-Poisson equation [49, 50].The paper is organized as follows. In Section II weintroduce our numerical methods, which are validatedagainst a suite of test cases in the Appendices. We theninvestigate spherical collapse of MDM in Section III andconclude in Section IV. https://github.com/axionyx α c α d α c α d α II. NUMERICAL METHODS
The dynamics of FDM is well described by theSchr¨odinger-Poisson equation i (cid:126) ∂ψ∂t = − (cid:126) ma ∇ ψ + mV ψ , (1) ∇ V = 4 πGa ρ , ρ = | ψ | . (2)Here, ψ is the FDM wave function, V denotes the grav-itational potential, and a is the scale factor. The con-stants (cid:126) , m , and G stand for the Planck’s constant, themass of the axion, and Newton’s gravitational constant,respectively.We employ the Nyx code [44], modified to handleFDM physics. Nyx is an adaptive mesh refinement codewith an MPI/OpenMP parallelization scheme, based onthe AMReX [45] library, that is known to scale well, evenover thousands of distributed CPUs. It encompasses rou-tines for hydrodynamics and an N-body solver for CDMphysics.If periodic boundary conditions are applicable, which isoften the case in cosmological scenarios, a pseudospectralsolver is the fastest and most accurate algorithm to evolvethe system of equations 1 and 2 on the Eulerian rootgrid. However, the refined regions will be aperiodic andwe thus use a finite-difference scheme on subgrids.The finite differencing scheme uses a fourth orderRunge-Kutta solver to evolve equations 1 and 2, and hadpreviously been integrated with
Nyx [39]. As an explicitalgorithm, it suffers from stringent time-step constraintsthat scale quadratically with resolution, motivating theimplementation of a pseudospectral solver.The public AMReX repository already includesthe stand-alone version of HACC’s distributed-memory,pencil-decomposed, parallel 3D FFT routines and wrap-pers linking those routines to the AMReX framework. We optimized the HACC libraries for OpenMP paral-lelization and obtain good weak scaling at up to 512 https://github.com/AMReX-Astro/Nyx https://github.com/AMReX-Codes/amrex https://xgitlab.cels.anl.gov/hacc/SWFFT
512 kpc 16 kpc4 kpc
FIG. 1. Slice through the final FDM density field with 30% FDM and 70% CDM. The central region is highly refined, andresolves the FDM granules which form after shell crossing. Rectangular boxes show regions with the same level of refinement– the box edges are coloured such that they are darker at higher levels of refinement. cores, quantified by the negligible dependence of the so-lution time with respect to the number of processors fora fixed problem size per processor. We implement twodifferent pseudospectral schemes, one at second orderand the other at sixth order. The second order schemematches that implemented in
PyUltraLight [51]. Thesixth order scheme replicates the method described inRef [48] ψ ( t + ∆ t ) = α (cid:89) e − imd α ∆ tV α ( x ) / (cid:126) e − i (cid:126) c α ∆ tk / ma ψ ( t ) , (3)with weights summarized in Table I. Other higher-orderalgorithms could easily be implemented. https://github.com/auckland-cosmo/PyUltraLight Eulerian grids used to deposit information from N-body particles in cloud-in-cell simulations appropriate forCDM are typically only adaptively refined in over-denseregions to correctly capture the strong dynamics in fila-ments and halos. Conversely, in FDM simulations the ve-locity is inferred from the gradient of the wave function’scomplex phase; even low density regions may have to berefined if they contain high-speed flows. Consequently,we employ the L¨ohner error estimator [52], which is alsoimplemented in Enzo [53] and
FLASH [54] and was alsoused with FDM cosmological simulations in Refs [17, 55].However, as in Ref. [17], we prohibit grid refinement inregions below a certain density, chosen so that their de-tailed behaviour does not significantly alter the overalldynamics.The power of adaptive refinement is demonstrated inFigure 1, which shows the final snapshot of a sphericalcollapse with 30% FDM and 70% CDM. Six levels of re-finement were used, increasing the resolution by a factor scale factor10 m a x . C D M o v e r d e n s i t y F D M f r a c t i o n FIG. 2. Evolution of the maximum CDM over-density duringspherical collapse for different CDM fractions (1 − f ). Thenumerically obtained power law growth (solid lines) is wellapproximated by Equation 6 (dashed lines). of 2 with each level.We conducted a number of tests of the validity andaccuracy of AxioNyx , as described in the Appendix.These begin with comparisons of the pseudospectralSchr¨odinger solver with static potentials for which exactsolutions are known, verifying the scaling with timestepsize. In A 2 we verify that isolated solitons in a staticbackground behave as expected, showing small depar-tures from perfectly static solutions on small grids. InA 3 we replicate simulations of boson star condensationfrom random Gaussian initial conditions [48]. Since allregions are equally important, this problem is particu-larly suitable for a pseudospectral solver which allowsthe use of significantly larger time steps and lower spa-tial resolution than finite difference methods.The cosmological dynamics are verified by followingthe linear evolution of a one-dimensional, self-gravitating,sinusoidal over-density as it crosses the FDM Jeans scalein a comoving box (A 4). Finally, in A 5 we investigatesupport against spherical gravitational collapse near theJeans scale. This setup is useful to test the L¨ohner re-finement criterion and subcycling; even extreme caseswith more than 100 sub-steps between root and first levelyields well converged results. Results from a representa- tive simulations are shown in Figure 1.
III. SPHERICAL COLLAPSE OF MIXED DARKMATTER
In this section we extend the spherical collapse to aMDM scenario comprised of FDM and CDM with thefraction of FDM defined as f = ρ FDM ρ FDM + ρ CDM . (4)For FDM, we use the sixth order Schr¨odinger-Poissonsolver on the root grid and the finite difference algorithm scale factor10 m a x . F D M o v e r d e n s i t y F D M f r a c t i o n FIG. 3. Evolution of the maximum FDM over-density duringspherical collapse for different f . The initially homogeneousFDM fields (solid line) contract according to linear theory(dotted lines). The power law growth at low redshifts is inde-pendent of FDM initial conditions as small initial FDM over-densities (dashed line) show the same late time behaviour. on higher levels, for the CDM we make use of the N-bodyscheme as implemented in Nyx . A. The linear regime
The linear evolution of MDM is governed by the fol-lowing system of coupled differential equations:¨ δ FDM + 2 H ˙ δ FDM + (cid:18) k (cid:126) m a − πGf ρ (cid:19) δ FDM = 4 πG (1 − f ) ρδ CDM , (5a)¨ δ CDM + 2 H ˙ δ CDM − πG (1 − f ) ρδ CDM = 4 πGf ρδ
FDM . (5b)We choose a box size an order of magnitude belowthe FDM Jeans scale. Accordingly, the growth of anylocal deviation from the critical density is highly sup- pressed in the FDM sector. It is therefore reasonableto start with δ FDM ( z ini ) = 0. In contrast, we intro-duce a Gaussian overdensity with maximum amplitude f = 0.01 10 f = 0.01 FDMCDMf = 0.10 10 f = 0.10f = 0.30 10 f = 0.30f = 0.50 10 f = 0.50f = 0.75 10 f = 0.750 20 40 60 80f = 1.00 10 f = 1.00velocity [km/s] n o r m a li z e d d i s t r i b u t i o n radius [kpc] d e n s i t y [ M / p c ] FIG. 4. (Left) Normalized velocity spectra of the FDM and CDM components for various f . For large f they are well fitted byMaxwell-Boltzmann distributions shown in black. Dashed horizontal lines indicate soliton velocities as defined in Equation 8.(Right) Radial density profiles of the FDM and CDM components for various f . For large f central FDM profiles are well fitby solitonic profiles shown in black. Dashed vertical lines indicate soliton radii while solid vertical lines represent virial radii. δ CDM ( z ini ) = 0 .
001 on top of a CDM background at aninitial redshift z ini = 99. The time evolution of δ CDM isshown in Figure 2. As expected, it is suppressed relativeto a pure CDM simulation [56], via δ CDM ( a ) ∝ a ( √ − f ) − / . (6)Equation 6 can be obtained by setting δ FDM = 0 inEquation 5b. In Figure 3 we show the numerically ob-tained maximum FDM overdensity. It is well approx-imated by the solution of Equation 5a for δ CDM ( a ) asgiven in Equation 6. The late time evolution of modes k (cid:29) k J ( a ) follows a single power law δ FDM ( a ) ∝ a ( √ − f )+3) / . (7)We verify that the growth is not an artifact of the ho-mogeneous FDM initial conditions by running the samesimulation with an added FDM overdensity with initial maximum amplitude δ FDM ( z ini ) = δ CDM ( z ini ) = 0 . B. The non-linear regime
Employing adaptive mesh refinement with a 1024 rootgrid and up to 6 levels of refinement, we investigate spher-ical collapse all the way into the highly non-linear regimeafter shell crossing. We ran simulations with 11 differentmixed dark matter fractions f in a 2 Mpc comoving box.The typical mass of the collapsed halo is 3 . × solarmasses. The AMR structure of one of our simulationscan be seen in Figure 1, which depicts slices through the m a x . F D M o v e r d e n s i t y F D M f r a c t i o n FIG. 5. Maximum central FDM densities normalized by ini-tial values. Soliton formation breaks the degenerate evolutionaround a scale factor a = 0 . final FDM density. Focusing on the central halo regionwe see the expected granular structure in FDM as a resultof interfering plane waves.In Figure 4 we show velocity distributions on the leftside and density profiles on the right. Due to FDM ex-hibiting oscillations in time, we show six different snap-shots, extracted from the final stages of our simulations.They are shown in transparent blue or red lines for CDMand FDM, respectively. The thick red (blue) lines showprofiles averaged over time.In the outer NFW-like halo f remains constant evenafter non-linear collapse and velocity spectra are onlymildly sensitive to different f . This demonstrates thatthe Schr¨odinger-Vlasov correspondence also holds formixed dark matter, i.e. that FDM averaged over multiplede Broglie wavelengths behaves as CDM. This excludesformation of solitons which have no CDM analog. Theintriguing conclusion is that analytic estimates of con-densation [48], heating and cooling [57], and relaxationprocesses [16] in the FDM sector which depend on the ef-fective granule mass can be straightforwardly generalizedto mixed dark matter scenarios by re-scaling the effectivegranule mass by f .Within the central region of a pure FDM halo, sim-ulations reveal the formation of solitonic ground statesolutions. These solutions are present only if f (cid:38) .
1. Inthis case the inner part of the FDM radial density profilescan be well approximated by a soliton profile as depictedin Figure 4. For lower f we see strong fluctuations in thecentral FDM density profiles, which cannot be fitted bythe original or the modified soliton profiles [43]. We thusconclude that for mixed dark matter soliton formationonly occurs if f (cid:38) .
1. In that case the central density ρ c significantly increases beyond the initial collapse, break-ing its initially degenerate evolution, as seen in Figure 5.The vertical dashed lines on the right side of Figure 4indicate the soliton radius r c at which the spherically F D M f r a c t i o n radius [kpc] v e l o c i t y [ k m / s ] FIG. 6. Peaks of the final FDM (diamonds) and CDM(crosses) velocity spectra in the inner halo region as presentedin Figure 4 as a function of soliton radii for different f . FDMvelocitiy peaks tightly follow Equation 8 indicated by thesolid black line. For comparison, the dashed line indicatesthe halo’s virial velocity. averaged FDM density drops to half its central value.The corresponding soliton velocities [40] v c = 2 π . (cid:126) mr c (8)are represented by the dashed lines on the left. Theyalign well with the peaks of the Maxwell-Boltzmann-likedistributed, normalized FDM velocity spectra [41] f ( v ) = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) d x exp [ − im v · x / (cid:126) ] ψ (x) (cid:12)(cid:12)(cid:12)(cid:12) . (9)We integrated the final FDM state in the central 8 kpccubed box, which encompasses at least dozens of gran-ules apart from the soliton. The tight correlation be-tween maxima in FDM velocity spectra and soliton ve-locities shown in Figure 6 suggests that the soliton is inkinetic equilibrium with its surrounding. During solitonformation central velocities both in the FDM and CDMcomponent decouple from the virial velocity of the haloindicated by the dashed, horizontal line in Figure 6. Incontrast, in simulations with insufficient FDM content forsoliton formation ( f (cid:46) . A ( t ) = A · ( t − t ) /τ gr + A f / (10)with τ gr = 0 . √ π m v c G ρ c Λ (cid:39) . t c Λ (11)grows linearly in time. Here, Λ = log( r vir /r c ) is theCoulomb logarithm, r vir the virial radius, and in the last m a x . F D M a m p li t u d e [( M / M p c ) / ] F D M f r a c t i o n FIG. 7. Linear growth of maximum FDM amplitudes welldescribed by Equation 10. step we used soliton identities as listed in the Appendix Bof [16]. Numerically, we found A = 1350 ( M (cid:12) / Mpc ) / , A = 10 ( M (cid:12) / Mpc ) / , and t = 1 . t ff with free-falltime t ff = 0 .
003 s Mpc/km. Since v c only mildly de-pends on the FDM fraction while density scales linearlywith f , growth in soliton amplitude is suppressed roughlyquadratically with f . IV. CONCLUSIONS
We developed the highly-parallelised numerical code
AxioNyx for simulations of self-gravitating mixed fuzzyand cold dark matter. The code features adaptive meshrefinement which allows tremendous effective resolutionin regions of interest, while remaining relatively cheap onthe largest scales. We used a 6th order pseudospectralsolver on the root grid and a 4th order finite-differencesolver on the refined subgrids. We employ the Poissonsolver available in
Nyx . An exhaustive test suite of pro-gressively more complex problems demonstrates the cor-rectness and efficiency of
AxioNyx and highlights itsvast range of potential applications.With mixed dark matter spherical collapse simulationswe show that below the Jeans scale fuzzy dark mat-ter overdensities δ FDM ( a ) are supported against gravita-tional collapse by gradient energy. In turn, cold darkmatter collapse is increasingly suppressed with higherfractions of fuzzy dark matter. The deepening gravita-tional potential shrinks the Jeans scale and results in latetime fuzzy dark matter collapse. We provide analytic es-timates for the evolution of δ FDM ( a ) and δ CDM ( a ) in thelinear regime.Adaptively refined spherical collapse simulations well in the non-linear regime after shell-crossing confirm thatthe Schr¨odinger-Vlasov correspondence holds even in thecase of mixed dark matter. Averaging over multiple deBroglie wavelengths we see a confluent evolution of coldand fuzzy dark matter due to both species respondingto the same gravitational potential. Analytic estimatesof fuzzy dark matter soliton condensation, heating andcooling, and relaxation processes all depend on the effec-tive granule mass. They can thus be straightforwardlygeneralized to mixed dark matter scenarios by re-scalingthe effective granule mass by the fuzzy dark matter frac-tion.The degeneracy between FDM and CDM is broken bythe soliton formation which occurs in the center of col-lapsing overdensities if the fraction of fuzzy dark matteris sufficiently large, f (cid:38) .
1. While staying in kineticequilibrium, the linearly growing soliton heats up its sur-rounding inner halo region.Our adaptively refined simulations demonstrate thatcosmological simulations with an effective resolution ∼ are feasible with AxioNyx . Assuming maximumvelocities of ∼
100 km/s comparable to those in Figure 4and a FDM mass ∼ − eV, we are able to resolvethe minimum de Broglie wavelength with at least 6 cellswhen allowing for grid resolutions ∼
100 pc in a ∼
10 Mpccomoving box.
ACKNOWLEDGMENTS
We thank Benedikt Eggemeier, Oliver Hahn, ShaunHotchkiss, Emily Kendall, Doddy Marsh, Nathan Mu-soke, and Jan Veltmaat for important discussions, andthe AMReX code development team, especially AnnAlmgren and Guy Moore, for their assistance. Computa-tions described in this work were mainly performed withresources provided by the North-German Supercomput-ing Alliance (HLRN). Additionally, the authors wish toacknowledge the use of New Zealand eScience Infrastruc-ture (NeSI) high performance computing facilities, con-sulting support and/or training services as part of thisresearch. New Zealand’s national facilities are providedby NeSI and funded jointly by NeSI’s collaborator institu-tions and through the Ministry of Business, Innovation &Employment’s Research Infrastructure programme. Weacknowledge the yt [58] toolkit that was used for theanalysis of numerical data. BS acknowledges supportby the Deutsche Forschungsgemeinschaft. JCN acknowl-edges funding by a Julius von Haast Fellowship Awardprovided by the New Zealand Ministry of Business, Inno-vation and Employment and administered by the RoyalSociety of New Zealand. We acknowledge support fromthe Marsden Fund of the Royal Society of New Zealand. [1] M. S. Turner, Phys. Rev. D , 1243 (1983). [2] M. I. Khlopov, B. A. Malomed, and I. B. Zeldovich, Monthly Notices of the Royal Astronomical Society ,
575 (1985).[3] W. H. Press, B. S. Ryden, and D. N. Spergel, Phys. Rev.Lett. , 1084 (1990).[4] W. Hu, R. Barkana, and A. Gruzinov, Physical ReviewLetters , 1158–1161 (2000).[5] S.-J. Sin, Physical Review D , 3650–3654 (1994).[6] V. Sahni and L. Wang, Physical Review D , 103517(2000), arXiv:astro-ph/9910097 [astro-ph].[7] T. Matos, F. S. Guzm´an, and L. A. Ure˜na-L´opez, Clas-sical and Quantum Gravity , 1707 (2000), arXiv:astro-ph/9908152 [astro-ph].[8] F. S. Guzm´an and T. Matos, Classical and QuantumGravity , L9 (2000), arXiv:gr-qc/9810028 [gr-qc].[9] J. Goodman, New Astronomy , 103–107 (2000).[10] P. J. E. Peebles, The Astrophysical Journal ,L127–L129 (2000).[11] L. Amendola and R. Barbieri, Physics Letters B ,192–196 (2006).[12] J.-C. Hwang and H. Noh, Physics Letters B , 1(2009), arXiv:0902.4738 [astro-ph.CO].[13] D. J. E. Marsh, Physics Reports , 1 (2016).[14] J. C. Niemeyer, Progress in Particle and Nuclear Physics , 103787 (2020).[15] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper,and J. March-Russell, Physical Review D , 123530(2010).[16] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys.Rev. D , 043541 (2017).[17] H.-Y. Schive, T. Chiueh, and T. Broadhurst, NaturePhysics , 496 (2014), arXiv:1406.6586.[18] L. M. Widrow and N. Kaiser, The Astrophysical Journal , L71 (1993).[19] C. Uhlemann, M. Kopp, and T. Haugg, Phys. Rev. D , 023517 (2014).[20] D. J. E. Marsh and J. Silk, Monthly Notices of the RoyalAstronomical Society , 2652 (2014).[21] D. J. E. Marsh and A.-R. Pop, Monthly Noticesof the Royal Astronomical Society , 2479 (2015),arXiv:1502.03456 [astro-ph.CO].[22] A. X. Gonz´alez-Morales, D. J. E. Marsh, J. Pe˜narrubia,and L. A. Ure˜na-L´opez,
Monthly Notices of the RoyalAstronomical Society , 1346 (2017), arXiv:1609.05856[astro-ph.CO].[23] T. Bernal, L. M. Fern´andez-Hern´andez, T. Matos, andM. A. Rodr´ıguez-Meza,
Monthly Notices of the Royal As-tronomical Society , 1447 (2018), arXiv:1701.00912[astro-ph.GA].[24] E. Kendall and R. Easther, arXiv e-prints ,arXiv:1908.02508 (2019), arXiv:1908.02508 [astro-ph.CO].[25] H. Deng, M. P. Hertzberg, M. H. Namjoo, andA. Masoumi,
Physical Review D , 023513 (2018),arXiv:1804.05921 [astro-ph.CO].[26] N. Bar, D. Blas, K. Blum, and S. Sibiryakov, PhysicalReview D , 083027 (2018), arXiv:1805.00122 [astro-ph.CO].[27] V. H. Robles, J. S. Bullock, and M. Boylan-Kolchin, Monthly Notices of the Royal Astronomical Society ,289 (2019), arXiv:1807.06018 [astro-ph.CO].[28] A. Burkert, arXiv e-prints , arXiv:2006.11111 (2020),arXiv:2006.11111 [astro-ph.GA].[29] V. Irsic, M. Viel, M. G. Haehnelt, J. S. Bolton, andG. D. Becker, Physical Review Letters , 1 (2017). [30] E. Armengaud, N. Palanque-Delabrouille, C. Y`eche,D. J. E. Marsh, and J. Baur,
Monthly Noticesof the Royal Astronomical Society , 4606 (2017),arXiv:1703.09126 [astro-ph.CO].[31] H.-Y. Schive, T. Chiueh, T. Broadhurst, and K.-W.Huang, The Astrophysical Journal , 89 (2016).[32] Y. Ni, M.-Y. Wang, Y. Feng, and T. Di Matteo,
MonthlyNotices of the Royal Astronomical Society , 5551(2019), arXiv:1904.01604 [astro-ph.CO].[33] X. Li, L. Hui, and G. L. Bryan,
Physical Review D ,063509 (2019), arXiv:1810.01915 [astro-ph.CO].[34] P. Mocz and S. Succi, Physical Review E , 053304(2015), arXiv:1503.03869 [physics.comp-ph].[35] J. Veltmaat and J. C. Niemeyer, Physical Review D ,123523 (2016), arXiv:1608.00802 [astro-ph.CO].[36] M. Nori and M. Baldi, Monthly Notices of the Royal As-tronomical Society , 3935 (2018), arXiv:1801.08144[astro-ph.CO].[37] P. F. Hopkins,
Monthly Notices of the Royal Astronom-ical Society , 2367 (2019), arXiv:1811.05583 [astro-ph.CO].[38] T.-P. Woo and T. Chiueh, The Astrophysical Journal , 850 (2009).[39] B. Schwabe, J. C. Niemeyer, and J. F. Engels, PhysicalReview D (2016), 10.1103/physrevd.94.043513.[40] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala,M. Boylan-Kolchin, A. Fialkov, and L. Hernquist,Monthly Notices of the Royal Astronomical Society ,4559 (2017).[41] J. Veltmaat, J. C. Niemeyer, and B. Schwabe, PhysicalReview D (2018), 10.1103/physrevd.98.043509.[42] P. Mocz, A. Fialkov, M. Vogelsberger, F. Becerra, M. A.Amin, S. Bose, M. Boylan-Kolchin, P.-H. Chavanis,L. Hernquist, L. Lancaster, F. Marinacci, V. H. Robles,and J. Zavala, Phys. Rev. Lett. , 141301 (2019).[43] J. Veltmaat, B. Schwabe, and J. C. Niemeyer, Phys.Rev. D , 083518 (2020).[44] A. S. Almgren, . B. Bell, M. J. Lijewski, Z. Luki´c, andE. V. Andel, The Astrophysical Journal , 39 (2013).[45] A. Almgren, V. Beckner, J. Blaschke, C. Chan, M. Day,B. Friesen, K. Gott, D. Graves, M. Katz, A. Myers,T. Nguyen, A. Nonaka, M. Rosso, S. Williams, W. Zhang,and M. Zingale, (2019), 10.5281/zenodo.2555438.[46] F. S. Guzm´an and L. A. Ure˜na L´opez, Phys. Rev. D ,124033 (2004).[47] B. Eggemeier and J. C. Niemeyer, Phys. Rev. D ,063528 (2019).[48] D. Levkov, A. Panin, and I. Tkachev, Physical ReviewLetters (2018), 10.1103/physrevlett.121.151301.[49] N. Musoke, S. Hotchkiss, and R. Easther, Phys.Rev. Lett. , 061301 (2020), arXiv:1909.11678 [astro-ph.CO].[50] J. C. Niemeyer and R. Easther, (2019), arXiv:1911.01661[astro-ph.CO].[51] F. Edwards, E. Kendall, S. Hotchkiss, and R. Easther,JCAP , 027 (2018), arXiv:1807.04037 [astro-ph.CO].[52] R. L¨ohner, Computer Methods in Applied Mechanics andEngineering , 323 (1987).[53] G. L. Bryan, M. L. Norman, B. W. O’Shea, T. Abel,J. H. Wise, M. J. Turk, D. R. Reynolds, D. C. Collins,P. Wang, S. W. Skillman, B. Smith, R. P. Harkness,J. Bordner, J.-h. Kim, M. Kuhlen, H. Xu, N. Goldbaum,C. Hummels, A. G. Kritsuk, E. Tasker, S. Skory, C. M.Simpson, O. Hahn, J. S. Oishi, G. C. So, F. Zhao, R. Cen, t a v e r a g e r e l a t i v e e rr o r FIG. 8. Spatially averaged error of the final state as a functionof the time step for the two solvers. Each point is a separatesimulation with a fixed time step. The fitted slopes of the lineaccurated reproduce the expected scalings. At sixth orderwith ∆ t > − the error is dominated by numerical noisedue to finite precision floating point numbers.Y. Li, and The Enzo Collaboration, apjs , 19 (2014),arXiv:1307.2265 [astro-ph.IM].[54] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zin-gale, D. Q. Lamb, P. MacNeice, R. Rosner, J. W. Tru-ran, and H. Tufo, Astrophysical Journal, Supplement , 273 (2000).[55] M. Mina, D. F. Mota, and H. A. Winther, “Scalar: anamr code to simulate axion-like dark matter models,”(2019), arXiv:1906.12160 [physics.comp-ph].[56] W. Hu and D. J. Eisenstein, The Astrophysical Journal , 497 (1998).[57] D. J. E. Marsh and J. C. Niemeyer, Phys. Rev. Lett. ,051103 (2019).[58] M. J. Turk, B. D. Smith, J. S. Oishi, S. Skory, S. W.Skillman, T. Abel, and M. L. Norman, The As-trophysical Journal Supplement Series , 9 (2011),arXiv:1011.3514 [astro-ph.IM].[59] D. J. E. Marsh, (2016), arXiv:1605.05973 [astro-ph.CO]. Appendix A: Code Verification
We tested the code with a suite of increasingly chal-lenging computations.
1. Static potential tests
Setting the gravitational potential V = 0, an initiallyGaussian FDM overdensity broadens via diffusion. Thedynamics can be solved exactly, ρ ( x, t ) = 1 (cid:112) πσ ( t ) exp (cid:20) − ( x − x ( t )) σ ( t ) (cid:21) , (A1) σ ( t ) = σ (cid:34) (cid:18) ( t − t ) (cid:126) mσ (cid:19) (cid:35) , (A2) r e l . c h a n g e [ % ] M W K E
FIG. 9. Relative changes of global quantities over time duringthe simulation of a single soliton. Due to discreteness errors,it is slightly ringing with a mass dependent frequency as canbe inferred from the potential energy W and kinetic energy K (see, e.g., [39] for the definition of conserved quantities).Total mass and energy are well conserved even in the depictedlow resolution run. where the initial width σ ( t ) = σ and x ( t ) = 0. Themaximum relative errors are smaller then 10 − until σ ( t f ) = 2 σ even for a low resolution grid spacing,∆ x = 0 . σ . Placing the Gaussian over-density in anexternal harmonic oscillator potential V ( x ) = − (cid:20) (cid:126) mσ (cid:21) x , (A3)we likewise verify that σ ( t f ) = σ with similar numericalerror. If the overdensity is initially displaced by x fromthe origin, as expected the center oscillates within thepotential well as x ( t ) = x cos (cid:20) ( t − t ) (cid:126) mσ (cid:21) . (A4)We confirm the ideal temporal convergence of both nu-merical schemes by comparing spatially averaged devi-ations to the analytical solution after one oscillation,shown in Figure 8. Note that ∆ t = 10 − correspondsto the maximum time step allowed by stability criteriafor a comparable finite difference simulation.
2. Isolated soliton
As the ground state of the Schr¨odinger-Poisson sys-tem, a soliton is stationary in time. However, any initialor numerically evolved configuration can only approachthe analytical ground state. Figure 9 shows that smallperturbations result in a periodic ringing of the soliton,which contracts and expands. While we recover the cor-0rect quasinormal frequency [41] ν = 10 . (cid:18) ρ c M (cid:12) kpc − (cid:19) / Gyr − , (A5)the relative change in potential and kinetic energy de-pends on the size of the perturbation and thus on reso-lution. Here, ρ c is the central soliton density.
3. Boson Star Condensation
The condensation of a Boson star from random Gaus-sian initial fluctuations, as first investigated in [48], pro-vides a stringent test of the code. The initial conditionswere set up on a 64 grid with a side length of L = 30 λ dB .The overall amplitude of the fluctuations was scaled suchthat the total mass in the simulation box is N = 10; werefer the reader to the original paper for details [48]. Fig-ure 10 shows the expected linear increase in maximumamplitude after one condensation time. m a x . a m p li t u d e datafit FIG. 10. Linear growth of the maximum amplitude due tothe condensation of a Boson star. scale factor10 m a x . o v e r d e n s i t y high kmedium klow k FIG. 11. Time evolution of different self-gravitating modes inan expanding background. Modes with wave numbers largerthan the Jeans scale are stabilized by wave coherence effects.
4. Linear mode evolution
A one-dimensional, self-gravitating, sinusoidal over-density in a universe containing 100% FDM obeys equa-tion 5a with the RHS set to zero. We observe the ex-pected evolution of different growing modes [13] δ + ( k, a ) = 3˜ k sin (cid:16) ˜ k (cid:17) + (cid:20) k − (cid:21) cos (cid:16) ˜ k (cid:17) , (A6)where ˜ k = k/ (cid:112) √ amH / (cid:126) ∝ k/k J ( a ) is inversely propor-tional to the Jeans scale k J ( a ). While modes with wavenumber k < k J grow linearly with scale factor a , modeswith k > k J oscillate in time. Intermediate modes justbelow k J oscillate until they cross the Jeans scale afterwhich they start to grow linearly. All three scenarios canbe seen in Figure 11.
5. Spherical collapse of fuzzy dark matter
We investigated the spherical collapse of a Gaussianover-density with initial maximum density contrast δ =0 . σ region of theGaussian perturbation. We thus scale the over-densityby changing the box size. The Jeans mass below whichcollapse is suppressed can be approximated as [59] M J ( a ) =3 . × (cid:16) m − eV (cid:17) − / (cid:18) Ω m h . (cid:19) / × h − a − / M (cid:12) (cid:39) . × a − / M (cid:12) , (A7)where we used Ω m = 1, h = 0 . m = 2 . × − eV.From Figure 12 we see that collapse occurs once the massinside the box starts to exceed M J ( a ). scale factor10 m a x . F D M o v e r d e n s i t y a c <0.01 a c =0.0123 a c =0.1965 a c >1.0 FIG. 12. Evolution of the maximum FDM over-density duringspherical collapse beginning as soon as the mass inside thesimulation box exceeds the Jeans mass at a scale factor a cc