Recompute the warm cutoff, the fuzzy cutoff and the soliton

Read How warm and fuzzy dark matter erase small ripples and Three ways to make a dark matter core first. You need Python, NumPy and SciPy, and you should be comfortable reading an ordinary differential equation. Everything here is at redshift zero and in linear theory unless a step says otherwise. The next chapter, Reading the Milky Way satellite limits on warm and fuzzy dark matter, puts the half-mode mass to work on a published satellite census.

What is in the package

All files sit in one folder on this site. The frozen reference grid holds the cold dark matter power spectrum and the warm dark matter transfer functions computed with the Boltzmann code CLASS version 3.4.2. The script that generated it is the only file that needs CLASS. Two benchmark scripts read the grid or need no input at all: the linear cutoff benchmark and the isolated-halo benchmark. Their saved outputs are linear-output.json and isolated-output.json. A third script, compute_teaching_values.py, recomputes every worked-example number in chapters 2 to 4 and writes teaching-values.json.

The cosmology is the Planck 2018 combination TT,TE,EE+lowE+lensing+BAO, last column of Table 2 in Planck Collaboration 2020, arXiv:1807.06209v4: $h=0.6766$, $\omega_{\mathrm b}=0.02242$, $\omega_{\mathrm c}=0.11933$, $A_{\mathrm s}=2.105\times10^{-9}$, $n_{\mathrm s}=0.9665$, $\tau=0.0561$. Neutrinos are massless, a stated simplification. The grid has 141 wavenumbers spaced evenly in logarithm from $0.05$ to $300\,h\,\mathrm{Mpc^{-1}}$. Units are $h\,\mathrm{Mpc^{-1}}$ for $k$ and $(h^{-1}\,\mathrm{Mpc})^3$ for $P$. The warm runs give all the dark matter to a thermal relic of 2, 3 or 6.5 keV, with its temperature tuned so that its density equals the cold value.

The versions used for the saved outputs were Python 3.9.6, NumPy 2.0.2 and SciPy 1.13.1. We have run the scripts with these versions only. Check that you have the same grid first. Its SHA-256 digest begins with 04e1c5d5 and ends with 3bd48e.

shasum -a 256 reference-grid.json
python3 -I benchmark_linear_cutoffs.py reference-grid.json linear-output.json
python3 -I benchmark_isolated_soliton.py isolated-output.json
python3 -I compute_teaching_values.py reference-grid.json linear-output.json teaching-values.json

The flag that follows python3 runs each script in isolated mode, so a stray file in the current folder cannot be imported in place of NumPy. Regenerating the grid itself needs the classy wrapper for CLASS.

python3 -I generate_frozen_reference_grid.py reference-grid.json

Check one: the warm half-mode against CLASS

The Viel fit gives the amplitude ratio $T(k)=[1+(\alpha k)^{2\nu}]^{-5/\nu}$ with $\nu=1.12$ and $\alpha$ from the particle mass. Viel et al. 2005, equations 6 and 7 The half-mode has a closed form, $k_{\mathrm{hm}}=(2^{\nu/5}-1)^{1/(2\nu)}/\alpha$. From CLASS we take the stored ratio $T=\sqrt{P_{\mathrm{WDM}}/P_{\mathrm{CDM}}}$ and find where it first drops below $1/2$, interpolating linearly in $\log k$ and $\log T$ between grid points. The half-mode mass then follows from

Nadler et al. 2021, equation 2 with $\Omega_{\mathrm m}=0.3096$ and $\rho_{\mathrm{crit}}=2.775\times10^{11}\,h^2\,\mathrm M_\odot\,\mathrm{Mpc^{-3}}$.

Wavenumbers are in $h\,\mathrm{Mpc^{-1}}$ and masses in solar masses. Because the mass goes as $k_{\mathrm{hm}}^{-3}$, a 12 per cent shift in wavenumber becomes a 45 per cent shift in mass at 6.5 keV. The gap grows with particle mass, and the CLASS calculation always cuts at smaller scales than the fit. Viel and collaborators chose $\nu=1.12$ as the best fit to Boltzmann-code output for $k<5\,h\,\mathrm{Mpc^{-1}}$. These half-modes lie at 20 to 86, outside that range, so per cent agreement was never promised.

A second comparison is with a number someone else published. Nadler and collaborators write $M_{\mathrm{hm}}=5\times10^8\,\mathrm M_\odot\,(m_{\mathrm{WDM}}/3\,\mathrm{keV})^{-10/3}$. Nadler et al. 2021, equation 5 Their simulations use $h=0.7$ and $\Omega_{\mathrm m}=0.286$. Putting that cosmology and a warm density of $0.24$ into the Viel fit gives $4.81\times10^8\,\mathrm M_\odot$ at 3 keV, 4 per cent below their value. The warm density of $0.24$ is our own assumption for their setup.

Check two: the fuzzy cutoff against a direct growth calculation

Hu, Barkana and Gruzinov fitted the fuzzy amplitude ratio as $T=\cos x^3/(1+x^8)$ with $x=1.61\,m_{22}^{1/18}\,k/k_{\mathrm{J,eq}}$ and $k_{\mathrm{J,eq}}=9\,m_{22}^{1/2}\,\mathrm{Mpc^{-1}}$, and stated that the power falls to one half near $k_{1/2}\approx4.5\,m_{22}^{4/9}\,\mathrm{Mpc^{-1}}$. Hu, Barkana and Gruzinov 2000, equations 8 and 9 Here $m_{22}$ is the particle mass in units of $10^{-22}\,\mathrm{eV}$, and these wavenumbers are in $\mathrm{Mpc^{-1}}$ without the factor $h$.

The independent check integrates the linear growth of one wave. In the fluid form of the wave equation the particle's wavelength acts as a pressure. With $\delta$ the density contrast and $\ln a$ as time,

$\Omega_{\mathrm{cl}}(a)$ is the fraction of the critical density in matter that clusters at that time. The script lets only dark matter cluster before recombination at $a=1/1100$ and all matter afterwards, on a background of radiation, matter and a cosmological constant. Each wave starts at horizon entry with $\delta=1$ and $d\delta/d\ln a=1$, and the same wave is run again without the quantum term. The ratio of the two final values is the amplitude ratio, and the script finds where its square crosses $1/2$.

The growth equation cuts power at 11 to 17 per cent smaller wavenumber than the fit. This calculation is a simplified model: it has no baryon pressure, no photon coupling and a sharp switch at recombination. So the table says the fit and this model differ by that much. It does not say which is closer to a full Boltzmann code for fuzzy dark matter. The next section shows that the result also depends on the starting slope, which is part of why.

Check three: the soliton as an exact ground state

The core of a fuzzy halo is the lowest-energy, self-gravitating state of the wave. In units where $\hbar/m=1$ and $G=1$ it satisfies

where $\rho=\psi^2$ and $u$ is the gravitational potential minus the energy. The script fixes $\psi(0)=1$ and adjusts $u(0)$ by bisection until $\psi$ falls to zero without crossing it. The equations keep their form under $r\to r/\lambda$, $\rho\to\lambda^4\rho$, so $\rho_0r_{\mathrm c}^4$ is the same number for every ground state, with $r_{\mathrm c}$ the radius where the density falls to half.

Schive and collaborators fitted their simulated cores with $\rho_0=1.9\,(m/10^{-23}\,\mathrm{eV})^{-2}(r_{\mathrm c}/\mathrm{kpc})^{-4}\,\mathrm M_\odot\,\mathrm{pc^{-3}}$ and the shape $[1+0.091(r/r_{\mathrm c})^2]^{-8}$. Schive et al. 2014, PRL, equation 3 Converted to these units, their prefactor is $0.2223$.

Both tolerances, $10^{-7}$ and $10^{-10}$, give the same numbers to six figures. A shape that halves exactly at $r_{\mathrm c}$ needs the coefficient $2^{1/8}-1=0.0905$, so $0.091$ is that value rounded, and the central part matches well. The error grows outward: $0.4$ per cent at $r_{\mathrm c}$, $4.0$ per cent at $2r_{\mathrm c}$, and the largest, $5.9$ per cent, at $2.7r_{\mathrm c}$. The Letter's sentence about 2 per cent does not say what the fit was compared with. Our 5.9 per cent is against the exact ground state, and we have not tested the fit against their simulated cores.

The mass inside the core radius times that radius is fixed by the same scaling. For $m=10^{-23}\,\mathrm{eV}$ the exact state gives $M_{\mathrm c}r_{\mathrm c}=5.41\times10^9\,\mathrm M_\odot\,\mathrm{kpc}$. Integrating the Schive formula from $0$ to $r_{\mathrm c}$ gives $5.30\times10^9$, 2 per cent lower, about the same as the gap in the prefactor.

The same script checks the isothermal sphere used in the self-interacting model. Its numerical solution agrees with the series $x^2/6-x^4/120+x^6/1890$ to $2.4\times10^{-8}$ for $x\le0.4$. Far out, the density times $x^2/2$ should approach one, and it reads $1.016$, $1.049$ and $0.997$ at $x=200$, $600$ and $2000$. That slow, swinging approach is expected for the isothermal sphere. The density halves at $x=2.2746$.

How each check can fail

Every one of these failures is saved in linear-output.json, so you can confirm you reproduce it.

The fuzzy fit changes sign. The factor $\cos x^3$ turns negative above $x=(\pi/2)^{1/3}$, which is $k=6.50\,\mathrm{Mpc^{-1}}$ at $m_{22}=1$. The fit is a cutoff shape, and only $T^2$ should be used beyond that point. Code that takes $\log T$ there will stop with an error.

A fixed-step integrator blows up. At $k=30\,\mathrm{Mpc^{-1}}$ and $m_{22}=1$ the quantum term is huge at early times, and the equation becomes stiff. An adaptive solver returns an amplitude ratio of $3.69\times10^{-4}$. A fourth-order Runge–Kutta solver with 40 fixed steps returns $-1.86\times10^{26}$. A result that large or negative means the step is too long, not that the physics changed.

The starting condition matters. With the default start, the half-power point at $m_{22}=1$ is $3.80\,\mathrm{Mpc^{-1}}$. Starting with zero initial slope gives $3.22$. Starting ten times earlier in scale factor gives $3.31$. A 15 per cent spread from the start alone is one more reason not to call Table 5.2 a precision result.

The warm fit fails deep in the tail. At the first grid point where CLASS has fallen to $T\approx0.017$, the Viel fit is already at $0.0027$ for 2 keV, $0.0019$ for 3 keV and $0.0011$ for 6.5 keV, six to seventeen times lower. Even while $T$ is above $0.25$, the fit is off by up to 16, 23 and 31 per cent for the three masses. Read CLASS directly when the tail matters.

A checklist for your own run

1. Confirm the grid digest before you start.
2. Run the three scripts and compare your output files with the saved ones. Agreement to four significant figures is the target.
3. Change one input at a time. Try a 4 keV particle in the generator or $m_{22}=2$ in the growth benchmark, and predict the direction of the change before running it.
4. Rerun with a tighter tolerance and report the number of digits that did not move.
5. Write down what each comparison does not test. A fit can agree with a calculation that shares its simplifications.

Lab tasks

1. Predict. Set the halo Lab to warm dark matter at 3 keV and read the half-mode mass. Is it the fit value or the CLASS value from Table 5.1?
2. Change one condition. Move the fuzzy particle mass from $m_{22}=1$ to $m_{22}=4$ and read the soliton core radius. At fixed halo mass it should shrink by a factor of 4, because $r_{\mathrm c}\propto m^{-1}$.
3. Explain. Using Table 5.3, explain why the Lab can draw the soliton with the Schive formula without visibly misplacing the core.

Exercises

1. Show from the closed form that the warm half-mode wavenumber scales as $m_{\mathrm{WDM}}^{1.11}$ at fixed cosmology, and that the half-mode mass scales as $m_{\mathrm{WDM}}^{-3.33}$.
2. Using Table 5.2, fit a power law to the growth-equation half-power points. Does it keep the exponent $4/9$ of Hu equation 9?
3. The 40-step solver fails at $k=30\,\mathrm{Mpc^{-1}}$. Estimate the early-time value of the quantum term relative to the gravity term, and explain why a fixed step in $\ln a$ is too long there.