Early-Universe Simulations of the Cosmological Axion

TL;DR

This paper presents high-resolution simulations of early-universe axion models, revealing the formation of ultracompact dark matter minihalos and providing insights into their primordial spectrum and implications for axion detection.

Contribution

It introduces detailed simulations of axion-induced minihalo formation from the PQ transition to matter-radiation equality, highlighting the generated perturbation spectrum and refining axion mass estimates.

Findings

Primordial perturbation spectrum characterized.

DM density consistent with axion mass around 25.2 μeV.

Implications for axion dark matter detection strategies.

Abstract

Ultracompact dark matter (DM) minihalos at masses at and below $10^{-12}$ $M_\odot$ arise in axion DM models where the Peccei-Quinn (PQ) symmetry is broken after inflation. The minihalos arise from density perturbations that are generated from the non-trivial axion self interactions during and shortly after the collapse of the axion-string and domain-wall network. We perform high-resolution simulations of this scenario starting at the epoch before the PQ phase transition and ending at matter-radiation equality. We characterize the spectrum of primordial perturbations that are generated and comment on implications for efforts to detect axion DM. We also measure the DM density at different simulated masses and argue that the correct DM density is obtained for $m_a = 25.2 \pm 11.0$ $\mu\mathrm{eV}$.

Equations

$$\begin{array}{rl}{\mathcal{L}}_{PQ}=& \frac{1}{2}|\partial \Phi {|}^2-\frac{\lambda }{4}{\left(|\Phi {|}^2-f_a^2\right)}^2-\frac{\lambda T^2}{6}|\Phi {|}^2\\ & -m_a(T{)}^2{f}_a^2[1-\cos \mathrm{Arg}(\Phi )],\end{array}$$

$$m_{a}(T)^{2}=\mathrm{min}\bigg{[}\frac{\alpha_{a}\Lambda^{4}}{f_{a}^{2}(T/\Lambda)^{n}},\,m_{a}^{2}\bigg{]},$$

$${\mathcal{L}}_{PQ}=\frac{1}{2}|\partial \Phi {|}^2-\frac{\lambda }{4}{\left(|\Phi {|}^2-f_a^2\right)}^2-\frac{\lambda T^2}{6}|\Phi {|}^2-m_a(T{)}^2{f}_a^2[1-\cos \mathrm{Arg}(\Phi )],$$

$$m_{a}(T)^{2}=\mathrm{min}\bigg{[}\frac{\alpha_{a}\Lambda^{4}}{f_{a}^{2}(T/\Lambda)^{n}},\,m_{a}\bigg{]},$$

$$m_a^2=\frac{m_{\pi }^2{f}_{\pi }^2}{f_a^2}\frac{m_u{m}_d}{(m_u+m_d{)}^2},$$

$$\displaystyle\ddot{\phi}_{1}+3H\dot{\phi}_{1}-\frac{1}{R^{2}}\nabla^{2}\phi_{1}+\frac{1}{3}\lambda\phi_{1}\bigg{[}3\left(\phi_{1}^{2}+\phi_{2}^{2}-f_{a}^{2}\right)+T^{2}\bigg{]}-\frac{m_{a}(T)^{2}\phi_{2}^{2}}{\left(\phi_{1}^{2}+\phi_{2}^{2}\right){}^{3/2}}=0$$

$$\displaystyle\ddot{\phi}_{2}+3H\dot{\phi}_{2}-\frac{1}{R^{2}}\nabla^{2}\phi_{2}+\frac{1}{3}\lambda\phi_{2}\bigg{[}3\left(\phi_{1}^{2}+\phi_{2}^{2}-f_{a}^{2}\right)+T^{2}\bigg{]}+\frac{m_{a}(T)^{2}\phi_{1}\phi_{2}}{\left(\phi_{1}^{2}+\phi_{2}^{2}\right)^{3/2}}=0\,.$$

$$\hat{\eta }=\frac{R}{R(T=T_1)}=\frac{R}{R_1}={\left(\frac{t}{t_1}\right)}^{1/2},$$

$$\displaystyle\psi_{1}^{\prime\prime}+\frac{2}{\hat{\eta}}\psi_{1}^{\prime}-\bar{\nabla}^{2}\psi_{1}+\frac{1}{H_{1}^{2}}\left[\lambda\psi_{1}\bigg{(}\hat{\eta}^{2}f_{a}^{2}\left(\psi_{1}^{2}+\psi_{2}^{2}-1\right)+\frac{1}{3}T_{1}^{2}\bigg{)}-m_{a}^{2}(T_{1})\hat{\eta}^{2}\mathrm{min}(\hat{\eta},\hat{\eta}_{c})^{n}\left(\frac{\psi_{2}^{2}}{\left(\psi_{1}^{2}+\psi_{2}^{2}\right){}^{3/2}}\right)\right]=0$$

$$\displaystyle\psi_{2}^{\prime\prime}+\frac{2}{\hat{\eta}}\psi_{2}^{\prime}-\bar{\nabla}^{2}\psi_{2}+\frac{1}{H_{1}^{2}}\left[\lambda\psi_{2}\bigg{(}\hat{\eta}^{2}f_{a}^{2}\left(\psi_{1}^{2}+\psi_{2}^{2}-1\right)+\frac{1}{3}T_{1}^{2}\bigg{)}+m_{a}^{2}(T_{1})\hat{\eta}^{2}\mathrm{min}(\hat{\eta},\hat{\eta}_{c})^{n}\left(\frac{\psi_{1}\psi_{2}}{\left(\psi_{1}^{2}+\psi_{2}^{2}\right)^{3/2}}\right)\right]=0\,,$$

$${\psi }_1^{\prime \prime }+\frac{2}{\tilde{\eta }}{\psi }_1^{\prime }-{\bar{\nabla }}^2{\psi }_1+\lambda {\psi }_1\left[{\tilde{\eta }}^2\left({\psi }_1^2+{\psi }_2^2-1\right)+\frac{T_1^2}{3f_a^2}\right]=0$$

$${\psi }_2^{\prime \prime }+\frac{2}{\tilde{\eta }}{\psi }_2^{\prime }-{\bar{\nabla }}^2{\psi }_2+\lambda {\psi }_2\left[{\tilde{\eta }}^2\left({\psi }_1^2+{\psi }_2^2-1\right)+\frac{T_1^2}{3f_a^2}\right]=0,$$

$${\left(\frac{T_1}{f_a}\right)}^2\approx 8.4\times 10^5\left(\frac{10^{12}\mathrm{GeV}}{f_a}\right).$$

$$m_{\mathrm{eff}}^2=\lambda \left(\frac{T^2}{3}-f_a^2\right).$$

$$⟨{\varphi }_i(x){\varphi }_j(y)⟩$$

$$⟨{\dot{\varphi }}_i(x){\dot{\varphi }}_j(y)⟩$$

$$⟨{\dot{\varphi }}_i(x){\varphi }_j(y)⟩$$

$$n_k=\frac{1}{e^{{\omega }_k/T}-1},{\omega }_k=\sqrt{k^2+m_{\mathrm{eff}}^2}.$$

$$⟨{\varphi }_i(k){\varphi }_j(k^{\prime })⟩$$

$$⟨{\dot{\varphi }}_i(k){\dot{\varphi }}_j(k^{\prime })⟩$$

$$⟨{\dot{\varphi }}_i(k){\varphi }_j(k^{\prime })⟩$$

$$⟨{\varphi }_i(k)⟩$$

$$⟨{\dot{\varphi }}_i(k)⟩$$

$${\psi }_1^{\prime \prime }+\frac{2}{\hat{\eta }}{\psi }_1^{\prime }-{\bar{\nabla }}^2{\psi }_1+\tilde{\lambda }{\hat{\eta }}^2{\psi }_1({\psi }_1^2+{\psi }_2^2-1)-\min (\hat{\eta },{\hat{\eta }}_c{)}^n{\hat{\eta }}^2\left(\frac{{\psi }_2^2}{({\psi }_1^2+{\psi }_2^2{)}^{3/2}}\right)=0$$

$${\psi }_2^{\prime \prime }+\frac{2}{\hat{\eta }}{\psi }_2^{\prime }-{\bar{\nabla }}^2{\psi }_2+\tilde{\lambda }{\hat{\eta }}^2{\psi }_2({\psi }_1^2+{\psi }_2^2-1)+\min (\hat{\eta },{\hat{\eta }}_c{)}^n{\hat{\eta }}^2\left(\frac{{\psi }_1{\psi }_2}{({\psi }_1^2+{\psi }_2^2{)}^{3/2}}\right)=0,$$

$$\tilde{\lambda }=\lambda {\left(\frac{f_a}{m_a(T_1)}\right)}^2$$

$$\mathcal{L}=\frac{1}{2}(\partial a{)}^2-m_a^2(T)f_a^2\left[1-\cos \left(\frac{a}{f_a}\right)\right],$$

$${\theta }^{\prime \prime }+\frac{2}{\hat{\eta }}{\theta }^{\prime }-{\bar{\nabla }}^2\theta +\min (\hat{\eta },{\hat{\eta }}_c{)}^n{\hat{\eta }}^2\sin \theta =0.$$

$${\theta }^{\prime \prime }+\frac{2}{\hat{\eta }}{\theta }^{\prime }-{\nabla }^2\theta +{\hat{\eta }}_c^n{\hat{\eta }}^2\theta =0.$$

$$\theta (\hat{\eta })=f(\hat{\eta })\exp (i\mathbf{k}\cdot \mathbf{x}).$$

$$f^{\prime \prime }(\hat{\eta })+\frac{2f^{\prime }(\hat{\eta })}{\hat{\eta }}+f(\hat{\eta })\left({\hat{\eta }}^2{\hat{\eta }}_c^n+{\mathbf{k}}^2\right)=0,$$

$$f(\hat{\eta })=\frac{\exp (-\frac{i}{2}{\hat{\eta }}^2{\hat{\eta }}_c^{n/2})}{\hat{\eta }}\left[C_1{H}_{-\frac{1}{2}{\hat{\eta }}_c^{-n/2}\left({\hat{\eta }}_c^{n/2}+i{\mathbf{k}}^2\right)}\left(\sqrt[4]{-1}\hat{\eta }{\hat{\eta }}_c^{n/4}\right)+C_2{}_1{F}_1\left(\frac{1}{4}{\hat{\eta }}_c^{-n/2}\left({\hat{\eta }}_c^{n/2}+i{\mathbf{k}}^2\right);\frac{1}{2};i{\hat{\eta }}^2{\hat{\eta }}_c^{n/2}\right)\right],$$

$$\mathcal{H}=f_a^2\left[\frac{1}{2}{\dot{\theta }}^2+\frac{1}{2R^2}(\nabla \theta {)}^2+m_a^2(1-\cos \theta )\right],$$

$$\mathcal{H}=m_a^2{f}_a^2\left[\frac{{\theta }^{\prime 2}+(\bar{\nabla }\theta {)}^2}{2{\hat{\eta }}_c^{6.68}{\hat{\eta }}^2}+(1-\cos \theta )\right],$$

$$\mathcal{H}\approx \frac{m_a^2{f}_a^2}{2}\left(\frac{{\theta }^{\prime 2}}{{\hat{\eta }}_c^n{\hat{\eta }}^2}+{\theta }^2\right),$$

Full text

††preprint: LCTP-19-08

The quantum chromodynamics (QCD) axion is a well-motivated dark-matter (DM) candidate capable of producing the present-day abundance of DM while also resolving the strong CP problem of the neutron electric dipole moment Peccei and Quinn (1977a, b); Weinberg (1978); Wilczek (1978); Preskill et al. (1983); Abbott and Sikivie (1983); Dine and Fischler (1983). The axion is an ultralight pseudo-scalar particle whose mass primarily arises from the operator $aG\tilde{G}/f_{a}$, with $a$ the axion field, $G$ the QCD field strength, $\tilde{G}$ its dual, and $f_{a}$ the axion decay constant. Below the QCD confinement scale, this operator generates a potential for the axion; when the axion minimizes this potential it dynamically removes the neutron electric dipole moment, thus solving the strong CP problem. In the process the axion acquires a mass $m_{a}\sim\Lambda_{\rm QCD}^{2}/f_{a}$, with $\Lambda_{\rm QCD}$ the QCD confinement scale. The standard ultraviolet completion of the axion low-energy effective field theory is that the axion is a pseudo-Goldstone boson of a symmetry, called the Peccei-Quinn (PQ) symmetry, which is broken at the scale $f_{a}$ Kim (1979); Shifman et al. (1980); Dine et al. (1981); Zhitnitsky (1980); Srednicki (1985).

The cosmology of the axion depends crucially on the ordering of PQ symmetry breaking and inflation. If the PQ symmetry is broken before or during inflation, then inflation produces homogeneous initial conditions for axion field and generically the cosmology is relatively straightforward Marsh (2016). In this work we focus on the more complex scenario where the PQ symmetry is broken after reheating. Immediately after PQ symmetry breaking, the initial axion field is uncorrelated on scales larger than the horizon, with neighboring Hubble patches coming into causal contact in the subsequent evolution of the Universe. This leads to complicated dynamical phenomena, such as global axion strings, domain walls, and non-linear field configurations called oscillons (also referred to as axitons) Hogan and Rees (1988); Kolb and Tkachev (1993, 1994a, 1994b); Zurek et al. (2007); Enander et al. (2017); Vaquero et al. (2018).

We perform numerical simulations to evolve the axion field from the epoch directly before PQ symmetry breaking to directly after the QCD phase transition. Once the field has entered the linear regime after the QCD phase transition, we analytically evolve the free-field axion to matter-radiation equality. The central motivations for this work are to (i) quantify the spectrum of small-scale ultracompact minihalos that emerges through the non-trivial axion self-interactions and initial conditions, and (ii) to determine the $m_{a}$ that leads to the correct DM density in this scenario.

The post-inflation PQ symmetry breaking cosmological scenario has been the subject of considerable numerical and analytic studies. It has been conjectured that this cosmology gives rise to ultra-dense compact DM minihalos with characteristic masses \sim$$10^{-13}-$10^{-11}$ $M_{\odot}$, though we show that the typical masses are actually smaller than this, and initial DM overdensities of order unity Kolb and Tkachev (1994a, 1993, b); Tinyakov et al. (2016); Davidson and Schwetz (2016); Fairbairn et al. (2018); Vaquero et al. (2018). In this work we compute the minihalo mass function precisely, combining state-of-the-art numerical simulations with a self-consistent cosmological picture. Understanding this mass function is important as it affects the ways that we look for axions in this cosmological scenario. For example, it has been claimed that microlensing by minihalos and pulsar timing surveys Dror et al. (2019) may constrain the post-inflation PQ symmetry breaking axion scenario Fairbairn et al. (2018), but these analyses rely crucially on the form of the mass function at high overdensities and masses. The axion minihalos may also impact indirect efforts to detect axion DM through radio signatures Pshirkov and Popov (2009); Huang et al. (2018); Hook et al. (2018); Safdi et al. (2018); Bai and Hamada (2018); Buschmann et al. .

A precise knowledge of the $m_{a}$ that gives the observed DM density is of critical importance for axion direct detection experiments Shokair et al. (2014); Du et al. (2018); Brubaker et al. (2017a); Al Kenany et al. (2017); Brubaker et al. (2017b); Caldwell et al. (2017); Kahn et al. (2016); Ouellet et al. (2019a, b); Chaudhuri et al. (2015); Silva-Feaver et al. (2016). We find $m_{a}=25.2\pm 11.0\,\,\mu\mathrm{eV}$, which is within range of e.g. the HAYSTAC program Brubaker et al. (2017b). Our axion mass estimate is similar to that found in recent simulations Klaer and Moore (2017) but disagrees substantially with earlier semi-analytic estimates Davis (1986); Davis and Shellard (1989); Battye and Shellard (1994); Wantz and Shellard (2010); Hiramatsu et al. (2011); Kawasaki et al. (2015); Ballesteros et al. (2017). The minihalo mass function is also important for interpreting the results of the laboratory experiments. If a large fraction of the energy density of DM is in compact minihalos, it is possible that the expected DM density at Earth is quite low or highly time dependent, which means that direct detection experiments would need to be more sensitive than previously thought or use an alternate observing strategy.

The original simulations that tried to estimate the minihalo mass function were performed in Kolb and Tkachev (1994a) on a grid of size 1003. Ref. Kolb and Tkachev (1994a) found oscillons (soliton-like oscillatory solutions) that contribute to the high-overdensity tail of the mass function. Note that oscillons are analogous to the breather solutions found in the Sine-Gordon equation (see e.g. Visinelli et al. (2018)). Recently Vaquero et al. (2018) performed updated simulations on a grid of size 81923. Our results expand on and differ from those presented in Vaquero et al. (2018) in many ways, such as through our initial state that begins before the PQ phase transition, measurement of the overall DM density, evolution to matter-radiation equality, and accounting of non-Gaussianities.

Simulation setup: We begin our simulations with a complex scalar PQ field $\Phi$, with Lagrangian

[TABLE]

with $T$ the temperature, $\lambda$ the PQ quartic coupling strength, and $m_{a}(T)$ the temperature-dependent axion mass generated by QCD Hiramatsu et al. (2012). We fix $\lambda=1$ for definiteness, though this does not affect our final results. The parametrization of the temperature-dependent mass is adopted from the leading-order term in the fit in Wantz and Shellard (2010). Explicitly, the axion mass is parametrized by

[TABLE]

for $\alpha_{a}=1.68\times 10^{-7}$, $\Lambda=400\,\mathrm{MeV}$ and $n=6.68$, though in the Supplementary Material (SM) we consider alternate parameterizations. The growth of the mass is truncated when it reaches its zero-temperature value, which occurs at $T\approx 100\,\mathrm{MeV}$ independent of the axion decay constant. The zero-temperature mass is given by $m_{a}\approx 5.707\times 10^{-5}(10^{11}\,\,{\rm GeV}/f_{a})$ eV Grilli di Cortona et al. (2016).

For the PQ-epoch simulations we begin well before the breaking of the PQ symmetry at a time when the PQ field is described by a thermal spectrum. The simulation is performed by evolving the equations of motion on a uniformly spaced grid of side-length $L_{PQ}=8000$ in units of $1/(a_{1}H_{1})$, with $a_{1}$ ($H_{1}$) the scale factor (Hubble parameter) at the temperature when $H_{1}=f_{a}$, at a resolution of $1024^{3}$ grid-sites. We use a standard leap-frog algorithm in the kick-drift-kick form with an adaptive time-step size and with the numerical Laplacian calculated by the seven-point stencil. It is convenient to use the rescaled conformal time $\tilde{\eta}=\eta/\eta_{1}$, where $\eta_{1}$ is the conformal time at which point $H(\eta_{1})\equiv H_{1}=f_{a}$. The simulation begins at $\tilde{\eta}_{i}=0.0001$ and proceeds with initial time-step $\Delta\tilde{\eta}_{i}=0.004$ until $\tilde{\eta}=250$, after which a variable time-step calculated by $\Delta\tilde{\eta}_{i}(250/\tilde{\eta})$ is used to maintain temporal resolution of the oscillating PQ fields. Convergence was tested by re-running small time intervals of the simulation at smaller time steps. The PQ fields evolve from their initial thermal configuration until the PQ phase transition occurs at $\tilde{\eta}\approx 280$, after which the radial mode $|\Phi/f_{a}|$ acquires its vacuum expectation value (VEV). We simulate until $\tilde{\eta}_{f}=800$ in order to proceed to a time at which fluctuations around the radial mode VEV have become highly damped.

Note that the difference in $\tilde{\eta}$ between $\tilde{\eta}=1$ and the PQ phase transition is proportional to $\sqrt{m_{\rm pl}/f_{a}}$, with $m_{\rm pl}$ the Planck mass. The actual choice of $f_{a}$ here does not play an important role since we evolve the axion-string network into the scaling regime. In the left panel of Fig. 1 we show the final state of our simulation at the completion of the PQ simulation. The string network is seen in yellow, with the blue colors indicating regions of higher than average axion density. The length of the simulation box at this point is around $8000/(a(\tilde{\eta}_{f})H(\tilde{\eta}_{f}))$, and we indeed find that there is around one string per Hubble patch as would be expected in the scaling regime.

We use the final state of the PQ-epoch simulation as the initial state in our QCD-epoch simulation. To do so we assume that the axion-string network remains in the scaling regime between the two phase transitions (see, e.g., Hiramatsu et al. (2011)). Recently Gorghetto et al. (2018) found evidence for a logarithmic deviation to the scaling solution and we confirm this behavior in the SM. However, we perform tests to show that this deviation to scaling likely has a minimal impact on both the minihalo mass function and on the DM density, though we still assign a systematic uncertainty to our DM density estimate from the scaling violation.

Anticipating requiring greater spatial resolution for late-times in our QCD simulation, we increased the resolution of our simulation to $2048^{3}$ grid-sites with a nearest-neighbor interpolation algorithm. We re-interpreted the physical dimensions of our box from side-length $L_{\rm PQ}=8000$ in PQ spatial units to $L_{\rm QCD}=4$ in units of $1/(a_{1}H_{1})$. These units are defined such that $H_{1}\equiv H(\eta_{1}^{\rm QCD})=m_{a}(\eta_{1}^{\rm QCD})$ at conformal time $\eta_{1}^{\rm QCD}$. Further, we use the dimensionless parameter $\hat{\eta}=\eta/\eta_{1}^{\rm QCD}$. While our PQ simulation ended at $\tilde{\eta}_{f}=800$ in PQ units, the start time in the QCD phase transition is taken to be $\hat{\eta}_{i}=0.4$ in the QCD units. Modes enter the horizon as their co-moving wavenumber becomes comparable to the co-moving horizon scale, which scales linearly with $\eta$. Therefore, by maintaining the ratio $L_{\rm PQ}/\tilde{\eta}_{f}=L_{\rm QCD}/\hat{\eta}_{i}$, we preserve the status of our modes with respect to horizon re-entry.

We then evolve the equations of motion with our initial step size now chosen to be $\Delta\hat{\eta}_{i}=0.001$. As before, we adaptively refine our time step size, using time-step $\Delta\hat{\eta}_{i}(1.8/\hat{\eta})^{3.34}$ after $\hat{\eta}=1.8$, to maintain resolution of the oscillating axion field. We simulate until $\hat{\eta}_{f}=7.0$, periodically checking if all topological defects have collapsed. When this occurs, we switch to axion-only equations of motion for computational efficiency, since past this point the radial mode does not play an important role.

The conformal time $\hat{\eta}_{c}$ at which the mass growth was cut off corresponds to the physical value of the axion decay constant since it relates the temperature $T_{1}$ at which the axion begins to oscillate and the cutoff temperature $T_{c}\approx 100$ MeV at which the axion reaches its zero-temperature mass. We performed simulations at five values of $\hat{\eta}_{c}$ uniformly spaced between $2.8$ and $3.6$. These values are chosen to access different values of $f_{a}$ while still preserving a hierarchy between $\hat{\eta}_{c}$ and our simulation end time in order to provide sufficient time for the field to relax. At each of the five values of $\hat{\eta}_{c}$, we performed simulations at five different values of the parameter $\tilde{\lambda}$, defined by $\tilde{\lambda}\equiv\lambda f_{a}^{2}/m_{a}(\hat{\eta}_{1})^{2}$. This parameter can be interpreted as the squared mass of the radial PQ mode relative to the axion mass, at conformal time $\hat{\eta}_{1}$. In order for excitations of the radial mode to be well-resolved in our simulation, we require that the resolution of our simulation $\Delta\bar{x}$, with $\bar{x}=a_{1}H_{1}x$ and $x$ the spatial coordinate, be such that $1/(\hat{\eta}\tilde{\lambda}^{1/2}\Delta\bar{x})>1$, making simulations for realistic axion parameters $\tilde{\lambda}$ impossible. We break the relation between $\tilde{\lambda}$ and $f_{a}$ and consider $\tilde{\lambda}=[1024,1448,3072,3584,5504]$ in order to study the impact of this parameter.

We illustrate three important phases of the QCD-epoch simulation in Fig. 1. The left-most panel is the initial state discussed previously in the context of the PQ-epoch simulation final-state. When $m_{a}(\hat{\eta})=3H(\hat{\eta})$ at $\hat{\eta}\approx 1.22$, strings grow longer and become less numerous, with domain walls forming on surfaces bounded by the strings. This is illustrated in the middle panel, with red colors indicating domain walls. As the temperature continues to decrease with increasing $\hat{\eta}$, strings and domain walls tighten and decrease in size until they collapse. By $\hat{\eta}\gtrsim 2.0$, the network collapses in its entirety. Shortly thereafter, we observe the formation of oscillons Kolb and Tkachev (1994a); Amin and Shirokoff (2010); Vaquero et al. (2018). We note that the oscillon field configuration is relativistic, so that near the origin of the oscillons the oscillation wavelength is \sim$$m_{a}(\hat{\eta})^{-1}, which is rapidly shrinking with increasing time. After the zero-temperature mass is reached, oscillons stop shrinking and slowly dissipate at varying rates until the full field enters the linear regime. White regions in the right-most panel of Fig. 1 denote regions of high axion energy density, which are mostly inhabited by oscillons.

At the end of the simulation, the field has relaxed into the linear regime (e.g., axion self-interactions are unimportant), but the field remains mildly relativistic because axion radiation is produced during the string-network collapse and during the oscillon collapse. It is therefore important to continue evolving the axion field until a time nearer to matter-radiation equality to allow the field to become non-relativistic everywhere and also to allow the compact but high-momentum overdensities to spread out. We perform this evolution analytically by exactly solving the linear axion equations of motion in Fourier space. We end this evolution shortly before matter-radiation equality ($T\sim$ keV), at which time proper velocities have frozen out but local radiation domination is preserved at all locations in our simulation box so that gravitational effects remain negligible.

Analysis and results.— We provide Supplementary Data Buschmann et al. (2019) containing the final state from our most realistic QCD-epoch simulation, after having performed the evolution to near matter-radiation equality. Note that the axion field after the QCD phase transition is highly non-Gaussian and phase-correlated at small scales and cannot accurately be reconstructed from the power spectrum. In fact considering that we find large overdensities $\delta$ ($\delta\sim 10$), with $\delta=(\rho-\bar{\rho})/\rho$ and $\bar{\rho}$ ($\rho$) the average (local) DM density, the field could not possibly be Gaussian at these scales, considering that Gaussian random fields have symmetric over and under-densities but under-densities with $\delta<-1$ would have negative DM density.

We may try to estimate the present-day mass function by performing a clustering analysis on the final states. In particular, we expect that the large overdensities will detach from the cosmic expansion, due to reaching locally matter-radiation equality before the rest of the Universe, and collapse onto themselves under gravity. Thus by clustering the 3-D spatial energy density distribution from the simulation slightly before matter-radiation equality and quantifying the distribution of masses and overdensities that we find, we can make predictions for the spectrum of minihalo masses and concentrations today.

From the final-state we construct an overdensity field $\delta(x)$, and we identify overdensities as closed regions of positive $\delta$. Under this definition 50% of the total mass is in overdensities. In practice, we identify these regions by first finding all positive local maxima, then recursively identifying all neighboring grid sites that are larger than 20% of the corresponding local maxima. This threshold is necessary to have a clear boundary between different overdensities, though the final mass function is not strongly dependent on the specific choice of 20%. Additionally, we discard overdensities that consist of less than 80 grid sites to avoid discretization issues in the final result. Note that we discard only about 0.8% of the total mass that would otherwise be assigned to an overdensity due to the 80 grid-site limit and 20% threshold. We assign to each overdensity a mass $M$ and a single mean concentration parameter $\delta$.

An illustration of our clustering procedure is shown in Fig. 2.

In that figure we show a 2-dimensional slice through the overdensity field for our most realistic simulation with $\hat{\eta}_{c}=3.6$ and $\tilde{\lambda}=5504$. Note that in the left panel we show the field at $\hat{\eta}=7$ at the end of the QCD simulation while in the right panel we show the same slice slightly before matter-radiation equality, denoted by $\hat{\eta}_{\rm MR}=10^{6}$ and corresponding to $T\sim$ keV. While a large overdensity left over from oscillon decay, along with corresponding rings of relativistic axion radiation, is visible in the left panel, that structure largely disperses in the subsequent evolution to $\hat{\eta}_{\rm MR}$. Two-dimensional slices through the boundaries of the clustered regions are shown in red in Fig. 2.

We characterize the minihalo mass function through the distribution $d^{2}f/d(\log M)/d\delta$, where $f$ represents the fraction of mass in overdensities of mean overdensity $\delta$ and mass $M$ with respect to the total mass in minihalos. We compute the mass function for all of the 25 simulations at varying $\tilde{\lambda}$ and $\hat{\eta}_{c}$. To perform the extrapolation to the physical $f_{a}$ ($\hat{\eta}_{c}$), we use the following procedure. First, we normalize the total DM density found in the simulation at $\hat{\eta}_{\rm MR}$ to the value that would give the observed DM density today. Then we perform the clustering algorithm to determine $d^{2}f/d(\log M)/d\delta$. We rescale all of the masses by $\left[(a_{1}H_{1})^{\text{sim}}/(a_{1}H_{1})^{\text{target}}\right]^{3}$, where $(a_{1}H_{1})^{\text{sim}}$ is the simulated horizon size at $\hat{\eta}=1$ and $(a_{1}H_{1})^{\text{target}}$ is the horizon size at the target $f_{a}$. The shift accounts for the fact that the characteristic scale of the overdensities is expected to be set by the horizon volume when the axion field begins to oscillate (see, e.g., Fairbairn et al. (2018); Vaquero et al. (2018) and the SM). The effect of this shift is to move all of the masses to lower values, since the target $m_{a}$ is larger than those we simulate. The resulting mass function for our most realistic simulation is shown in Fig. 3. As we show in the SM, after applying the mass shift the mass functions appear to give relatively consistent results between the different $\hat{\eta}_{c}$, though the agreement is not perfect at high $M$. As a result, we cannot exclude the possibility that simulating to the target $\hat{\eta}_{c}$ would give different results, especially at high masses, compared to our extrapolations. On the other hand, the effect of $\tilde{\lambda}$ appears to be minimal, since this parameter only affects the decay of the string network.

We may also compare our determinations of the total DM density produced during the QCD phase transition to previous analyses (see e.g. Davis (1986); Davis and Shellard (1989); Battye and Shellard (1994); Wantz and Shellard (2010); Hiramatsu et al. (2011); Kawasaki et al. (2015); Ballesteros et al. (2017); Klaer and Moore (2017)). Our results are summarized in Fig. 4, where we show the DM density today that we find for our top four $\hat{\eta}_{c}$, converted to $f_{a}$, for our most physical $\tilde{\lambda}$. The uncertainties in our $\rho_{a}$ measurements are determined from the variance between the different $\tilde{\lambda}$ simulations, and while some small dependence on $\tilde{\lambda}$ is expected, we find that this dependence is subdominant to statistical noise and no trend is detectable in our data. We also include a conservative $10\%$ systematic uncertainty that accounts for our unphysical fixing of the effective number of degrees of freedom $g_{*}$ throughout our simulation, a 15% systematic uncertainty from violations to scaling between the PQ and QCD phase transitions, and the uncertainty on the measured value of $\Omega_{a}$ in our Universe Aghanim et al. (2018) (see the SM for details).

In Fig. 4 we compare our results to the best-fit simulation result from Klaer and Moore (2017), which like us numerically evolved the axion-string system through the QCD phase transition, albeit with a different formalism, and also the semi-analytic calculations from Kawasaki et al. (2015). Our results are in reasonable agreement with those in Klaer and Moore (2017) and significantly disagree with those in Kawasaki et al. (2015). Note that we self-consistently account for all production mechanisms for axion DM in our simulation, including string decay in the few decades before the QCD phase-transition. It is the late-time axion production, right before the QCD phase transition, which is most important since it is the least redshifted Kawasaki et al. (2015). The source of the discrepancy could be due in part to the fact that by artificially separating the production mechanisms, Kawasaki et al. (2015) over-counted the DM density produced (see Klaer and Moore (2017)). Additionally, the highly non-linear axion dynamics at the QCD epoch likely violate the number-conserving assumptions made by Kawasaki et al. (2015).

We may estimate the $f_{a}$ that gives the correct DM density by fitting our results to a power-law $\Omega_{a}\sim f_{a}^{\alpha}$. We find the best-fit index $\alpha=1.24\pm 0.04$, only including statistical uncertainties, which is marginally compatible with the analytic calculations in Kawasaki et al. (2015); Klaer and Moore (2017) that predict $\alpha=(n+6)/(n+4)\approx 1.187$. Fixing $\alpha$ to the theoretical value, we find $\Omega_{a}=(0.102\pm 0.02)\times(f_{a}/10^{11}\mathrm{GeV})^{1.187}$, now incorporating the correlated systematic uncertainties, which leads to the prediction that the correct DM density is achieved for $f_{a}=(2.27\pm 0.33)\times 10^{11}\,\,\mathrm{GeV}$ ($m_{a}=25.2\pm 3.6\,\,\mu\mathrm{eV}$) in agreement with Klaer and Moore (2017). Note that if we fit for $\alpha$ instead of fixing $\alpha$ to the theoretical value we find $m_{a}=17.4\pm 4.5\,\mu\mathrm{eV}$; the difference between the two $m_{a}$ estimates could be due to a systematic difference between the theoretically predicted $\alpha$ and the actual dependence of $\Omega_{a}$ on $f_{a}$. In light of this we use the difference between the two $m_{a}$ estimates as an estimate of the systematic uncertainty from the extrapolation to $f_{a}$ below those simulated. We additionally include a $\sim$27% uncertainty on $m_{a}$ from uncertainties in the mass growth of the axion (see the SM for details), leading to the prediction $m_{a}=25.2\pm 11.0$ $\mu$eV.

Discussion.— We performed high-resolution simulations of axion DM in the cosmological scenario where the PQ symmetry is broken after inflation, starting from the epoch before the PQ phase transition and evolving the field until matter-radiation equality. After matter-radiation equality one should still evolve the axion field gravitationally down to lower redshifts, which we plan to do in future work. Our mass function is an estimate of the resulting mass function one would find after simulating the gravitational collapse. It is possible that the true halos will be slightly larger in mass due to e.g. accretion of surrounding DM.

We may try to estimate the halo sizes based upon when we expect the halos to collapse gravitationally. Under the assumption, for example, that the final density profile is a constant-density sphere of radius $R$ (which is likely not a good approximation but still is useful to get a sense of the halo sizes), then the halo density today was argued to be approximately $\rho\approx 140\rho_{\rm eq}\delta^{3}(\delta+1)$, where $\rho_{\rm eq}$ is the DM density at matter-radiation equality Kolb and Tkachev (1994b). This implies, for example, that a $M=10^{-14}$ $M_{\odot}$ subhalo with an initial average overdensity $\delta=3$ will have a characteristic size of \sim$$1\times 10^{6} km. The implications for direct and indirect axion detection efforts (e.g., non-trivial time dependence) are likely substantial and will be the subject of future work. One immediate implication, however, is that the axion minihalos are likely out of reach for microlensing and pulsar timing surveys Dror et al. (2019), given the small minihalo masses.

*We thank Jiji Fan for collaboration on early stages of this project, and we thank Asimina Arvanitaki, Masha Baryakhtar, Gus Evrard, Andrew Long, Nicholas Rodd, Jesse Thaler, Ken Van Tilburg, and Kathryn Zurek for useful comments and discussion. This work was supported in part by the DOE Early Career Grant de-sc0019225. JF received additional support from a Leinweber Graduate Fellowship. This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. *

I Simulation Equations of Motion

Our phenomenological Lagrangian describing the PQ field is adopted from the construction of Hiramatsu et al. (2012) and is of the form

[TABLE]

where $\Phi$ is the complex PQ scalar, $T$ is the temperature, $\lambda$ is the PQ quartic coupling strength, $f_{a}$ is the PQ-scale identified as the axion decay constant, and $m_{a}(T)$ is the temperature-dependent axion mass Hiramatsu et al. (2012). The parametrization of the temperature-dependent mass is adopted from the leading order term in the fit in Wantz and Shellard (2010). Explicitly, the axion mass is parametrized by

[TABLE]

for $\alpha=1.68\times 10^{-7}$, $\Lambda=400\,\mathrm{MeV}$ and $n=6.68$. The growth of the mass is truncated at $T\approx 100\,\mathrm{MeV}$. The zero-temperature mass is given by

[TABLE]

where $m_{\pi}$ is the pion mass, $f_{\pi}$ is the pion decay constant, $m_{u/d}$ is the up/down quark mass. Details of the temperature-dependent axion mass, or equivalently, the topological susceptibility, remain uncertain, especially at low temperatures. Note that we do not explore here how our results are affected by uncertainties in the temperature-dependent axion mass, though doing so is a worthwhile direction for future work.

Decomposing the complex scalar as $\Phi=\phi_{1}+i\phi_{2}$, and assuming a radiation-dominated cosmological background, leads to equations of motion in metric coordinates of the form

[TABLE]

Over temperatures $T\gtrsim 100\,\mathrm{MeV}$, the number of relativistic degrees of freedom $g_{*}$ in the Standard Model is expected to vary only mildly. For simplicity, we therefore assume $g_{*}=81$, which is a typical value adopted at high temperatures (though later in the SM we explore the systematic uncertainty introduced by this assumption). It is useful to define a dimensionless conformal time $\hat{\eta}$ such that

[TABLE]

where $R$ is the scale factor and the time $t_{1}$ (with $T(t_{1})\equiv T_{1}$) is a reference time that will be defined differently in the PQ and QCD epoch simulations.

The axion-mass term is not included in our PQ-epoch simulations. In our QCD-epoch simulations, on the other hand, the mass term is included and drives the dynamics. In this case, the mass grows until the cutoff temperature $T_{c}$ at which point the axion mass reaches its zero-temperature value; the corresponding conformal time is given by $\hat{\eta}_{c}=R(T=T_{c})/R_{1}$. Rewriting (S5) with the dimensionless coordinates, we then find

[TABLE]

where $\phi=f_{a}\psi$, primes denote derivatives with respect to $\hat{\eta}$, and the spatial gradient is taken with respect to $\bar{x}=a_{1}H_{1}x$.

I.1 The PQ Epoch

Simulations in the PQ epoch occur at $T\gg\Lambda_{\rm QCD}$ and so the temperature-dependent axion mass may be neglected. We therefore take our equations of motion to be

[TABLE]

and we fix $\tilde{\eta}=1$ to be the time at which $H_{1}=f_{a}$. Note that for our PQ-epoch simulations we refer to $\hat{\eta}$, defined in (S6), as $\tilde{\eta}$ in order to avoid confusion with the dimensionless conformal time $\hat{\eta}$ used in the QCD-epoch simulations. The ratio $(T_{1}/f_{a})^{2}$ is determined by

[TABLE]

In principle, it would seem that axions of different decay constants would require different simulations in the PQ epoch. However, this ratio is degenerate with our choice of physical box size and dynamical range in $\tilde{\eta}$ in a particular simulation, allowing us to perform only one PQ simulation and interpret its output as the initial state of the axion field for several different values of $f_{a}$. The key assumption behind this, however, is that at late times after the PQ phase transition the field enters the scaling regime so that we may reinterpret the output of the PQ simulation in the appropriately rescaled box as the initial state of the QCD simulation at much lower temperatures. Note that the value of $\lambda$ is a free parameter, which we naturally choose to be $\lambda=1$ though it has little effect.

I.1.1 Initial Conditions for a PQ Scalar

We generate initial conditions for our PQ scalar by taking it to be described by a thermal distribution characterized by the temperature $T$ at the initial early time. As can be read off from the Lagrangian, each of the two fields has an effective mass of the form

[TABLE]

Correlation functions of the initially-free massive scalar fields are given by

[TABLE]

where overdots denote differentiation with respect to time, and we have defined

[TABLE]

In momentum space, these correlation functions take the form

[TABLE]

Our simulations occur on a discrete lattice of finite size, so the correlation functions above lead to initial conditions set by a realization of a Gaussian random field specified in Fourier space by

[TABLE]

Note that we include the 50 lowest $k$-modes in each of the three directions when constructing the initial conditions, and we have verified that including more modes does not affect our results.

I.2 Early Times in the QCD Epoch

During the QCD epoch, $T\sim\Lambda_{\rm QCD}$, and so the axion mass is non-negligible. Here, we define $\hat{\eta}=1$ to be the time at which $H_{1}=m_{a}(T_{1})$, with the axion field beginning to oscillate shortly thereafter when $m_{a}=3H$. The equations of motion are then given by

[TABLE]

where we have neglected the $T_{1}$ contribution to the PQ scalar mass as it is small compared to $f_{a}$. The parameter $\tilde{\lambda}$ is defined by

[TABLE]

and can be interpreted as the squared mass of the radial mode $|\Phi/f_{a}|$. For physical parameters we expect $\tilde{\lambda}\gg 1$, though in practice we find that the final results are relatively independent of $\tilde{\lambda}$ for moderately sized values of the parameter, as described in the main text and later in the SM. Indeed, our choices for $\tilde{\lambda}$ allow us to resolve the radial mode mass by more than a few grid-spacings, satisfying the requirement of Gorghetto et al. (2018) to accurately study the axion spectrum from string radiation, unlike Vaquero et al. (2018). There exist additional criteria on the largeness of $\tilde{\lambda}$ such that the metastability of topological defects is preserved despite the unphysical smallness of simulated $\tilde{\lambda}$ in comparison with the rapidly increasing axion mass. At all times prior to expected defect collapse, our choices of $\tilde{\lambda}$ satisfy the simplest construction of these conditions Vaquero et al. (2018), with our choice of $\tilde{\lambda}=5504$ satisfying the most stringent criteria established in Fleury and Moore (2016). We note that we are largely unable to differentiate between simulations at any two particular values of $\tilde{\lambda}$, and that our choice of values appear to have minimal impact, as illustrated further below.

I.3 Late Times in the QCD Epoch

The presence of topological defects in the axion field at early times during the QCD epoch requires that we fully simulate both degrees of freedom of the PQ field. Once the topological defects have collapsed, however, we are free to use the axion-only equations of motion. Our axion is defined by $a=f_{a}\mathrm{arctan2}(\phi_{1},\phi_{2})$ and has the Lagrangian

[TABLE]

along with corresponding equations of motion

[TABLE]

Above, we define $\theta=a/f_{a}$. Evolving these equations of motion is formally equivalent to freezing out excitations of the radial mode by taking $\tilde{\lambda}\rightarrow\infty$, which more accurately recovers the true physics of the evolution of the axion field for realistic values of $f_{a}$. Note that the coordinate $\bar{x}$ and $\hat{\eta}$ here are identical to those used in evolving the two degrees of freedom of the complex scalar performed prior to defect collapse.

I.4 Analytically Evolving in the Fixed-Mass Small-Field Limit

At late times when the axion mass has reached its zero-temperature value and the axion field has redshifted considerably so that $|\theta|\ll 1$, the equations of motion are linear and well-approximated by

[TABLE]

We may solve this equation analytically by going to Fourier space and adopting an ansatz for the solution as

[TABLE]

This ansatz leads to the equation

[TABLE]

which has the general solution

[TABLE]

for coefficients $C_{1}$ and $C_{2}$ determined by boundary conditions, and where $H_{n}$ and $\,{}_{1}F_{1}$ are the analytic continuations of the Hermite polynomials and the confluent hypergeometric function of the first kind, respectively. From this analytic solution, we can transfer late-time field configurations from our simulation to arbitrary large $\hat{\eta}$. Differentiation with respect to $\hat{\eta}$ may be straightforwardly performed to find $f^{\prime}(\hat{\eta})$ at large $\hat{\eta}$ as well. The computation of the analytically continued Hermite polynomials and hypergeometric functions was performed with the python package mpmath.

We directly compare the differential mass spectrum at $\hat{\eta}=7$ with the same field analytically evolved to $\hat{\eta}=\eta_{\text{MR}}$ in Fig. S1. While the basic differential shape is the same, the $\hat{\eta}=7$ results have a much wider distribution in $\delta$. In particular, all overdensities above $\delta>10$ have vanished by the time matter-radiation equality is reached. However, the peak of the distribution is still around $\delta=1$. Evolving the fields down to matter-radiation equality is important because many of the modes are generated with high momentum at the QCD epoch, causing the large overdensities to disperse by the time of matter-radiation equality.

II Studying the (Over)Density Field

Our interest in this work is studying the energy density field $\rho$ and the overdensity field $\delta=(\rho-\bar{\rho})/\rho$ realized in the axion field from our simulations. The axion energy density for the axion field $a=f_{a}\theta$ is computed by the Hamiltonian density

[TABLE]

which can be rewritten in simulation units as

[TABLE]

assuming $\hat{\eta}>\hat{\eta}_{c}$. At late times, the Hamiltonian is approximately

[TABLE]

when all modes in the simulation are non-relativistic and the field values are small.

II.1 Oscillons

Large overdensities right after the QCD phase transition are caused by oscillons. Oscillons are, in contrast to strings and domain walls, not topological defects but arise due to non-linearities in the equation of motions, forming at locations where the the axion self-interaction dominates the Hubble friction. As a result, the first oscillons form at the location of collapsed strings and domain walls, where the axion remains excited and reaches large field values. However, at later times, oscillons are observed forming throughout the simulation box. The dynamics of the oscillons are highly non-trivial, especially as the axion self-interaction increases in strength with the growing axion mass.

Oscillons decrease in size over time following the oscillation wavelength \sim$$m_{a}(T)^{-1}, as axions in the core are relativistic. Good spatial resolution is therefore needed to resolve them. In order to study their behavior we perform a 2D (two spatial dimensions, one time) simulation using the same simulation setup in the PQ- and QCD-epoch as in 3D. We find that there is no qualitative difference between 2D and 3D simulations regarding oscillons, but going to 2D allows us to increase the spatial resolution to $4096^{2}$ grid sites and to subsequently increase $\hat{\eta}_{c}$.

We illustrate the evolution of an oscillon in Fig. S2. Two scenarios are considered with different truncation points of the mass growth, $\hat{\eta}_{c}=4.0$ and $\hat{\eta}_{c}=6.0$. Note how the radius of the oscillon decreases as long as $m_{a}(T)$ is increasing. The circles in S2 have radius $m_{a}(T)^{-1}$, and the oscillon cores are seen to track this scale. Subsequently, if the mass growth is truncated at $\hat{\eta}_{c}=4.0$, the radius of the oscillon is constant as well. When the mass growth is cut-off, the density contrast at the core of the oscillon slowly decreases over time and the oscillons dissipate, as can be seen in the two lower right panels in Fig. S2.

II.2 Calculating the Axion Relic Abundance

To calculate the axion DM abundance as a function of $m_{a}$, we first need to understand the relationship between the mass cutoff conformal time $\hat{\eta}_{c}$ and the decay constant $f_{a}$. Here we use the relation $T_{1}/\hat{\eta}_{c}=T_{c}$, with $T_{c}\approx 100$ MeV. This allows us to solve for $f_{a}$ in terms of $\hat{\eta}_{c}$. The energy densities are calculated from the axion field and its derivatives according to (S33) after numerically evolving until $\hat{\eta}=7$, then analytically evolving until $\hat{\eta}_{\rm MR}=10^{6}$, at which point the contribution of the gradient term to the energy density is negligible. As a side note, our definition of $\hat{\eta}_{\rm MR}$ actually puts us at slightly earlier times than global matter-radiation equality. This us because matter-radiation equality is, locally, reached earlier for the largest overdensities and because we want to make sure that gravitational interactions can be neglected. In particular, note that the temperature corresponding to $\hat{\eta}_{\rm MR}$ is given by $T_{\rm MR}=T_{c}\eta_{c}/\hat{\eta}_{\rm MR}$. For our most realistic simulation with $\eta_{c}=3.6$ this corresponds to $T_{\rm MR}\approx 0.5$ keV. However, if we reinterpret the final state for a more realistic axion with $m_{a}\approx 25$ $\mu$eV, which has a higher $\eta_{c}$, then $T_{\rm MR}\approx 4$ keV. In practice, though, the exact value of $T_{\rm MR}$ is not important because by these temperatures the proper motions in the axion field are frozen out and the field is thus not evolving non trivially. As a consequence our results (both for the DM density and for the spectrum of overdensities) are not sensitive to small (or even relatively large) changes to the exact value of $\hat{\eta}$ that we evolve to.

Note that we present our results in terms of the DM density fraction today $\Omega_{a}$, which is defined as the ratio of the average energy density today in DM relative to the observed critical energy density. We compute statistical error bars at each value of $f_{a}$ from the variance as a function of $\tilde{\lambda}$ at fixed $\hat{\eta}_{c}$. We note that no trend is visible in the data for the dependence of $\Omega_{a}$ on $\tilde{\lambda}$, as is shown in Fig. S3.

The statistical noise is inferred from the spread in $\Omega_{a}$ values, which are determined from the output at $\eta_{\rm MR}$, between different $\tilde{\lambda}$. The observed variations are consistent with the expected noise from Poisson counting statistics due to having a finite number of overdensities within the simulation box.

In Fig. S4 we show our results for $\Omega_{a}$ as a function of $f_{a}$, compared to earlier predictions in Kawasaki et al. (2015) and Klaer and Moore (2017).

For reference, we also include predictions for the relic abundance based on the field value and the time derivative at $\hat{\eta}=7$. Here it is less straightforward to determine the DM axion abundance, relative to taking the results at $\hat{\eta}_{\rm MR}$, as some of the modes in the simulation are still relativistic. This introduces an additional systematic uncertainty, since the field is not completely red-shifting like radiation at this time. For these reasons it is important to evolve the field until it is completely non-relativistic before measuring the DM density.

Because the ratio of the axion mass density to entropy density is constant after the axions have become non-relativistic and the number of axions is conserved, we can redshift our energy density from our matter-radiation equality $\hat{\eta}_{\rm MR}$ to today. Then, we compare this energy density to the most up-to-date measurement of the average DM density in the Universe today $\rho_{\rm DM}=33.5\pm 0.6$ $M_{\odot}/{\rm kpc}^{3}$ Aghanim et al. (2018). Note that we have propagated all cosmological uncertainties other than those on $N_{\rm eff}$, which we have fixed to the Standard Model value. These cosmological uncertainties introduce an approximately $3\%$ correlated uncertainty across the results of our simulations. We additionally have an approximately $8\%$ uncertainty due to our assumption of fixed $g_{*}$, which is examined in greater detail later in this Supplement. These uncertainties are the dominant ones in our results, and we emphasize that they have not been typically considered in determinations of the DM axion mass. From the $\Omega_{a}$ data, for the various $f_{a}$ values simulated, we may extrapolate to predict the $f_{a}$ for an axion which produces the observed DM relic abundance by fitting a simple power law relation of the form

[TABLE]

as discussed in the main body of this work. Note that we expect $\alpha=(6+n)/(4+n)$, where $n$ is the index of the axion mass growth. We assume this scaling is valid to make our estimate for the $m_{a}$ that gives the correct DM abundance. The relation between $\alpha$ and $n$ is expected to arise for the following reason. Let us estimate the axion DM density from an axion with a constant initial misalignment angle $\theta_{i}$. The present-day axion abundance as produced by the misalignment mechanism can be estimated by

[TABLE]

where $T_{0}$ is the present-day temperature, $T_{3}$ is the temperature at which the axion began to oscillate ($m_{a}(T_{3})=3H(T_{3})$), and $g_{*}(T)$ the number of effective degrees of freedom at temperature $T$.

The initial axion abundance $\rho_{a}(T_{1})$ is given

[TABLE]

Anharmonicity factors can be included, but have no temperature or $f_{a}$ dependence. The temperature $T_{3}$ depends on $f_{a}$ through the relation $T_{3}\propto f_{a}^{-2/(4+n)}$. Substituting these relations in and keeping only terms which depend on $f_{a}$, we have

[TABLE]

We thus expect the relic abundance to scale with $f_{a}$ like $\rho_{a}\propto f_{a}^{(6+n)/(4+n)}$. Note that the DM abundance from string and domain wall production is calculated similarly in Kawasaki et al. (2015), and although our results are not consistent with those presented in that work, the abundance calculation they present proceeds similarly, yielding string and domain wall production that scale like $f_{a}^{(6+n)/(4+n)}$ as well.

On the other hand, we may also calculate the the $m_{a}$ that gives the correct DM abundance by using our fit value for $\alpha$, as defined in (S34), instead of the theoretical value. Doing so leads to a slightly lower $m_{a}$ estimate, as described in the main text.

II.3 Tests of the Overdensity Field Gaussianity

In typical cosmological contexts, overdensity fields are treated under the assumption that they are Gaussian random fields. For a real-space Gaussian field, we may Fourier transform the field and find that the squared magnitude of each mode is independently exponentially distributed with mean set by the power-spectrum and with the phase of each mode independently uniformly distributed on $[0,2\pi)$ Adler (2010). For reference, in Fig. S5 we show our power spectra $\Delta_{k}^{2}$ at fixed $\tilde{\lambda}$ across our various choices for $\hat{\eta}_{c}$. Note that we construct the power spectra from the fields that have been evolved until $\hat{\eta}=\hat{\eta}_{\rm MR}$. However, as we demonstrate below, the power spectrum fails to accurately describe the overdensity field we realize in our simulations because the field is highly non-Gaussian at small scales. As a result, standard tools for predicting structure formation that rely upon an underlying Gaussian overdensity field, such as the Press-Schechter formalism, cannot be applied to predict the spectrum of structures that form from the overdensities in the axion field, at least on the very smallest scales.

First, we note that the largest field values taken within the overdensity fields at the state realized by the analytic evolution until $\hat{\eta}=\hat{\eta}_{\rm MR}$ are $\mathcal{O}(10)$, whereas the minimum value the overdensity field can take is $-1$ by construction. This is trivially incompatible with the interpretation of the overdensity field as a Gaussian random field, which would have symmetric variance about its mean of [math]. For our overdensity fields to realize $\mathcal{O}(10)$ maxima with $-1$ as a construction-imposed minimum, there must exist considerable phase-correlations between Fourier modes, contrary to the uncorrelated phases of a Gaussian random field.

We also may inspect the distribution of power at each mode in the Fourier transformed overdensity field. If the overdensity field were Gaussian, then the power in each mode would be exponentially distributed with mean set by the value of the mean power spectrum. To test this, we plot the probability distribution $dP/dx$ of $x=|\hat{\delta}(\mathbf{k})|^{2}/\langle|\hat{\delta}(\mathbf{k})|^{2}\rangle_{|\mathbf{k}|=k}$, with $\hat{\delta}({\bf k})$ the Fourier transformed overdensity field at momentum ${\bf k}$, as measured in the final states of our field at $\hat{\eta}=\hat{\eta}_{\rm MR}$. We compare the observed distributions with the expected Gaussian random field assumption of an exponential distribution with unit mean in Fig. S6. Dramatic deviations from the expected behavior are observed for large $|\mathbf{k}|$. We stress, however, that in addition to these distributions departing from the expected exponential distributions, the real and imaginary components across modes are also highly phase correlated on small scales.

II.4 Minihalo Mass Spectrum

In this subsection we give additional details and results for the minihalo mass and density spectrum. In addition to the technical difficulties associated with a non-Gaussian overdensity field, computational limitations prevent us from performing realistic simulations of $f_{a}\sim 10^{11}\,\mathrm{GeV}$ axions, which would require us to simulate until $\hat{\eta}_{c}\approx 15$. We instead interpret our simulation results at smaller $\hat{\eta}_{c}$ in appropriate units to rescale these results to the target $f_{a}\approx 2\times 10^{11}\,\mathrm{GeV}$. We do so with the following methods. The total axion mass contained within some set of grid-sites in our simulation can be computed from the Hamiltonian as

[TABLE]

where $\bar{\rho}$ is computed by the average of our Hamiltonian in (S33) in the simulation box. We calculate $H_{1}$ from $T_{1}$ based on our choice of $f_{a}$, then rescale $\bar{\rho}$ to the value of the axion energy density at the time $\hat{\eta}$ such that the correct relic abundance is realized today. In this manner, we aim to rescale all dimension-full quantities related to $f_{a}$ to our target $f_{a}$. In particular, we rescale the DM density $\bar{\rho}$ to give the correct DM density realized in our Universe, and we also rescale the minihalo masses by the factor $\propto(a_{1}H_{1})^{-3}$ appearing in (S38) to those for the target $f_{a}$.

We illustrate the rescaling procedure in Fig. S7. In that figure we show the differential mass distribution of minihalos $df/d\log M$ as a function of minihalo mass $M$. These mass distributions have been rescaled such that $\bar{\rho}$ matches the actual DM density. However, the solid curves do not have the $a_{1}H_{1}$ Hubble volume rescaling included. The dashed curves, on the other hand, apply the Hubble volume rescaling factor but for a target $\hat{\eta}_{c}$ of $\hat{\eta}_{c}=3.6$, which is that corresponding to the black curve. Clearly there are still differences between the rescaled dashed curves and the black curve, which tells us that there are dynamical effects that arise from changing $\hat{\eta}_{c}$ that are not captured by the simple rescaling. This should not be too surprising considering that e.g. the mass growth affects the oscillon stability, which determines the high-mass part of the distribution. In our work we rescale the mass function to the target $f_{a}$ as described above, but it is important to keep in mind that this almost certainly results in a systematic uncertainty from the fact that we do not capture the full oscillon dynamical range in doing so. Also note that all of the mass functions abruptly drop off at low halo masses. This is due to our resolution limit on the finite lattice. We also cannot rule out the possibility that the low-mass tail continues down to much smaller masses.

We compare the single-differential mass fractions for different values of $\hat{\eta}_{c}$ and $\tilde{\lambda}$ as a function of $\delta$ and $M$ in Fig. S8. Note that here we have applied the rescaling factors for the masses to our true target $f_{a}$, which is that which gives the correct DM density. First of all we note that there is no dependence on $\tilde{\lambda}$ visible in our parameter range within than statistical scatter. As for the differential distribution as a function of $\delta$, there is also no clear dependence on $\hat{\eta}_{c}$ visible. The only place where a clear dependence on $\hat{\eta}_{c}$ is visible is in the mass fraction as a function of $M$. Here, the peak values shift to smaller masses upon increasing $\hat{\eta}_{c}$, even after including the rescaling factors.

It is useful to have an approximate analytic formula for the differential mass fraction. We find that the differential mass fraction as function of $\delta$ can be accurately described by a Crystal Ball function based on a generalized Gaussian and a power-law high-end tail together with a suppression factor at high-$\delta$:

[TABLE]

The parameters $B$ and $C$ are given by

[TABLE]

and they are chosen such that $df/d\delta$ and its first derivative are continuous. $A$ is not a free parameter as $\int_{0}^{\infty}d\delta(df/d\delta)=1$ must hold. The fit parameters from our most realistic simulation with $\hat{\eta}_{c}=3.6$ and $\tilde{\lambda}=5504$ are given by

[TABLE]

This fit allows us to make a precise comparison with previous work by Kolb and Tkachev Kolb and Tkachev (1996). We present in Fig. S9 the cumulative mass fraction that is in overdensities larger than $\delta_{0}$,

[TABLE]

Unsurprisingly, we find considerably less mass in highly concentrated overdensities relative to Kolb and Tkachev (1996). Whereas Kolb and Tkachev (1996) predicts roughly $10\%$ of the mass is in overdensities with $\delta=10$ or more, we find a similar result only when using the simulation output at $\hat{\eta}=7$. Once evolved to matter-radiation equality, that percentage falls to \sim$$0.1%.

III Equations of Motion for Varying Relativistic Degrees of Freedom

In this section we investigate the systematic effect on our results from the assumption of fixed $g_{*}$. In truth, the value of $g_{*}$ is not fixed at $g_{*}\approx 81$ but instead varies rather sharply in the temperature range of interest; it varies from as large as roughly $100$ to as little as roughly $10$ around the time of the QCD phase transition. This does not represent a dire shortcoming of our simulation procedure, however, as varying $g_{*}$ should only nontrivially affect the dynamics of the axion field during times when axion number density is not a conserved quantity. By $\hat{\eta}\approx 3$ most of the field has become linear, except for the isolated oscillon configurations, which means that the axion number density is mostly conserved at this time and beyond. The variation in $g_{*}$ before $\hat{\eta}\approx 3$ for our target $f_{a}$ is relatively minor. To quantify this impact, however, we perform 2D simulations accommodating the varying $g_{*}$.

With a change of variable we may rewrite the axion equations of motion, in the two-field formalism during the QCD epoch, as

[TABLE]

where we define $\tilde{\lambda}=\lambda f_{a}^{2}/H_{1}^{2}$ as before. Citing standard references Kolb and Turner (1990), we have

[TABLE]

Using these relations, we may compute

[TABLE]

Above, we have defined

[TABLE]

where the right hand side is evaluated at the time $t$ corresponding to the conformal time $\hat{\eta}$. Similarly, we evaluate

[TABLE]

for

[TABLE]

Finally, we define $f_{3}(\hat{\eta})=m_{a}(\hat{\eta}(T))^{2}/H_{1}^{2}$. Combining these results, the equations of motion take the form

[TABLE]

In the single-field formalism, these results are analogously applied to obtain

[TABLE]

In Fig. S10 we show the functions $f_{1}$, $f_{2}$, and $f_{3}$ entering into the equations of motion as functions of $\hat{\eta}$. Note that for $f_{3}$, we normalize against $\tilde{f}_{3}(\hat{\eta})$, which we define to be $f_{3}$ but with a fixed $g_{*}$. In the absence of a varying $g_{*}$, all of the curves appearing in Fig. S10 would be identically one.

To test the impact of an evolving $g_{*}$, we adopt the parametrization of $g_{*}$ from Wantz and Shellard (2010) and simulate, in two spatial dimensions, for an axion with $f_{a}=4.83\times 10^{15}\,\mathrm{GeV}$. When we assumed $g_{*}$ was constant, this axion reached its zero-temperature mass at $\hat{\eta}=3.6$, but accounting for the changing value of $g_{*}$, the axion now reaches its zero-temperature mass at $\hat{\eta}\approx 5.5$. As before, we conclude our simulation at $\hat{\eta}=7$, and we calculate a relic abundance that is $7.7\%$ smaller than it is in the fixed $g_{*}$ case. We note that this scenario represents something of a worst-case scenario for the impact of $g_{*}$ on the dynamics because $g_{*}$ varies significantly during the epoch where axion number density is not conserved for this choice of $f_{a}$, and so we adopt this as a quantification of the systematic error associated with adopting a fixed $g_{*}$. For our target $f_{a}$, $g_{*}$ varies less, relative to the example illustrated, when the axion is in the non-linear regime and so we expect the effect of varying $g_{*}$ to be less important in this case.

IV Testing the Impact of the Mass Parametrization

Precise details regarding the temperature dependence of the axion mass remain uncertain. While we have chosen to use the parametrization of Wantz and Shellard (2010) with index $n=6.68$ as done in Hiramatsu et al. (2011); Kawasaki et al. (2015), an alternative result is provided in Borsanyi et al. (2016). In that work, an index of $n\approx 8.2$ is found at high temperatures, though we do note that an increasingly shallow dependence on $T$ is realized at lower temperatures. Motivated by power-law fits to this numerical result and informed by considerations of the changing number of degrees of freedom, recent works have taken an index of $n=7.6$ in Klaer and Moore (2017) and $n=7.3$ in Vaquero et al. (2018) to study the axion field. In this section, we use the extreme value $n=8.2$ to estimate the maximal effect that uncertainties in the mass growth may have on the determination of the DM density.

We perform simulations for different $n$ in two spatial dimensions. This is done for computational efficiency, and while we do not expect such a simplification to significantly affect the main conclusions we caution some care should be taken when interpreting these results for this reason. We fix $\tilde{\lambda}=5504$ and $f_{a}\approx 4.8\times 10^{14}\,\mathrm{GeV}$. This choice of $f_{a}$ corresponds to $\hat{\eta}_{c}=3.6$ in the $n=6.68$ parametrization. However, the value of $\hat{\eta}_{c}$ depends on our choice of $n$ and is $\hat{\eta}_{c}\approx 3.1$ for the choice of $n=8.2$, since the mass grows faster in that case. We re-simulate with this alternative choice of index until $\hat{\eta}=7$ and then recompute the present-day axion abundance by analytically transferring the simulation fields to the same late time physical temperature. We find that there is a \sim$$10\% enhancement in the expected relic abundance with $n=8.2$ versus $n=6.68$. This is somewhat surprising, considering that the analytic estimate predicts that higher $n$ should result in a lower DM abundance at fixed $f_{a}$. To understand how this result affects the final determination of the axion mass, we fit the predicted scaling $\Omega_{a}\sim f_{a}^{(n+6)/(n+4)}$ for the DM abundance using $n=8.2$ and find the $m_{a}$ that gives the correct DM abundance. The result is that with $n=8.2$ we find that the $m_{a}$ that gives the correct DM abundance is enhanced by $\sim$27% compared to the $n=6.68$ case. We account for this 27% uncertainty as an additional systematic uncertainty in our final determination of the axion mass.

V Testing Deviations from the String Scaling Regime

While our simulation was performed in two stages, it can be understood as a single simulation in which the PQ phase transition and the beginning of the QCD phase transition are separated by approximately an order of magnitude in temperature. By comparison, for a physically motivated hierarchy, we would expect these two epochs in our simulation to be separated by at least 11 orders of magnitude in temperature. As a result, our simulations might be expected to be highly unphysical. However, it has been conjectured that the axion field and associated defect network enters a scaling regime some time after the PQ phase transition (see, e.g., Hiramatsu et al. (2011)). If this conjecture is true and our field configuration has entered the scaling regime before the axion begins to oscillate, our simulation should be expected to give a good description of the physics of interest despite the abbreviated hierarchy.

Recent work has found evidence for logarithmic deviations to the number of strings per Hubble patch in the scaling regime Gorghetto et al. (2018). In this section we confirm that we also observe such deviations. This implies that we are not fully justified in taking the final state of our PQ-epoch simulation, fast-forwarding through the rest of the radiation dominated epoch to the QCD phase transition, and then restarting our simulation directly before the QCD phase transition. This is because the axion-string network should change logarithmically during the evolution between the phase transitions. Below, we provide evidence for the logarithmic deviation to scaling and then perform simulations to address the impact of this deviation on our determination of the axion mass $m_{a}$ and the spectrum of DM minihalos.

The average number of strings per Hubble patch is commonly defined by Hiramatsu et al. (2011); Fleury and Moore (2016); Gorghetto et al. (2018)

[TABLE]

where $l(\tilde{\eta})$, $L(\tilde{\eta})$, and $t(\tilde{\eta})$ are the physical total string length in the box, the physical box length, and the physical time, each a function of $\tilde{\eta}$, respectively. We measure the string length by first identifying grid sites that are next to a string. This is achieved by forming a loop in each of the three dimensions around a test grid site. The grid site is flagged once at least one change larger than $\pi$ in the axion field between consecutive grid sites is found. In a 2D slice this implies the 4 closest grid sites that surround the string core are tagged, such that we use the number of tagged grid sites divided by 4 as a measure for the string length. Note that this is a rough estimate for the string length and more sophisticated methods exist Hiramatsu et al. (2011).

We compute $\xi(\tilde{\eta})$ at 13 points in $\tilde{\eta}$ in our PQ simulation, with results illustrated in Fig. S11. As in Gorghetto et al. (2018), we find that $\xi$ depends logarithmically on $\tilde{\eta}$ after the PQ phase transition. Note that the shaded region denotes $\tilde{\eta}$ before the PQ phase transition, where it does not make sense to talk about axion strings. We fit the model

[TABLE]

to the $\{\tilde{\eta},\xi\}$ data, where $T_{PQ}$ denotes the temperature of the PQ phase transition, and we find $\alpha\approx-2.60$ and $\beta\approx 1.27$. Note that our values for $\xi$ at comparable $\tilde{\eta}$ are significantly larger than those found in Gorghetto et al. (2018). Part of this discrepancy could be due to the way in which we measure string length versus in that work, which may introduce an overall rescaling between our two results. We are in good agreement, however, with Gorghetto et al. (2018) regarding the distribution of string length in long and short strings. As in that work, we find that approximately $80\%$ of the string length resides in long strings, much larger than a Hubble length, at all times in our simulation as seen in Fig. S12.

Since we do observe logarithmic scaling violations, it is important to determine the impact of these corrections on the minihalo mass spectrum and the DM relic abundance. In particular, we find that $\xi$ should be around a factor of 15 higher at the QCD phase transition than it is for the final state of our most realistic PQ-epoch simulation. The string density $\xi$ at the beginning of the QCD simulation depends on the simulation box size as $\xi_{\rm QCD}\propto L_{\rm QCD}^{-2}$, where $L_{\rm QCD}$ is the box size in units of $1/(a_{1}H_{1})$ when $H=m_{a}$, for a fixed initial state. Thus by performing new simulations with the same initial conditions and run parameters as our fiducial analysis (namely $\hat{\eta}_{c}=3.6,\tilde{\lambda}=5504$, and starting at $\hat{\eta}_{i}=0.4$) but changing $L_{\rm QCD}$ from 4 to $L_{\rm QCD}=3$ and $L_{\rm QCD}=2$, we may enhance $\xi_{\rm QCD}$ by a factor of $16/9$ and $4$, respectively, compared to our fiducial simulation. While these value still fall short of the physically motivated enhancement \sim$$15, such simulations still allow us to see if there is a trend in how $\xi$ affects observables such as the DM density. We do caution that modifying $L_{\rm QCD}$ in this way is somewhat unphysical as it changes the horizon entry status of modes in the simulation box from the end of the PQ simulation as compared to the beginning of the QCD simulation.

The results of varying $L_{\rm QCD}$ in order to modify $\xi$ are shown in Fig. S13, where we see no discernible trends in the dependence of the relic abundance on $\xi$. Note that the uncertainties in Fig. S13 are estimates of the statistical uncertainty. As the box size gets smaller the statistical uncertainty increases. However, we caution that these are estimates only, as we have not run multiple independent simulations for each box size due to computational limitations. It is possible that the true uncertainties at small box sizes are larger than indicated due to the fact that there are a small number of e.g. domain walls that form these cases. Still, to be maximally conservative given the available datasets, we estimate the difference between the $L_{\rm QCD}=2$ and $L_{\rm QCD}=3$ values for $\Omega_{a}$ as a systematic uncertainty induced from the deviation to scaling. However, we cannot be sure that this difference is not a result of statistics or from the way in which we artificially mock up initial conditions with higher $\xi$ values. The systematic uncertainty we assign from these tests is $15\%$ correlated between different $f_{a}$ points.

We show in Fig. S14 the impact on the minihalo mass spectrum. Again, all simulations largely agree within their error bands (estimated from statistical uncertainties), indicating that an increase in $\xi$ has only a marginal effect on the late-time axion field. Note that computational resources limit us to just these three additional simulations, and we leave a more detailed investigation of the dependence of $\Omega_{a}$ on $\xi$ to future work.