# Cosmological simulations of mixed ultralight dark matter

Alex Laguë<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA
19104-6396, USA  Dunlap Institute for Astronomy & Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada  Canadian Institute
for Theoretical Astrophysics, University of Toronto, 60 St. George Street,
Toronto, ON M5S 3H8, Canada  David A. Dunlap Department of Astronomy &
Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S
3H4, Canada  Bodo Schwabe  Departamento de Física Teórica, Universidad de
Zaragoza, 50009 Zaragoza  Institut für Astrophysik, Universität Göttingen,
Germany  Renée Hložek  Dunlap Institute for Astronomy & Astrophysics,
University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada
David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto,
50 St. George Street, Toronto, ON M5S 3H4, Canada  David J. E. Marsh
Theoretical Particle Physics and Cosmology, King’s College London, Strand,
London, WC2R 2LS, United Kingdom  Keir K. Rogers  Dunlap Institute for
Astronomy & Astrophysics, University of Toronto, 50 St. George Street,
Toronto, ON M5S 3H4, Canada

###### Abstract

The era of precision cosmology allows us to test the composition of the dark
matter. Mixed ultralight or fuzzy dark matter (FDM) is a cosmological model
with dark matter composed of a combination of particles of mass $m\leq
10^{-20}$ eV, with an astrophysical de Broglie wavelength, and particles with
a negligible wavelength sharing the properties of cold dark matter (CDM). In
this work, we simulate cosmological volumes with a dark matter wave function
for the ultralight component coupled gravitationally to CDM particles. We
investigate the impact of a mixture of CDM and FDM in various proportions
($0\%,\;1\%,\;10\%,\;50\%,\;100\%$) and for ultralight particle masses ranging
over five orders of magnitude ($2.5\times 10^{-25}\;\mathrm{eV}-2.5\times
10^{-21}\;\mathrm{eV}$). To track the evolution of density perturbations in
the non-linear regime, we adapt the simulation code AxioNyx to solve the CDM
dynamics coupled to a FDM wave function obeying the Schrödinger-Poisson
equations. We obtain the non-linear power spectrum and study the impact of the
wave effects on the growth of structure on different scales. We confirm that
the steady-state solution of the Schrödinger-Poisson system holds at the
center of halos in the presence of a CDM component when it composes 50% or
less of the dark matter but find no stable density core when the FDM accounts
for 10% or less of the dark matter. We implement a modified friends-of-friends
halo finder and find good agreement between the observed halo abundance and
the predictions from the adapted halo model axionHMCode.

††preprint: APS/123-QED

##  I Introduction

The cold dark matter (CDM) model explains the formation of cosmic structures
exceptionally well on large scales. Observations suggest that dark matter may
consist of multiple components, with varying properties affecting the growth
of structure on different scales [1]. Among the dark matter candidates which
cluster distinctly from CDM are very light scalar bosons with negligible non-
gravitational interactions known as fuzzy dark matter (FDM) [2, 3]. Due to
their low particle mass, FDM condensates are subject to coherence effects
which impede clustering on scales below their de Broglie wavelength $\hbar/mv$
(for a particle of mass $m$ and velocity $v$), but behave as CDM on large
scales. Ideal candidates for FDM from particle physics are ultralight axions
which arise naturally in quantum chromodynamics [4, 5] and in high energy
physics extensions to the Standard Model such as string theory [6, 7, 8, 9].
In this scenario, it is postulated that they form in a plenitude of
$\mathcal{O}(100)$ axion fields with logarithmically distributed masses [6].
Other extensions of the Standard Model not based on string theory, such as the
pi-axiverse [10], can also lead to a set of ultralight bosons. In both cases,
the lightest of these axions could have a mass lower than $10^{-19}$ eV,
implying the existence of particles with de Broglie wavelengths on galactic
scales. In the high energy physics models where they arise, the ultralight
particles’ relic density, which we will denote $\Omega_{\mathrm{FDM}}$, is not
necessarily equal to the total dark matter density of the Universe. For
instance, Ref. [11] finds an axion with a mass of $2.5\times 10^{-22}$ eV and
a cosmological density about one-tenth of the predicted total dark matter
density from the $\Lambda$CDM model. Throughout this work, we make the
simplifying assumption that only one of the particles produced in this
scenario is ultralight while the others have negligible de Broglie
wavelengths. In this case, we group their combined relic density into
$\Omega_{\mathrm{CDM}}$ since they can be modeled as CDM by virtue of the
Schrödinger-Vlasov correspondence [12].

Figure 1: Side-by-side slice plot of CDM and FDM in the same 1 Mpc/$h$ box
(comoving) at redshift $z=4$. The dark matter is composed of 10% FDM, with the
rest being CDM. The inner region illustrates the difference in clustering
between the two dark matter species inside the central halo.

With the total relic density of the non-ultralight particles captured by
$\Omega_{\mathrm{CDM}}$, we define the FDM fraction as

$\displaystyle
f\equiv\frac{\Omega_{\mathrm{FDM}}}{\Omega_{\mathrm{CDM}}+\Omega_{\mathrm{FDM}}}.$
(1)

The mixed cold and fuzzy dark matter model has been studied in numerical
simulations for single halos [13, 14], but cosmological $N$-body and
hydrodynamical simulations (accounting for FDM dynamics) up to date have
assumed $f=1$ [15, 16, 17, 18]. In this work, we investigate the full non-
linear behaviour of FDM particles with cosmological initial conditions while
foregoing this assumption. Previews of the simulation with multiple dark
matter components are shown in Figs. 1 and 2.

Constraints on the FDM particle mass and fraction have been reached with CMB
[19, 20, 3], galaxy clustering [21, 3], galaxy weak lensing [22], and
Lyman-$\alpha$ forest data [23, 24, 25]. With the Lyman-$\alpha$ forest,
accounting for the full evolution of the FDM wave function is crucial to
arrive at unbiased constraints in the presence of baryons [26]. While
constraints using the Lyman-$\alpha$ forest were obtained for mixed and pure
FDM scenarios, only the case where $f=1$ has been verified with $N$-body
simulations which included a non-linear treatment of the wave effects of FDM
[27]. Studies of the density profiles of ultra-faint dwarf galaxies indicate a
preference for a FDM mass $m\sim 3.7-5.6\times 10^{-22}$ eV [28] while dwarf
spheroidal galaxies suggest $m\leq 1.1\times 10^{-22}$ eV [29], and
Lyman-$\alpha$ forest results give $m\geq 3.8\times 10^{-21}$ eV [23] and
$m>2\times 10^{-20}$ eV [25], respectively. Furthermore, the suppression of
smaller halos conflicts with the subhalo mass function for ultralight dark
matter masses below $m\sim 2.1\times 10^{-21}$ eV [30]. It is also in tension
with measurements of stellar streams in the Milky Way [31]. However, these
arguments are based on simulations of a single dark matter component.
Investigating the behaviour and scaling relations of cored density profiles
(which are the ground state of the wave function and often referred to as
solitons) in more general scenarios is crucial to establish if an internal
tension with FDM exists when considering the full phenomenology of the
axiverse. At the higher mass end, the lack of detection of black hole super-
radiance constrains ultralight masses between $10^{-19}\;\mathrm{eV}\lesssim
m\lesssim 10^{-16}\;\mathrm{eV}$ and $10^{-13}\;\mathrm{eV}\lesssim m\lesssim
10^{-11}\;\mathrm{eV}$ [32, 33]. The latter constraints still hold in the
presence of multiple axion fields.

Figure 2: Density slices of 1 Mpc/$h$ simulation boxes (comoving) at redshift
$z=4$. The three panels (starting from the left) illustrate the total (CDM
plus FDM) density around the central halo for a simulation of pure CDM, an
even mixture of CDM and FDM, and of pure FDM for a FDM mass of $2.5\times
10^{-22}$ eV.

Axions are described by a scalar field $\phi$ with mass $m$ obeying the Klein-
Gordon equation [34, 35]

$\displaystyle\Box\phi-\frac{\partial V_{\phi}}{\partial\phi}=0,$ (2)

where the box denotes the d’Alembertian operator and $V_{\phi}$ is the field
potential. The general form of this potential is the periodic

$\displaystyle
V_{\phi}=m^{2}f_{a}^{2}\left[1-\cos\left(\frac{\phi}{f_{a}}\right)\right],$
(3)

where $f_{a}$ is the field decay constant. In the limit of small
displacements, $\phi\ll f_{a}$, the potential can be expressed as a quadratic:
$V_{\phi}=m^{2}\phi^{2}/2$. In the non-relativistic limit, the field has the
ansatz [36]

$\displaystyle\phi=\frac{\hbar}{m\sqrt{2}}\left(\psi
e^{-imt/\hbar}+\psi^{*}e^{imt/\hbar}\right),$ (4)

where $|\psi|^{2}$ is proportional to the FDM density in that limit.

##  II Mixed Dark Matter Simulations

###  II.1 Wave Function Evolution

Following the non-relativistic approximation of Eq. 4, the dynamics of FDM are
governed by the Schrödinger equation

$\displaystyle i\hbar\frac{\partial\psi}{\partial
t}=\left(-\frac{\hbar^{2}}{2ma^{2}}\nabla^{2}+m\Phi_{\rm N}\right)\psi\,,$ (5)

where $m$ is the FDM particle mass and $\Phi_{\rm N}$ is the Newtonian
gravitational potential. In a mixed dark matter cosmology, the evolution of
the FDM component is coupled to CDM via the Poisson equation

$\displaystyle\nabla^{2}\Phi_{\rm N}=\frac{4\pi
G}{a}\left(\rho_{\mathrm{CDM}}+\rho_{\mathrm{FDM}}-\bar{\rho}_{\mathrm{tot}}\right)\,,$
(6)

where $\rho_{\mathrm{FDM}}=|\psi|^{2}$ and $\bar{\rho}_{\mathrm{tot}}$ is the
mean of the total dark matter density.

To solve the above system of equations, we use AxioNyx [13] which is an
extension of the cosmological simulation code Nyx [37]. While the $N$-body
solver of the Nyx code is used to evolve the CDM component, the pseudo-
spectral solver described in Ref. [13] is used for the FDM component. This
method has been shown to solve the Schrödinger-Poisson system accurately and
to resolve the small-scale features of the wave function [13]. However, the
latter is only possible if the simulation grid spacing is sufficiently small
to resolve the de Broglie wavelength throughout the simulation. This means
that the Schrödinger-Poisson solver is subject to the time-step criterion [38,
39]

$\displaystyle\Delta t\leq 4\min\left[\frac{m}{\hbar}\frac{(\Delta
x)^{2}}{2\pi},\frac{\pi\hbar}{4m|\Phi_{\rm N}|_{\mathrm{max}}}\right]\,,$ (7)

because the time step needs to capture the coherence time scale
$t_{c}\sim(\Delta x)^{2}m/\hbar\sim\hbar/(mv^{2})$ of the field. The
constraints in both spatial and temporal resolution are the reason why FDM
simulations cannot reach the spatial extent of pure CDM $N$-body simulations
and most FDM simulations are stopped around redshift $z\lesssim 3$.

Multiple algorithms have been adapted or modified to solve the Schrödinger-
Poisson system of equations. Gamer [40] was notably one of the first followed
by AREPO [38, 41], ENZO [42], NYX [43], GIZMO [18] and RAMSES [44]. A Python-
based pseudo-spectral solver named PyUltralight [39] has also been recently
developed. Some codes have been extended to simulate multiple axion fields
[45, 14, 46]. There are also algorithms which take advantage of the fact that
the Schrödinger-Poisson system admits a steady-state solution at the center of
virialized structures. These express the wave function as a sum of eigen-
functions which can be evolved forward in time at a reduced computational cost
[47].

We exploit the efficiency of the pseudo-spectral solver to numerically evolve
the FDM component and, given the size of the boxes and resolution used in this
work, we stop the evolution of the wave function at redshift $z=4$ before all
the modes become nonlinear. Since the pseudo-spectral solver relies on
periodic boundary conditions, it cannot be used on higher levels of refinement
where instead, a finite-difference solver is generally used. Finite-difference
is used to approximate the Laplacian appearing in the Schrödinger-Poisson
equations while, in the pseudo-spectral method, the Laplacian is computed
without numerical approximations. This difference in the interplay of the two
solvers can result in the lagging of the finite-difference solver with respect
to the pseudo-spectral solver. For this reason, we postpone the implementation
and testing of adaptive mesh refinement (AMR) in mixed dark matter simulations
to future work.

###  II.2 Initial Conditions

Given the FDM wave effects, the CDM and FDM components will not have the same
initial distributions of densities and velocities. One also has to use higher-
order perturbation theory to avoid the formation of transient features in the
simulations [48]. A second-order Lagrangian perturbation theory (LPT) scheme
for mixed ultralight axions was developed in Ref. [49]. The scale-dependent
suppression in the FDM growth factor reproduces the effects of the wave
diffusion at high redshifts since wave effects are captured at the linear
level by the effective sound speed [50, 51]

$\displaystyle c^{2}_{s,\mathrm{eff}}=\frac{\hbar^{2}k^{2}}{4m^{2}a^{2}},$ (8)

where $k$ is the comoving Fourier mode and $a$ is the cosmological scale
factor. The modified LPT approach is, therefore, well suited for generating
initial conditions for FDM cosmological simulations at high redshift (here we
choose $z_{\mathrm{ini}}=100$ to be well into the linear regime). The main
difference for the FDM component is that the velocity field at first order
obeys

$\displaystyle\mathbf{v}(\mathbf{k})=-\mathcal{H}\frac{d\ln D(k,a)}{d\ln
a}\frac{\mathbf{k}}{k^{2}}\delta(\mathbf{k}),$ (9)

where $\mathcal{H}=aH$, $\delta(k)$ is the Fourier transform of the
overdensity field $\delta\equiv\rho/\bar{\rho}-1$, and $D$ is the linear
growth factor of the density perturbations (for the CDM component, the growth
factor is scale-independent). The density field used is to compute the
velocities is the total density, including both components. Therefore, the
presence of FDM also affects the initial velocity of CDM particles by
partially suppressing the initial power spectrum on small scales. In Ref.
[49], it was shown that one can approximate

$\displaystyle D(k,a)\approx L(k)D_{\mathrm{CDM}}(a),$ (10)

where the prefactor $L$ is a monotonically decreasing function which
asymptotically tends to zero on small sales. We use the public code MUSIC [52]
along with the modified Boltzmann code axionCAMB [19]. We create a transfer
function input file for each of the components and then generate initial
positions and velocities for each, given the same total gravitational
potential. The suppressed FDM transfer function ensures that the FDM particles
are not given inconsistent velocities. The resulting initial density and
velocity fields for both the CDM and FDM components of the simulations are
shown for $f=0.1$ and $m=2.5\times 10^{-22}$ eV in Fig 3. The difference in
the number and placements of the velocity field arrows is due to the fact that
the velocity field of the FDM is evaluated on a grid rather than with
particles. To initialize the wave function, we use the Madelung change of
variables $\psi=Re^{i\theta}$ where the magnitude and phase can be obtained
with

$\displaystyle R(\mathbf{x})$
$\displaystyle\equiv\sqrt{\frac{\rho_{\mathrm{FDM}}(\mathbf{x})}{m}},$ (11)
$\displaystyle\nabla\theta(\mathbf{x})$
$\displaystyle\equiv\frac{m}{\hbar}\mathbf{v}_{\mathrm{FDM}}(\mathbf{x}).$
(12)

To obtain the phase of the wave function, Eq. (12) can be solved in Fourier
space given the Fourier transform of the velocity field which was calculated
with Eq. (9). In Fig. 3, we can observe the difference in clustering between
the two species at high redshift. Many simulations have used the modified
initial conditions of FDM combined with CDM-like evolution to approximate the
behaviour of FDM at a lower computational cost (by not solving the Schrödinger
equation). This has been known as the warm dark matter (WDM) or classical FDM
approach [16, 53]. We will test the validity of this approximation for low FDM
concentrations in Sec. III.

Figure 3: Initial density field for the CDM (top) and FDM (bottom) components
of our simulation box 3 (see Table 1) at $z=100$. The initial velocity fields
are plotted and rescaled to the maximal velocity for each component. The ratio
of the FDM to CDM velocity magnitudes goes to one on large scales.

##  III Results

We simulate the evolution of the CDM and FDM components from the initial
redshift $z_{\mathrm{ini}}=100$. We generate a set of 12 simulation boxes
listed in Table 1. First, we fix the mass at $m=2.5\times 10^{-22}$ eV and
select five FDM fractions: $f\in\\{0,\;0.01,\;0.1,\;0.5,\;1\\}$. We then
consider four other masses $m\in
2.5\times\\{10^{-25},\;10^{-24},\;10^{-23},\;10^{-21}\\}$ eV, but keep the
fraction at $f=0.1$ for these runs. We also create two simulations where the
initial conditions account for 1% and 10% FDM, but the system is evolved using
only CDM particles and a simulation with twice the resolution using a grid of
$1024^{3}$ points for a convergence analysis (and Appendix B). Snapshots of
the total matter density in boxes 1, 4, and 5 at redshift $z=4$ are displayed
in Fig. 2. The visualizations are generated using the analysis code yt [54].

Box | FDM | FDM | Box |  $N$ Grid/
---|---|---|---|---
Number | Mass (eV) | Fraction | Length | Particles
1 | $2.5\times 10^{-22}$ | $0.00$ | 1.0 Mpc$/h$ | $512^{3}$
2 | $2.5\times 10^{-22}$ | $0.01$ | 1.0 Mpc$/h$ | $512^{3}$
3 | $2.5\times 10^{-22}$ | $0.10$ | 1.0 Mpc$/h$ | $512^{3}$
4 | $2.5\times 10^{-22}$ | $0.50$ | 1.0 Mpc$/h$ | $512^{3}$
5 | $2.5\times 10^{-22}$ | $1.00$ | 1.0 Mpc$/h$ | $512^{3}$
6 | $2.5\times 10^{-21}$ | $0.10$ | 0.3 Mpc$/h$ | $512^{3}$
7 | $2.5\times 10^{-23}$ | $0.10$ | 1.0 Mpc$/h$ | $512^{3}$
8 | $2.5\times 10^{-24}$ | $0.10$ | 10.0 Mpc$/h$ | $512^{3}$
9 | $2.5\times 10^{-25}$ | $0.10$ | 30.0 Mpc$/h$ | $512^{3}$
10 | $2.5\times 10^{-22}$ | $0.10$ | 1.0 Mpc$/h$ | $1024^{3}$
11 | FDM ICs | 0.01 | 1.0 Mpc$/h$ | $512^{3}$
12 | FDM ICs | 0.10 | 1.0 Mpc$/h$ | $512^{3}$

Table 1: List of completed simulations with FDM masses and fractions. The box
size is given in comoving coordinates. The number of particles for the CDM
component is the same as the number of grid points for the FDM in all cases.
FDM ICs under particle mass indicates that the particles were evolved as CDM
with modified initial conditions.

###  III.1 Non-Linear Matter Power Spectrum

Figure 4: Power spectra in the simulation boxes as a function of redshift and
FDM fraction. (Top) Combined (CDM+FDM) power spectra normalized by the time-
dependent growth factor. The dashed line represents the prediction from linear
theory for pure CDM. (Bottom) Ratio of the combined power spectra with respect
to the non-linear pure CDM power spectrum at the same redshift.  Figure 5:
Combined (CDM+FDM) power spectra normalized by the time-dependent growth
factor in the simulation boxes as a function of redshift and FDM mass for a
fraction of $f=0.1$. We note the formation of non-linear structures at
redshift $z=7$ for all masses which is not found in simulations where $f=1$.

We extract the total power spectra from our simulations using the public code
nbodykit [55]. The spectra and relative difference of the FDM and CDM power
spectra for a series of FDM fractions is shown as a function of redshift in
Figs 4-5. We have factored out the homogeneous (large-scale) growth factor
$D(z)\propto 1/(1+z)$ in the matter-dominated Universe to highlight the growth
of structure on small scales. For high FDM fractions (redder curves), the FDM
and CDM spectra are relatively close together at high redshift. In this case,
the CDM follows the distribution of the FDM and will be similarly suppressed
on small scales. However, the converse does not hold. If the CDM is the main
component of the dark matter, the FDM still does not cluster on scales below
its Jeans scale (where the FDM sound speed equates to the gravitational
potential) [34], which in the pure FDM case is given by

$\displaystyle k_{\rm
J}=66.5a^{1/4}\left(\frac{m}{10^{-22}\;\mathrm{eV}}\right)^{1/2}\left(\frac{\Omega_{\mathrm{FDM}}h^{2}}{0.12}\right)^{1/4}\;\mathrm{Mpc}^{-1}.$
(13)

This is visible in the pale blue and grey curves in the $z=15$ panel of Fig.
4.

In the bottom row of Fig. 4, we also plot the ratio of the total matter power
spectrum to the non-linear matter power spectrum of the pure CDM case ($f=0$).
As observed in Ref. [21], even a FDM fraction of a few percents leads to a
strong suppression of clustering. This explains why the 10% FDM fraction (pink
line) is damped by over 50% at redshift $z=15$. We note however, that the FDM
rapidly falls into the CDM potential well, and the total non-linear power
spectrum approaches that of the pure CDM case as structure forms at lower
redshifts. Comparing the results shown in Fig. 5 with the simulations of Ref.
[56], we can see that non-linearities form sooner in the $f=0.1$ case than in
the $f=1$ case as a power spectrum excess (compared to the linear prediction)
is visible at redshift $z\sim 7$. This is mostly attributable to the lower
level of suppression in the linear spectrum for $f<1$, but also to the scale-
dependent structure growth of FDM which delays the in-fall of FDM in halos
[49].

Figure 6: Impact of wave effects on the power spectrum at redshift $z=4$. We
increase the amplitude of $\mathcal{Q}$ by an order of magnitude for $f=0.01$
to compare the shape of the curves across FDM fractions. The late-time FDM
wave effects not captured by the modified initial conditions are more apparent
when $f=0.1$. The dotted black line denotes the Jeans scale for a mass of
$2.5\times 10^{-22}$ eV.

To measure the impact of the wave effects and differentiate them from the
structure suppression caused by the modified initial conditions, we run two
additional simulations using the WDM approximation as used in Refs [16, 53].
The WDM approach consists in modifying the initial conditions to account for
the presence of FDM but to evolve the system for $z<z_{\mathrm{ini}}$ using
only the CDM dynamics. For this, we use the initial conditions of boxes 2 and
3 with $f=0.01-0.10$, we label the power spectra of such simulations
$P^{\mathrm{w}}_{f}$. Taking the ratio of this power spectrum to the power
spectrum of the simulation boxes which account for both the modified initial
conditions and the FDM dynamics, we have

$\displaystyle\mathcal{Q}(k)\equiv\frac{P_{f}(k)}{P^{\mathrm{w}}_{f}(k)}-1.$
(14)

The role of the function $\mathcal{Q}(k)$ is to isolate the impact of the
late-time effect of the wave behaviour of FDM. If this function is exactly
zero, this tells us the wave effects play no role in clustering and that the
all the suppression observed in Figs 4-5 is due solely to the initial
conditions. We observe increased suppression of the non-linear matter power
spectrum when accounting for the wave effects, as shown in Fig. 6. We note
that the shape of the function $\mathcal{Q}$ remains the same as we increase
the FDM fraction, and only its amplitude changes. Our results agree with the
finding of Ref. [15] with the added property that difference between the WDM
and FDM treatments becoming negligible at low fractions. In other words, we
find $\mathcal{Q}(k)\to 0$ as $f\to 0$ roughly linearly.

###  III.2 Density Profiles

A well-known prediction of FDM models is the formation of solitonic cores at
the center of halos [40]. We can obtain the shape of this core by noting that
the Schrödinger-Poisson system has an equilibrium solution with a complex
phase. The wave function takes the form $\psi(\mathbf{x},t)=e^{-i\gamma
t/\hbar}\phi(\mathbf{x})$, where $\gamma$ is a constant. Assuming spherical
symmetry (with $r=|\mathbf{x}|$), the system of equations (5)-(6) becomes [57,
58]

$\displaystyle\frac{\partial^{2}(r\phi)}{\partial r^{2}}$
$\displaystyle=2r\left(\frac{m^{2}}{\hbar^{2}}V-\frac{m}{\hbar^{2}}\gamma\right)\phi,$
(15) $\displaystyle\frac{\partial^{2}(rV)}{\partial r^{2}}$
$\displaystyle=4\pi Gr\left(\phi^{2}+\rho_{\mathrm{CDM}}\right).$ (16)

We fit the density profile of the FDM by solving for the equilibrium solution
given the shape of the CDM density profile and the ratio of their central
densities. We develop a more stable alternative to the shooting method used in
previous studies to find the eigenvalue $\gamma$ [57] which we describe in
Appendix A. We confirm that the equilibrium state of the scalar function
$\phi$ does satisfy the Schrödinger-Poisson system but with a different
eigenvalue than the one found for the zero-node solution in the pure FDM case.

Figure 7: CDM and FDM density profiles in simulation box 4 with 50% FDM (see
Table 1). The dark dashed line represents a traditional soliton profile fit to
the FDM density and the green line denotes the fit found by solving Eq. (16)
in the presence of CDM. The arrows denote the half-density radius of the
curves of the corresponding color. Figure 8: CDM and FDM density profiles in
simulation box 3 with 10% FDM (see Table 1). The dark dashed line represents a
traditional soliton profile fit to the FDM density and the green line denotes
the fit found by solving Eq. (16) in the presence of CDM. The arrows denote
the half-density radius of the curves of the corresponding color.

As with a baryon-FDM mixture, we observe that a secondary contributor to the
potential of the FDM causes the core radius to shrink [58]. In Fig. 7 and Fig.
8, we display the CDM and FDM density alongside the canonical soliton fit of
Ref. [40]. The fit can be obtained with two free parameters giving

$\displaystyle\rho_{\mathrm{sol}}(r)=\frac{\rho_{0}}{\left[1+0.091(r/r_{c})^{2}\right]^{8}},$
(17)

where $\rho_{0}$ is the central density and where $r_{c}$ is the radius at
which the density reaches half of its central value. The FDM mass and the two
parameters of the fit obey the scaling relation

$\displaystyle r_{c}=1.0$ $\displaystyle\left(\frac{\rho_{0}}{3.1\times
10^{15}\;M_{\odot}/\mathrm{Mpc}^{3}}\right)^{1/4}$
$\displaystyle\times\left(\frac{m}{2.5\times
10^{-22}\;\mathrm{eV}}\right)^{1/2}\;\mathrm{kpc}.$ (18)

We find in the mixed CDM-FDM case that the relation for the core radius no
longer holds. In the case where $f=0.5$, the FDM dominates the inner region of
the halo and the soliton fit still provides a good description of the inner
FDM density with a slightly different slope at high $r$. However, when CDM
composes 90% of the dark matter, we do not detect a proper soliton core. We
test to see if the core radius falls below the resolution of the simulation by
solving the time-independent ground state of the Schrödinger-Poisson system in
the presence of CDM. The CDM creates a steeper gravitational potential
gradient on the FDM wave function causing it to be radially compressed. In our
simulations, we find that the FDM density profile’s inner slope is too steep
to be well modeled by a proper soliton core for the range of radii we can
resolve. Our results agree with previous studies on mixed dark matter using
spherical collapse, which found no core formation for fractions below
$f\lesssim 0.3$ [13].

On the outskirts of the halos, we find that the total density profile
approaches the Navarro-Frenk-White (NFW) fit [59]:

$\displaystyle\rho(r)=\frac{\rho_{s}}{(r/r_{s})(1+r/r_{s})^{2}},$ (19)

where $r_{s}$ is the NFW scale radius and $\rho_{s}$ is its characteristic
density. The NFW profile has a diverging mass as $r\to\infty$, which is
unphysical. It is customary to denote the limit of the halo with its virial
radius, which is defined as a function of the virial overdensity
$\Delta_{\mathrm{vir}}$ given by [60]

$\displaystyle\Delta_{\rm
vir}=18\pi^{2}+82\left[\Omega_{\mathrm{m}}(z)-1\right]-39\left[\Omega_{\mathrm{m}}(z)-1\right]^{2}$
(20)

and which we can relate to the halo mass through

$\displaystyle M(r_{\rm vir})=\frac{4\pi}{3}r_{\rm vir}^{3}\Delta_{\rm
vir}\bar{\rho}_{\rm tot}.$ (21)

We will use this definition in the following section.

Figure 9: Tri-dimensional density isocontours around a CDM-dominated sub-halo
with an FDM mass of $m=2.5\times 10^{-22}$ eV. Despite being in a simulation
with $f=0.1$, $f\ll 0.1$ within the halo, indicating it is close to the cut-
off where no FDM is clustered. The three isocontour panels show the separate
components and the total density. The lower left panel shows the measured
density of CDM and FDM as a function of the distance from the sub-halo center
up to 40 kpc which corresponds to the white dashed circle in the other three
panels. The colored text indicates the density of the isocontours in units of
$M_{\odot}/\mathrm{Mpc}^{3}$.

A key assumption in the halo model of Ref. [61] is that CDM could form bound
structures which contain little to no FDM. This can happen when a halo has a
radius smaller than the halo Jeans radius ($r_{\mathrm{hJ}}$) which is defined
as the radius where

$\displaystyle\rho(r_{\mathrm{hJ}})\approx\frac{\bar{\rho}_{\rm
tot}\Delta_{\rm vir}c^{2}r_{\rm vir}}{3f(c)r_{\mathrm{hJ}}},$ (22)

where $c$ is the concentration of the halo which is defined as the ratio of
the virial radius to the scale radius (from the NFW fit) $c\equiv r_{\rm
vir}/r_{s}$ and $f(x)=\ln(x+1)-\frac{x}{x+1}$. However, Eq. 22 is only valid
in the limit where $r_{\rm vir}\geq r_{\mathrm{hJ}}$. For a FDM particle mass
of $2.5\times 10^{-22}$ eV at $f=0.1$, this translates to a halo mass of about
$10^{6}\;M_{\odot}/h$. Halos for which this does not hold are taken to be
devoid of FDM in the halo model of Ref. [61].

In Fig. 9, we show a halo of simulation box 3 and its surrounding region. By
plotting the density isocontours, we notice that the FDM component is very
diffuse while the CDM has a steep density gradient. This is shown
quantitatively in the lower left panel of Fig. 9 where the CDM is more than
100 times denser in the center of the region than the FDM. When calculating
the Jeans radius for this halo, we find no value of $r_{\mathrm{hJ}}$
satisfying Eq. (22) and conclude that $r_{\rm vir}<r_{\mathrm{hJ}}$. The
findings shown in Fig. 9 supports the approximations of the halo model as the
halo is composed at $\approx 97.5$% of CDM.

In the limit of high halo mass with low FDM fraction, we find that the total
(combined) density is largely unaffected by the presence of FDM. This is shown
for fractions $f\leq 0.1$ in Fig. 10. We fit the radial density of the most
massive halo in each of the boxes 1-3 using a NFW profile and found the
concentration parameter to be unchanged by small amounts of FDM even if the
halo central density is slightly reduced. Moreover, the FDM fraction
calculated within the three halos of Fig. 10 are $(0.000,0.107,0.011)$ for
cosmological FDM fractions $f=(0.00,0.10,0.01)$, respectively. Thus, we find
that massive halos with a radius larger than the Jeans radius accumulate a
concentration of FDM equal to the cosmological FDM fraction. This matches the
prediction of the biased tracer model of Ref. [61].

Figure 10: (Top) Radial density profile of the most massive halos in low FDM
fraction simulation boxes for a FDM mass of $m=2.5\times 10^{-22}$ eV and a
redshift $z=4$. The dotted lines represent the NFW profile fits to the total
density. (Bottom) Relative difference in density with respect to the CDM-only
($f=0$) case.

###  III.3 Halo Mass Function

In this section, we measure the halo abundance and compare our results with
the theoretical predictions of a halo model for mixed dark matter. In a
mixture of CDM and FDM, the CDM is free to coalesce in small halos while the
FDM substructure is still washed out by wave effects. Therefore, unlike
cosmologies with a FDM fraction of $f=1$, halos with a radius below the FDM
Jeans scale will still form in mixed dark matter. To investigate this
quantitatively, we measure the halo mass function (HMF) in simulation boxes 1
and 3 (see Table 1).

We use the CDM particles as a tracer of the halos and use the friends-of-
friends (FOF) algorithm [62] to identify the halo centers. To correct for the
presence of FDM, we create a spherical profile around the halo center and vary
its radius until the total mean density of the halo reaches the virial
overdensity given in Eq. (20). We only use this approach in simulations where
CDM composes at least 90% of the dark matter and where FDM acts as a tracer of
the CDM. The CDM particles have a smaller separation between them in high
density regions, and not all halos that are resolved with particles have a
high enough number of grid cells to apply the FDM mass correction. About 80
halos had a radius sufficiently large to accurately account for the presence
of FDM. The final halo masses used are found by summing over the mass of the
particles in the halos and applying the FDM halo mass correction obtained from
the spherical profile. The other halos identified by the FOF algorithm had a
radius smaller or comparable to the FDM de Broglie wavelength. One such halo
which had a sufficiently large radius to encompass many FDM grid cells is
shown in Fig. 9. Most of the halos, however, did not have a sufficient number
of grid points to resolve the FDM density profile and the FDM mass correction
assumed a flat density profile equal to the mean FDM density.

We compare our results with the halo model predictions of axionHMCode [61],
which is a an adaptation of HMCode [63] accounting for ultralight particles
composing part of the dark matter. It is based on the halo model first
implemented in Ref. [22] to compute weak lensing shear statistics with FDM. We
show our results along with the model prediction in Fig. 11. We split the
simulation volume in four boxes and compute the error bars using a jackknife
method. To improve the accuracy of the comparison, we also rescale the HMF
prediction. This rescaling is necessary due to the fact that the HMF model
assumes a different mass definition than the output of the FOF halo finder. We
use the python package hmf [64, 65] to compute the ratio of HMFs using the FOF
and virial mass definitions (for a more in-depth discussion of halo mass
definitions see [66]). The rescaling is roughly mass-independent and has a
value of $\sim 0.83$. This difference is negligible compared to the error
margin on the HMF measurement and does not affect our results. We find a
reasonable agreement between the model and the simulations. We expect the halo
mass functions for different FDM fractions to converge at halo masses larger
than what can be captured by our box size as shown by the model prediction.
Given our box size, the results for the HMF are somewhat inconclusive,
although in the largest mass bins, the measured HMF shows a similar decrement
to the theoretical HMF (see the inset of Fig. 11).

Figure 11: Halo mass function (number density of halos above a certain mass
threshold) for two FDM fractions at redshift $z=4$. The colored regions are
measurements from simulations with the 68% confidence intervals and the solid
lines correspond to the predictions from axionHMCode. (Inset) Relative
difference in the number density of halos in the mass range
$10^{7}-10^{8}\;M_{\odot}/h$ compared to the model prediction (red dashed).

###  III.4 Halo Shapes

Halo triaxiality has been suggested as a way to detect dark matter properties
beyond CDM, such as dark matter self-interactions [67, 68, 69]. In this
section, we investigate the impact of FDM on the shape of dark matter halos.
Given the presence of a spherical solitonic core, we expect that the FDM will
lead to more spherical halo centers than pure CDM. Further away from the halo
center, a suppression of the initial power spectrum can also give more
spherical halos in the absence of wave effects as shown in Ref. [53]. To test
this prediction, we compare the shape of halos for FDM fractions
$f=\\{0.00,\;0.01,\;0.10,\;0.50\\}$ and use the central halo in simulation
boxes 1, 2, 3, and 4.

Halo shapes are parameterized by looking at the length of the halo’s principal
axes which are labelled in increasing size $a\geq b\geq c$. To compute the
direction and relative length of the principal axes, we diagonalize the
reduced inertia tensor $\mathcal{I}$ of the halo. The inertia tensor can be
computed using different methods [70], but we follow the approach of Refs [71,
67, 72] giving

$\displaystyle\mathcal{I}_{ij}=\sum_{n\;\rm
in\;halo}\left[1+\delta_{\mathrm{tot}}(\mathbf{x}_{n})\right]\frac{x_{n,i}x_{n,j}}{r_{n}^{2}},$
(23)

where $i,\;j=1,2,3$ are the indexes of the axes, $x_{n}$ is the distance from
the center of the halo to the $n^{\rm th}$ grid point in the ellipsoid axis
frame, $r_{n}=\sqrt{x_{n,1}^{2}+x_{n,2}^{2}/q^{2}+x_{n,3}^{2}/s^{2}}$ is the
ellipsoidal radius. The ellipsoidal radius is obtained from the axes ratios
$q=b/a$ and $s=c/a$. Normally, the inertia tensor is computed as a sum over
particles of the same mass. In the case at hand, we take a weighted sum over
grid points located within the central halo. We interpolate the particle mass
density of the CDM using a cloud-in-cell scheme on the same grid as the FDM
and we weigh the grid points by their total matter overdensity $\delta_{\rm
tot}=f\delta_{\rm FDM}+(1-f)\delta_{\rm CDM}$. We then compute this tensor
iteratively following the procedure of Ref. [72]. We begin by calculating the
principal axes for a halo centered at the point of highest density. The radius
of the sphere is taken to be the virial radius. Then, we create an ellipsoid
of the same volume as the initial sphere, with the principal axes aligned with
the eigenvectors of the inertia tensor and with the proper axes ratios. The
lengths of the axes simply correspond to the square root of the eigenvalues of
the tensor. This ellipsoid defines a new halo for which we repeat the same
procedure until we have reached ten iterations (after which the inferred shape
of the halo converges). For a numerical implementation of this procedure, we
make use of the public repositories inertia-tensors [73] and rotations [74].

To quantify the ellipticity of halos, we calculate the triaxiality parameter
[75, 72, 76]

$\displaystyle T\equiv\frac{a^{2}-b^{2}}{a^{2}-c^{2}}.$ (24)

Using this parameter, we can categorize halos as oblate ($T\leq 0.33$),
prolate ($T\geq 0.66$), or triaxial ($0.33\leq T\leq 0.66$).

Figure 12: Triaxial parameter of the central halo for various FDM fractions.
We note a more spherical inner halo region with high fraction due to the
presence of the spherical soliton core. The $f=1$ halo was too diffuse and the
ellipsoid algorithm did not converge.

In Fig. 12, we show the impact of having a high fraction at various radii. We
note that the effects on halo shape only manifest at around $f=0.5$, where the
central region within $0.5\;r_{\rm vir}$ is more spherical. Contrary to past
studies on halo triaxiality, we are not comparing halo populations but
identical halos in simulations with different dark matter contents. We also
only study halos at high redshift $z=4$. Since the formation of a soliton core
is universal in halos containing ultralight particles, we expect our
observations to hold for halo populations at lower redshifts. We note that our
algorithm did not converge for the FDM-only $f=1$ case as the structure was
too diffuse and contained a significant number of interference fringes.

##  IV Discussion

In this study, we have explored the non-linear behaviour of dark matter
composed of a mixture of cold and ultralight (fuzzy) dark matter. We adapted a
simulation algorithm for dark matter models with gravitationally coupled CDM
and FDM components evolved simultaneously.

We first examined the non-linear power spectrum of each component up to
redshift $z=4$. While the resulting spectra proved too noisy for a comparison
with halo model predictions, we confirmed that the presence of wave effects in
the non-linear regime leads to enhanced suppression of the power spectrum when
compared to simulations where only the initial conditions had been modified.
This agrees with the findings of AXIREPO and AX-GADGET codes [41, 15]. In
simulations with $f=1$ and $m\sim 10^{-22}$ eV, it was found that the wave
effects can cause an increase in power on scales approaching $k\sim
1000\;h/$Mpc [41] which we do not resolve sufficiently well with a $512^{3}$
grid (see Appendix B). In the low fraction limit, we found that the
suppression of the power spectrum was mostly attributable to the change in the
initial conditions rather than the wave effects. We also observed that the
excess suppression due to the wave effects scales roughly linearly with $f$ if
$f\lesssim 0.1$.

Next, we investigated the density profiles of halos with mixed dark matter. We
found that, for high fractions $f\geq 0.50$, the halo exhibited a cored
density profile matching closely the fitting formula of Ref. [40] even in the
presence of CDM. For a fraction of $f=0.5$, the CDM created a steep potential
well exerting a pressure on the FDM soliton and caused its core radius to
shrink. For lower fractions, we found that no stable core could form due to
the enhanced gravitational potential in the presence of CDM. We compared the
predicted size of the core (accounting for the CDM density) with the measured
density profile of the FDM and concluded that the ground-state solution would
have a radius large enough to be detected given our numerical resolution. Our
findings that the soliton core does not form at a fraction of $f=0.1$ are in
agreement with the results of spherical collapse simulations [13].

We ran a customized halo finder and measured the halo mass function for pure
CDM and 10% FDM with a mass of $2.5\times 10^{-22}$ eV. Our implementation of
the halo finder is a combination of the friends-of-friends algorithm to
identify the halo centers with CDM particles and the spherical overdensity
finder to correct for the FDM density on a grid. The final halo finder gave us
a halo catalog with mass measurements accounting for both dark matter
components. We could not be conclusive in our axionHMCode comparison, but
found a qualitative agreement given the limited statistics. We also identified
a low mass halo below the predicted cut-off introduced in the biased tracer
halo model of Ref. [61] and found that it had a very small FDM fraction, as
expected.

As is the case for the power spectra, testing non-linear halo models further
will benefit from cosmological simulations with a larger halo population and a
lower final redshift. The time-step requirements for the Schrödinger-Poisson
system make this task challenging and such simulations may require other
computational approaches to evolve the FDM dynamics. This is exacerbated in
the presence of a steep external gravitational potential created by the CDM
which shrinks the de Broglie wavelength of the FDM thus increasing the
resolution needed to study the halo cores. We leave more involved analyses of
halo statistics in mixed dark matter for future work.

Finally, we investigated the possibility that FDM could reduce the halo
ellipticity around the soliton cores since a similar effect has been shown to
occur for models of self-interacting dark matter. We expected the wave effects
and modified initial conditions to isotropize the halo density on scales
comparable to the FDM de Broglie wavelength and found this to be the case for
the high FDM fraction ($f=0.50$). However, we found no trace of this
phenomenon in halos with $f\leq 0.10$. It is possible that the impact of the
wave effects manifested on scales which escaped the simulation resolution. In
a mixed dark matter cosmology, we found that the FDM is not distributed
equally across halos of different masses as displayed in Fig. 1.

In this study, we have run the first cosmological simulations of a mixture of
cold and fuzzy dark matter. We have found that the resulting dark matter model
combined features from both its constituents but exhibited a unique behaviour
amongst known dark matter models. These findings will have a profound impact
in the study of ultralight particles and the modelling of their behaviour on
non-linear scales.

###### Acknowledgements.

AL acknowledges support from NASA grant 21-ATP21-0145. AL would like to thank
Jens Niemeyer, Benedikt Eggemeier, Mateja Gosenca, Tibor Dome, Simon May,
Mathew Madhavacheril, J. Richard Bond, and Daniel Grin for useful discussions.
RH acknowledge support from the Natural Sciences and Engineering Research
Council of Canada. The Dunlap Institute is funded through an endowment
established by the David Dunlap family and the University of Toronto. The
authors at the University of Toronto acknowledge that the land on which the
University of Toronto is built is the traditional territory of the
Haudenosaunee, and most recently, the territory of the Mississaugas of the New
Credit First Nation. They are grateful to have the opportunity to work in the
community, on this territory. R. H. additionally acknowledges support from
CIFAR, and the Azrieli and Alfred. P. Sloan Foundations. DJEM is supported by
an Ernest Rutherford Fellowship from the STFC, Grant No. ST/T004037/1 and by a
Leverhulme Trust Research Project (RPG-2022-145). Computations described in
this work were performed with resources provided by the North-German
Supercomputing Alliance (HLRN). The analysis was performed on the Niagara
supercomputer at the SciNet HPC Consortium. SciNet is funded by Innovation,
Science and Economic Development Canada; the Digital Research Alliance of
Canada; the Ontario Research Fund: Research Excellence; and the University of
Toronto [77, 78].

##  Appendix A Pade Approach to SP System

We propose an alternative method for solving the system of equations (16). The
system is a boundary value problem with conditions at $r\to\infty$ and is
therefore potentially unstable numerically. Moreover, it allows an infinite
number of solutions, only one of which (the zero node solution) is stable
[79]. The usual method is to first transform the system in dimensionless
variables

$\displaystyle r$ $\displaystyle\to\frac{m}{\hbar}r,$ (25) $\displaystyle\phi$
$\displaystyle\to\frac{\hbar\sqrt{4\pi G}}{m^{2}}\phi,$ (26)
$\displaystyle\rho_{\mathrm{CDM}}$ $\displaystyle\to\frac{4\pi
G\hbar^{2}}{m^{2}}\rho_{\mathrm{CDM}},$ (27) $\displaystyle\gamma$
$\displaystyle\to\frac{\hbar}{m}\gamma.$ (28)

Then, we can use the scaling of the system using a scalar $\lambda$ [57, 79]
$\\{r,\phi,\rho_{\mathrm{CDM}},\gamma\\}\to\\{\lambda^{-1}\tilde{r},\lambda^{2}\tilde{\phi},\lambda^{4}\tilde{\rho}_{\mathrm{CDM}},\lambda^{2}\tilde{\gamma}\\}$
where we choose $\lambda=\phi(r=0)^{1/2}$ so that $\tilde{\phi}(r=0)=1$. Our
new rescaled system of equations is then

$\displaystyle\frac{\partial^{2}(\tilde{r}\tilde{\phi})}{\partial\tilde{r}^{2}}$
$\displaystyle=2\tilde{r}\left(\tilde{V}-\tilde{\gamma}\right)\tilde{\phi},$
(29)
$\displaystyle\frac{\partial^{2}(\tilde{r}\tilde{V})}{\partial\tilde{r}^{2}}$
$\displaystyle=\tilde{r}\left(\tilde{\phi}^{2}+\tilde{\rho}_{\mathrm{CDM}}\right),$
(30)

with boundary conditions

$\displaystyle\tilde{\phi}(0)$ $\displaystyle=1,$ (31)
$\displaystyle\tilde{\phi}(\tilde{r}\to\infty)$ $\displaystyle=0$ (32)
$\displaystyle\tilde{\phi}^{\prime}(0)$ $\displaystyle=0,$ (33)
$\displaystyle\tilde{V}^{\prime}(0)$ $\displaystyle=0$ (34)
$\displaystyle\tilde{V}(r\to\infty)$ $\displaystyle=0,$ (35)

where the prime denotes differentiation with respect to $r$.

Some tricks have been suggested to approximate the problematic boundary
conditions at infinity [79], but here we consider a complementary approach to
this problem which satisfies the exact boundary conditions and closely
approximates $\phi$. Namely, we define the inverse polynomial estimator

$\displaystyle\hat{\phi}(x)=\left(1+\sum_{k=1}^{n}a_{k}x^{k}\right)^{-1}.$
(36)

This resembles the Padé approximant with the condition that the polynomial in
the numerator is set to unity. This ensures that the condition at $r=0$ is
met. We note that the boundary conditions at infinity are trivially satisfied
given that the inverse polynomial is monotonically decreasing. Given the
success of the soliton fit with $\phi\propto r^{-8}$ at large $r$, we pick
$n=8$.

From Ref. [57], we can solve for the potential through

$\displaystyle
V(r)=V_{0}+\int_{0}^{r}y\left(\phi^{2}+\rho_{\mathrm{CDM}}\right)dy-\frac{M(r)}{r},$
(37)

where $M$ is the number density

$\displaystyle
M(r)\equiv\int_{0}^{r}y^{2}\left(\phi^{2}+\rho_{\mathrm{CDM}}\right)dy.$ (38)

For simplicity, let us consider only the case where the CDM is absent. Then
the potential $\tilde{V}^{\prime}\propto r^{1-N}$ for large $r$ and is equal
to zero at $r=0$. The same goes for $\tilde{\phi}$. Finally, the boundary
condition $\tilde{V}(r\to\infty)=0$ can be satisfied by an appropriate choice
of $V_{0}$. The problem is now to find the coefficients $a_{k}$ so that
$\hat{\phi}\approx\tilde{\phi}$ without having an analytic representation for
$\tilde{\phi}$. For this, we borrow the definition of the loss function
implemented in physics-informed neural networks [80]. We define

$\displaystyle\mathcal{L}\equiv\frac{1}{N}\sum_{r_{i}\in
R}\left[\partial^{2}_{\tilde{r}}\left(\tilde{r}_{i}\hat{\phi}\right)\bigg{|}_{r=r_{i}}-2\tilde{r}_{i}\left(\tilde{V}_{i}-\tilde{\gamma}\right)\hat{\phi}_{i}\right]^{2},$
(39)

where $\hat{\phi}_{i}=\hat{\phi}(r_{i})$, $\tilde{V}_{i}=\tilde{V}(r_{i})$,
and where $R$ is a set of $N$ sample points between 0 and $r_{\rm max}$. We
note here that we can be quite flexible in our choice of $r_{\rm max}$ given
that our sample solution satisfies all the boundary conditions. To
approximately solve this system, we then minimize the loss function over our
coefficients $a_{k}$. We also have two other unknowns in $\tilde{V}_{0}$ and
$\tilde{\gamma}$. However, given the form of Eq. (29), it is immediate that
there is a degeneracy between the two parameters and that we can simply
combine them in one $\tilde{\gamma}_{0}\equiv\tilde{V}_{0}-\tilde{\gamma}$.
This gives a total $N+1$ free parameters over which to minimize. We employ the
basinhopping method implemented with the SciPy Python package [81] for this
procedure. Having found the $a_{k}$ coefficient, we compare our approximated
solution to the numerical result when using the shooting method. The two
solutions are plotted in Fig. 13 where we note the remarkable agreement
between the two. The advantage of the method we develop here is that it allows
us to solve the system for a variety of gravitational potential shapes without
having to venture a guess about the value of $\gamma$. In other instances,
this could be problematic as the wrong guess would lead to unstable solutions
with a non-zero number of nodes.

Figure 13: Comparison between the numerical solution using the shooting method
and the Padé approximant approach in the pure FDM case. We note an agreement
to about $0.5\%$.

##  Appendix B Impact of Box Length and Resolution

In Fig. 5, we observe a change in the wavenumber at which the matter power
spectrum deviates form the linear theory prediction as a function of FDM mass.
This is particularly visible for high masses $m\sim 10^{-21}$ eV. However, we
expect to recover the CDM non-linear power spectrum in the limit of
$m\to\infty$. The main difference in configuration between the simulations of
different FDM masses is the choice of box length as listed in Table 1. We
investigate if this change in box size is responsible for the varying non-
linear scales.

First, we define what we mean by the non-linear scale which we label $k_{\rm
NL}$. For this we use he approach of Ref. [82] and define

$\displaystyle\sigma^{2}_{d}(k,k_{\rm min})\equiv\frac{4\pi}{3}\int_{k_{\rm
min}}^{k}dq\;P(q),$ (40)

where $P$ is the linear matter power spectrum. From the above definition, we
consider that non-linearities arise when $k\sigma_{d}\gg 1$. In analytic
calculations, we omit the lower bound on the integral since we can generate
the linear matter power spectrum to very small wavenumbers. In our case, since
some of our simulation volumes are very small, the coupling between large and
small-scale modes is limited to the scales below the fundamental scale of the
box. To account for this, we solve for the non-linear scale using

$\displaystyle k_{\rm NL}\sigma_{d}(k_{\rm NL},2\pi/L_{\rm box})=5,$ (41)

where $L_{\rm box}$ is the size of the box. Note that the factor of five on
the RHS of Eq. (41) is arbitrary, but our conclusion are largely unchanged for
any number greater than unity but less than order 10. We display the
calculated $k_{\rm NL}$ as a function of FDM mass and box size in Fig. 14. All
the calculations involved in this appendix are done assuming a redshift $z=4$.

Figure 14: Non-linear and fundamental (minimum) scales for simulation boxes 3,
6, 7, 8, 9 (see Table 1) as a function of FDM mass and box size. The grey line
denotes the non-linear scale in the limit of an infinitely large simulation
volume, matching the theoretical calculation. We observe a strong dependence
of the non-linear scale on box size which matches the result of Fig. 5 for
$z=4$.

We also compare the density fields of boxes 3 and 10 which have the same
cosmology, but different resolutions. We take the Fourier transform of the
(total) density fields in both boxes at redshifts $z=15,\;7,\;4$ and compute
the correlation coefficient

$\displaystyle r(k)\equiv\frac{\langle\delta_{\rm LR}(\mathbf{k})\delta_{\rm
HR}(\mathbf{k})\rangle}{\sqrt{\langle\delta_{\rm
LR}^{2}(\mathbf{k})\rangle\langle\delta_{\rm HR}^{2}(\mathbf{k})\rangle}},$
(42)

where $\delta_{\rm LR,\;HR}(\mathbf{k})$ are the density fields of the low and
high-resolution boxes (respectively boxes 3 and 10). We plot the results of
this calculation in Fig. 15.

Figure 15: Cross-correlation coefficient of the $1024^{3}$ and $512^{3}$
resolution boxes as a function of redshift.

At early times, the system is still very linear and the difference in
resolution does not impact the density field. However, we note a degradation
of the correlation coefficient at scales above $\sim 100\;h/$Mpc which we
attribute to small-scale wave fluctuations around the mean density.

## References

  * Abdalla  _et al._ [2022] E. Abdalla _et al._ , Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies, Journal of High Energy Astrophysics 34, 49 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Hu  _et al._ [2000] W. Hu, R. Barkana, and A. Gruzinov, Fuzzy Cold Dark Matter: The Wave Properties of Ultralight Particles, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 [astro-ph] .
  * Rogers  _et al._ [2023] K. K. Rogers, R. Hložek, A. Laguë, M. M. Ivanov, O. H. E. Philcox, G. Cabass, K. Akitsu, and D. J. E. Marsh, Ultra-light axions and the S 8 tension: joint constraints from the cosmic microwave background and galaxy clustering, JCAP 06, 023, arXiv:2301.08361 [astro-ph.CO] .
  * Kim and Marsh [2016] J. E. Kim and D. J. E. Marsh, An ultralight pseudoscalar boson, Phys. Rev. D 93, 025027 (2016), arXiv:1510.01701 [hep-ph] .
  * Davoudiasl and Murphy [2017] H. Davoudiasl and C. W. Murphy, Fuzzy Dark Matter from Infrared Confining Dynamics, Phys. Rev. Lett. 118, 141801 (2017), arXiv:1701.01136 [hep-ph] .
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String axiverse, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Mehta  _et al._ [2021] V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. McAllister, and M. J. Stott, Superradiance in string theory, Journal of Cosmology and Astroparticle Physics 2021, 033 (2021), arXiv:2103.06812 [hep-th] .
  * Cicoli  _et al._ [2022] M. Cicoli, V. Guidetti, N. Righi, and A. Westphal, Fuzzy Dark Matter candidates from string theory, Journal of High Energy Physics 2022, 107 (2022), arXiv:2110.02964 [hep-th] .
  * Gendler  _et al._ [2023] N. Gendler, D. J. E. Marsh, L. McAllister, and J. Moritz, Glimmers from the Axiverse, arXiv e-prints , arXiv:2309.13145 (2023), arXiv:2309.13145 [hep-th] .
  * Alexander  _et al._ [2023] S. Alexander, H. Gilmer, T. Manton, and E. McDonough, The $\pi$-axion and $\pi$-axiverse of dark QCD, arXiv e-prints , arXiv:2304.11176 (2023), arXiv:2304.11176 [hep-ph] .
  * Bachlechner  _et al._ [2018] T. C. Bachlechner, K. Eckerle, O. Janssen, and M. Kleban, Multiple-axion framework, Phys. Rev. D 98, 061301 (2018), arXiv:1703.00453 [hep-th] .
  * Mocz  _et al._ [2018] P. Mocz, L. Lancaster, A. Fialkov, F. Becerra, and P.-H. Chavanis, Schrödinger-Poisson-Vlasov-Poisson correspondence, Phys. Rev. D 97, 083519 (2018), arXiv:1801.03507 [astro-ph.CO] .
  * Schwabe  _et al._ [2020] B. Schwabe, M. Gosenca, C. Behrens, J. C. Niemeyer, and R. Easther, Simulating mixed fuzzy and cold dark matter, Phys. Rev. D 102, 083518 (2020), arXiv:2007.08256 [astro-ph.CO] .
  * Gosenca  _et al._ [2023] M. Gosenca, A. Eberhardt, Y. Wang, B. Eggemeier, E. Kendall, J. L. Zagorac, and R. Easther, Multifield ultralight dark matter, Phys. Rev. D 107, 083014 (2023), arXiv:2301.07114 [astro-ph.CO] .
  * Nori and Baldi [2018] M. Nori and M. Baldi, AX-GADGET: a new code for cosmological simulations of Fuzzy Dark Matter and Axion models, Monthly Notices of the Royal Astronomical Society 478, 3935 (2018), arXiv:1801.08144 [astro-ph.CO] .
  * Mocz  _et al._ [2019] 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, First Star-Forming Structures in Fuzzy Cosmic Filaments, Phys. Rev. Lett. 123, 141301 (2019), arXiv:1910.01653 [astro-ph.GA] .
  * Veltmaat  _et al._ [2018] J. Veltmaat, J. C. Niemeyer, and B. Schwabe, Formation and structure of ultralight bosonic dark matter halos, Phys. Rev. D 98, 043509 (2018), arXiv:1804.09647 [astro-ph.CO] .
  * Hopkins [2019] P. F. Hopkins, A stable finite-volume method for scalar field dark matter, Monthly Notices of the Royal Astronomical Society 489, 2367 (2019), arXiv:1811.05583 [astro-ph.CO] .
  * Hlozek  _et al._ [2015] R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, A search for ultralight axions using precision cosmological data, Phys. Rev. D 91, 103512 (2015), arXiv:1410.2896 [astro-ph.CO] .
  * Hložek  _et al._ [2018] R. Hložek, D. J. E. Marsh, and D. Grin, Using the full power of the cosmic microwave background to probe axion dark matter, Monthly Notices of the Royal Astronomical Society 476, 3063 (2018), arXiv:1708.05681 [astro-ph.CO] .
  * Laguë  _et al._ [2022] A. Laguë, J. R. Bond, R. Hložek, K. K. Rogers, D. J. E. Marsh, and D. Grin, Constraining ultralight axions with galaxy surveys, Journal of Cosmology and Astroparticle Physics 2022 (1), 049, arXiv:2104.07802 [astro-ph.CO] .
  * Dentler  _et al._ [2022] M. Dentler, D. J. E. Marsh, R. Hložek, A. Laguë, K. K. Rogers, and D. Grin, Fuzzy dark matter and the Dark Energy Survey Year 1 data, Monthly Notices of the Royal Astronomical Society 515, 5646 (2022), arXiv:2111.01199 [astro-ph.CO] .
  * Iršič  _et al._ [2017] V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, First Constraints on Fuzzy Dark Matter from Lyman-$\alpha$ Forest Data and Hydrodynamical Simulations, Phys. Rev. Lett. 119, 031302 (2017), arXiv:1703.04683 [astro-ph.CO] .
  * Armengaud  _et al._ [2017] E. Armengaud, N. Palanque-Delabrouille, C. Yèche, D. J. E. Marsh, and J. Baur, Constraining the mass of light bosonic dark matter using SDSS Lyman-$\alpha$ forest, Monthly Notices of the Royal Astronomical Society 471, 4606 (2017), arXiv:1703.09126 [astro-ph.CO] .
  * Rogers and Peiris [2021] K. K. Rogers and H. V. Peiris, Strong Bound on Canonical Ultralight Axion Dark Matter from the Lyman-Alpha Forest, Phys. Rev. Lett. 126, 071302 (2021), arXiv:2007.12705 [astro-ph.CO] .
  * Zhang  _et al._ [2018] J. Zhang, J.-L. Kuo, H. Liu, Y.-L. Sming Tsai, K. Cheung, and M.-C. Chu, The Importance of Quantum Pressure of Fuzzy Dark Matter on Ly$\alpha$ Forest, Astrophys. J. 863, 73 (2018), arXiv:1708.04389 [astro-ph.CO] .
  * Nori  _et al._ [2019] M. Nori, R. Murgia, V. Iršič, M. Baldi, and M. Viel, Lyman $\alpha$ forest and non-linear structure characterization in Fuzzy Dark Matter cosmologies, Monthly Notices of the Royal Astronomical Society 482, 3227 (2019), arXiv:1809.09619 [astro-ph.CO] .
  * Calabrese and Spergel [2016] E. Calabrese and D. N. Spergel, Ultra-light dark matter in ultra-faint dwarf galaxies, Monthly Notices of the Royal Astronomical Society 460, 4397 (2016), arXiv:1603.07321 [astro-ph.CO] .
  * Safarzadeh and Spergel [2020] M. Safarzadeh and D. N. Spergel, Ultra-light Dark Matter Is Incompatible with the Milky Way’s Dwarf Satellites, Astrophys. J. 893, 21 (2020), arXiv:1906.11848 [astro-ph.CO] .
  * Schutz [2020] K. Schutz, Subhalo mass function and ultralight bosonic dark matter, Physical Review D 101, 123026 (2020), arXiv:2001.05503 [astro-ph.CO] .
  * Dalal  _et al._ [2020] N. Dalal, J. Bovy, L. Hui, and X. Li, Don’t cross the streams: caustics from Fuzzy Dark Matter, arXiv e-prints , arXiv:2011.13141 (2020), arXiv:2011.13141 [astro-ph.CO] .
  * Arvanitaki and Dubovsky [2011] A. Arvanitaki and S. Dubovsky, Exploring the string axiverse with precision black hole physics, Physical Review D 83, 044026 (2011), arXiv:1004.3558 [hep-th] .
  * Stott and Marsh [2018] M. J. Stott and D. J. E. Marsh, Black hole spin constraints on the mass spectrum and number of axionlike fields, Physical Review D 98, 083006 (2018), arXiv:1805.02016 [hep-ph] .
  * Marsh [2016] D. J. Marsh, Axion cosmology (2016).
  * Hui  _et al._ [2017] L. Hui _et al._ , Ultralight scalars as cosmological dark matter, Physical Review D 95, 10.1103/PhysRevD.95.043541 (2017).
  * Widrow and Kaiser [1993] L. M. Widrow and N. Kaiser, Using the Schrodinger equation to simulate collisionless matter, Astrophys. J. 416, L71 (1993).
  * Almgren  _et al._ [2013] A. S. Almgren, J. B. Bell, M. J. Lijewski, Z. Lukić, and E. Van Andel, Nyx: A Massively Parallel AMR Code for Computational Cosmology, Astrophys. J. 765, 39 (2013), arXiv:1301.4498 [astro-ph.IM] .
  * Mocz  _et al._ [2017] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, Galaxy formation with BECDM - I. Turbulence and relaxation of idealized haloes, Monthly Notices of the Royal Astronomical Society 471, 4559 (2017), arXiv:1705.05845 [astro-ph.CO] .
  * Edwards  _et al._ [2018] F. Edwards, E. Kendall, S. Hotchkiss, and R. Easther, PyUltraLight: A Pseudo-Spectral Solver for Ultralight Dark Matter Dynamics, JCAP 10, 027, arXiv:1807.04037 [astro-ph.CO] .
  * Schive  _et al._ [2014] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Cosmic structure as the quantum interference of a coherent dark wave, Nature Physics 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * May and Springel [2023] S. May and V. Springel, The halo mass function and filaments in full cosmological simulations with fuzzy dark matter, Monthly Notices of the Royal Astronomical Society 524, 4256 (2023), arXiv:2209.14886 [astro-ph.CO] .
  * Li  _et al._ [2019] X. Li, L. Hui, and G. L. Bryan, Numerical and perturbative computations of the fuzzy dark matter model, Phys. Rev. D 99, 063509 (2019), arXiv:1810.01915 [astro-ph.CO] .
  * Schwabe  _et al._ [2016] B. Schwabe, J. C. Niemeyer, and J. F. Engels, Simulations of solitonic core mergers in ultralight axion dark matter cosmologies, Phys. Rev. D 94, 043513 (2016), arXiv:1606.05151 [astro-ph.CO] .
  * Mina  _et al._ [2020] M. Mina, D. F. Mota, and H. A. Winther, SCALAR: an AMR code to simulate axion-like dark matter models, Astronomy and Astrophysics 641, A107 (2020), arXiv:1906.12160 [physics.comp-ph] .
  * Glennon  _et al._ [2023] N. Glennon, N. Musoke, and C. Prescod-Weinstein, Simulations of multifield ultralight axionlike dark matter, Phys. Rev. D 107, 063520 (2023), arXiv:2302.04302 [astro-ph.CO] .
  * Nhan Luu  _et al._ [2023] H. Nhan Luu, P. Mocz, M. Vogelsberger, S. May, J. Borrow, S. H. H. Tye, and T. Broadhurst, Nested solitons in two-field fuzzy dark matter, arXiv e-prints , arXiv:2309.05694 (2023), arXiv:2309.05694 [astro-ph.CO] .
  * Zagorac  _et al._ [2023] J. L. Zagorac, E. Kendall, N. Padmanabhan, and R. Easther, Soliton formation and the core-halo mass relation: An eigenstate perspective, Phys. Rev. D 107, 083513 (2023), arXiv:2212.09349 [astro-ph.CO] .
  * Crocce  _et al._ [2006] M. Crocce, S. Pueblas, and R. Scoccimarro, Transients from initial conditions in cosmological simulations, Monthly Notices of the Royal Astronomical Society 373, 369 (2006), arXiv:astro-ph/0606505 [astro-ph] .
  * Laguë  _et al._ [2021] A. Laguë, J. R. Bond, R. Hložek, D. J. E. Marsh, and L. Söding, Evolving ultralight scalars into non-linearity with lagrangian perturbation theory, Monthly Notices of the Royal Astronomical Society 10.1093/mnras/stab601 (2021), arXiv:2004.08482 [astro-ph.CO] .
  * Hwang and Noh [2009] J.-C. Hwang and H. Noh, Axion as a cold dark matter candidate, Physics Letters B 680, 1 (2009), arXiv:0902.4738 [astro-ph.CO] .
  * Park  _et al._ [2012] C.-G. Park, J.-c. Hwang, and H. Noh, Axion as a cold dark matter candidate: Low-mass case, Phys. Rev. D 86, 083535 (2012), arXiv:1207.3124 [astro-ph.CO] .
  * Hahn and Abel [2011] O. Hahn and T. Abel, Multi-scale initial conditions for cosmological simulations, Monthly Notices of the Royal Astronomical Society 415, 2101 (2011), arXiv:1103.6031 [astro-ph.CO] .
  * Dome  _et al._ [2023] T. Dome, A. Fialkov, P. Mocz, B. M. Schäfer, M. Boylan-Kolchin, and M. Vogelsberger, On the cosmic web elongation in fuzzy dark matter cosmologies: Effects on density profiles, shapes, and alignments of haloes, Monthly Notices of the Royal Astronomical Society 519, 4183 (2023), arXiv:2208.03827 [astro-ph.CO] .
  * Turk  _et al._ [2011] M. J. Turk, B. D. Smith, J. S. Oishi, S. Skory, S. W. Skillman, T. Abel, and M. L. Norman, yt: A Multi-code Analysis Toolkit for Astrophysical Simulation Data, The Astrophysical Journal Supplement Series 192, 9 (2011), arXiv:1011.3514 [astro-ph.IM] .
  * Hand  _et al._ [2018] N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Seljak, and Z. Slepian, nbodykit: An Open-source, Massively Parallel Toolkit for Large-scale Structure, The Astronomical Journal 156, 160 (2018), arXiv:1712.05834 [astro-ph.IM] .
  * May and Springel [2021] S. May and V. Springel, Structure formation in large-volume cosmological simulations of fuzzy dark matter: impact of the non-linear dynamics, Monthly Notices of the Royal Astronomical Society 506, 2603 (2021), arXiv:2101.01828 [astro-ph.CO] .
  * Guzmán and Ureña-López [2004] F. S. Guzmán and L. A. Ureña-López, Evolution of the Schrödinger-Newton system for a self-gravitating scalar field, Phys. Rev. D 69, 124033 (2004), arXiv:gr-qc/0404014 [gr-qc] .
  * Veltmaat  _et al._ [2020] J. Veltmaat, B. Schwabe, and J. C. Niemeyer, Baryon-driven growth of solitonic cores in fuzzy dark matter halos, Phys. Rev. D 101, 083518 (2020), arXiv:1911.09614 [astro-ph.CO] .
  * Navarro  _et al._ [1996] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Structure of Cold Dark Matter Halos, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph] .
  * Bryan and Norman [1998] G. L. Bryan and M. L. Norman, Statistical Properties of X-Ray Clusters: Analytic and Numerical Comparisons, Astrophys. J. 495, 80 (1998), arXiv:astro-ph/9710107 [astro-ph] .
  * Vogt  _et al._ [2023] S. M. L. Vogt, D. J. E. Marsh, and A. Laguë, Improved mixed dark matter halo model for ultralight axions, Phys. Rev. D 107, 063526 (2023), arXiv:2209.13445 [astro-ph.CO] .
  * Davis  _et al._ [1985] M. Davis, G. Efstathiou, C. S. Frenk, and S. D. M. White, The evolution of large-scale structure in a universe dominated by cold dark matter, Astrophys. J. 292, 371 (1985).
  * Mead  _et al._ [2021] A. J. Mead, S. Brieden, T. Tröster, and C. Heymans, HMCODE-2020: improved modelling of non-linear cosmological power spectra with baryonic feedback, Monthly Notices of the Royal Astronomical Society 502, 1401 (2021), arXiv:2009.01858 [astro-ph.CO] .
  * Murray  _et al._ [2013] S. G. Murray, C. Power, and A. S. G. Robotham, HMFcalc: An online tool for calculating dark matter halo mass functions, Astronomy and Computing 3, 23 (2013), arXiv:1306.6721 [astro-ph.CO] .
  * Murray  _et al._ [2021] S. G. Murray, B. Diemer, Z. Chen, A. G. Neuhold, M. A. Schnapp, T. Peruzzi, D. Blevins, and T. Engelman, THEHALOMOD: An online calculator for the halo model, Astronomy and Computing 36, 100487 (2021), arXiv:2009.14066 [astro-ph.CO] .
  * More  _et al._ [2011] S. More, A. V. Kravtsov, N. Dalal, and S. Gottlöber, The Overdensity and Masses of the Friends-of-friends Halos and Universality of Halo Mass Function, The Astrophysical Journal Supplement Series 195, 4 (2011), arXiv:1103.0005 [astro-ph.CO] .
  * Davé  _et al._ [2001] R. Davé, D. N. Spergel, P. J. Steinhardt, and B. D. Wandelt, Halo Properties in Cosmological Simulations of Self-interacting Cold Dark Matter, Astrophys. J. 547, 574 (2001), arXiv:astro-ph/0006218 [astro-ph] .
  * Peter  _et al._ [2013] A. H. G. Peter, M. Rocha, J. S. Bullock, and M. Kaplinghat, Cosmological simulations with self-interacting dark matter - II. Halo shapes versus observations, Monthly Notices of the Royal Astronomical Society 430, 105 (2013), arXiv:1208.3026 [astro-ph.CO] .
  * Johana Gonzalez  _et al._ [2023] E. Johana Gonzalez, A. Rodriguez-Medrano, L. A. Pereyra, and D. García Lambas, Cluster halo shapes in CDM and SIDM models: Unveiling the DM particle nature using a weak lensing approach, arXiv e-prints , arXiv:2308.13461 (2023), arXiv:2308.13461 [astro-ph.CO] .
  * Zemp  _et al._ [2011] M. Zemp, O. Y. Gnedin, N. Y. Gnedin, and A. V. Kravtsov, On Determining the Shape of Matter Distributions, The Astrophysical Journal Supplement Series 197, 30 (2011), arXiv:1107.5582 [astro-ph.IM] .
  * Dubinski and Carlberg [1991] J. Dubinski and R. G. Carlberg, The Structure of Cold Dark Matter Halos, Astrophys. J. 378, 496 (1991).
  * Allgood  _et al._ [2006] B. Allgood _et al._ , The shape of dark matter haloes: dependence on mass, redshift, radius and formation, Monthly Notices of the Royal Astronomical Society 367, 1781 (2006), arXiv:astro-ph/0508497 [astro-ph] .
  * Campbell [2019] D. Campbell, Inertia tensors, https://github.com/duncandc/inertia_tensors (2019).
  * Campbell and Lanusse [2020] D. Campbell and F. Lanusse, Rotations, https://github.com/duncandc/rotations (2020).
  * Franx  _et al._ [1991] M. Franx, G. Illingworth, and T. de Zeeuw, The Ordered Nature of Elliptical Galaxies: Implications for Their Intrinsic Angular Momenta and Shapes, Astrophys. J. 383, 112 (1991).
  * Schneider  _et al._ [2012] M. D. Schneider, C. S. Frenk, and S. Cole, The shapes and alignments of dark matter halos, Journal of Cosmology and Astroparticle Physics 2012, 030 (2012), arXiv:1111.5616 [astro-ph.CO] .
  * Loken  _et al._ [2010] C. Loken _et al._ , SciNet: Lessons Learned from Building a Power-efficient Top-20 System and Data Centre, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 256 (2010) p. 012026.
  * Ponce  _et al._ [2019] M. Ponce _et al._ , Deploying a Top-100 Supercomputer for Large Parallel Workloads: the Niagara Supercomputer, arXiv e-prints , arXiv:1907.13600 (2019), arXiv:1907.13600 [cs.DC] .
  * Marsh and Pop [2015] D. J. E. Marsh and A.-R. Pop, Axion dark matter, solitons and the cusp-core problem, Monthly Notices of the Royal Astronomical Society 451, 2479 (2015), arXiv:1502.03456 [astro-ph.CO] .
  * Raissi  _et al._ [2017] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations, arXiv e-prints , arXiv:1711.10561 (2017), arXiv:1711.10561 [cs.AI] .
  * Virtanen and SciPy 1.0 Contributors [2020] P. Virtanen and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Blas  _et al._ [2013] D. Blas, M. Garny, and T. Konstandin, On the non-linear scale of cosmological perturbation theory, Journal of Cosmology and Astroparticle Physics 2013, 024 (2013), arXiv:1304.1546 [astro-ph.CO] .