Dark matter from axion strings with adaptive mesh refinement

An outstanding puzzle of the Standard Model is the extreme smallness of the neutron electric dipole moment. Axions, originally introduced to solve this strong-CP problem, may also constitute the dark matter. Laboratory searches are hindered by the poorly constrained axion mass. In post-inflation Peccei–Quinn (PQ) symmetry breaking, a unique axion mass yields the observed dark matter abundance, but accurately determining this mass is challenging due to the formation and evolution of axion strings, which radiate axions as the network evolves. The role of string-induced axions in the total dark matter density has been debated. Numerical studies since the 1980s, culminating in modern static-lattice simulations with up to ~8000^3 sites, have produced disparate predictions for the axion mass. This work addresses the problem using adaptive mesh refinement (AMR) to efficiently resolve string cores and extend the dynamical range, aiming to robustly determine the axion radiation spectrum from strings and the resulting dark matter abundance.

Early simulations of axion strings were limited to small static lattices (~150^3) and suggested differing contributions of string-produced axions to dark matter. Subsequent analytic and numerical studies developed network evolution models (e.g., velocity-dependent one-scale) and larger simulations. Recent large static-grid works (e.g., Gorghetto, Hardy, Villadoro) pushed to ~4500^3 sites and log(m_r/H) ~ 7.9, reporting linear growth in the string length per Hubble volume parameter ξ with log(m_r/H) and suggesting a time-evolving emission index q > 1. However, predictions for the axion mass remained widely spread (25–>500 μeV). Other studies examined domain wall dynamics and axion production near the QCD phase transition. AMR methodologies have been proposed for global string simulations to better capture string cores while reducing computational cost elsewhere. This work builds on these developments, directly comparing to and testing claims from static-grid simulations, especially regarding the growth of ξ and the axion emission spectrum index q.

The axion is the phase of a complex PQ field Φ=(r+f_a)/√2 e^{ia/f_a}. The heavy radial mode r has mass m_r=√2 f_a (setting λ=1) and is non-dynamical for T< m_r; the axion is effectively massless until the QCD epoch. Simulations are performed in a radiation-dominated background using conformal time η=R/R1=(t/t1)^{1/2}, with H1=f_a, starting at η_i=0.1 and ending at η_f=75.7. The comoving periodic box has side length L=120/(R1 H1), corresponding to ~1200 Hubble lengths at η_i and ~1.6 at η_f. Equations of motion for the rescaled fields are integrated with a third-order strong-stability-preserving Runge–Kutta (SSPRK3) method, satisfying the CFL condition. Spatial derivatives use second-order finite differences. The AMR implementation is based on AMReX (block-structured AMR). The simulation starts on a uniform 2048^3 grid (coarse level). Additional refinement levels are added when the comoving string width l (set by m_r^{-1}) falls below a target number of cells at the finest level: levels added at η≈3, 6, 12, 24 resolve l by four cells; a sixth level at η≈64 resolves l by three cells (corresponding to log(m_r/H) ≈ 2.6, 3.9, 5.3, 6.7, and 8.7). Each level halves Δx and Δt relative to the coarser level (Δx_ℓ=Δx_0/2^ℓ, Δt_ℓ=Δt_0/2^ℓ; Δt_0≈0.02). Subcycling in time is used with fourth-order spatial and second-order temporal interpolation for synchronization. Refinement is localized using tagging criteria: (i) string cores identified via plaquette piercing (Fleury & Moore algorithm), (ii) large gradients in the axion field (Δx ∇^2ψ > 0.04), and (iii) control of potential instabilities from short-wavelength radiation in ψ and r (Δx(∇^2ψ+∇^2r) > 4 at coarse level). Buffer zones of 11 cells around tagged regions ensure strings remain away from coarse-fine boundaries; regridding occurs every Δt_ℓ=0.2^ℓ Δt_0. Thermal initial conditions are generated with modes up to k_max=25 (details per Buschmann et al. 2020). The computation used NERSC Cori (1024 KNL nodes; ~5.2 million CPU hours), with auxiliary tests on Perlmutter GPUs. String network diagnostics: The total string length l is obtained by counting string-pierced plaquettes; ξ = l t^2 / V is evaluated at Hubble-spaced times Δlog(m_r/H)=log 2. ξ versus log(m_r/H) is fit with ξ = c_{−2}/log^2 + c_{−1}/log + c_0 + c_1 log, with data-driven uncertainties modeled to scale with the number of Hubble patches. Effective string tension μ is extracted from the partition of total energy density into axion, radial, and string components (masking highest AMR regions), and compared to μ_th ≈ π f_a^2 log(m_r/H) up to offsets around AMR level additions. Axion radiation spectrum: The instantaneous emission spectrum F(k/H) ∝ d log ρ_a / d(k/H) is computed from the time derivative of a screened axion field â_scr(x)=f(x) ȧ(x), with screening choices f(x)= (1+r/f_a)^2 (fiducial), 1+r/f_a, or 1. Fourier transforms use HACC SWFFT; spectra are binned in 1774 k-bins up to k/(2π)=1024√3/L. Instantaneous spectra are estimated with finite differences in time using Δlog(m_r/H)=0.25 (and tested up to log 2). Power-law fits F(k/H) ∝ (k/H)^{−q} are performed over k_IR=χ_IR H (fiducial χ_IR=50; varied 30–100) to k_UV=m_r/χ_UV (fiducial χ_UV=16; varied down to 8) and analyzed with Gaussian likelihoods including nuisance terms for combined statistical/systematic uncertainties. The evolution q(t)=q_0+q_1 log(m_r/H) is tested versus constant q. From F(k/H), moments ⟨H/k⟩ and ⟨(H/k)^2⟩ are computed numerically up to k/H=50 and analytically extended using the fitted power law to the UV cutoff. These define parameters δ via ⟨H/k⟩=δ/ξ (or δ log dependence for exact q=1) and β via ⟨(H/k)^2⟩=(β/ξ^2). The axion number density from strings at the onset of the QCD epoch (log ξ ~ 60–70) is estimated using number conservation arguments, with conservative treatment of possible additional radiation during string collapse after m_a turns on. Convergence tests compare AMR to an over-resolved static grid on a smaller box; sensitivity tests vary initial k_max (5 vs 25) and analysis cuts, showing robustness.

- AMR enables >3 orders of magnitude improvement in dynamical range relative to prior static-grid approaches; achieving equivalent finest resolution on a uniform grid would require ~65,536^3 sites. - String length per Hubble volume: ξ grows linearly with log(m_r/H) at late times, indicating a logarithmic violation of perfect scaling. Fit over log(m_r/H) ∈ (4,9): c1 = 0.254 ± 0.002; a late-time linear fit yields c1 ≈ 0.252, consistent with Gorghetto et al. 2021. Extrapolating to the QCD onset log(m_r/H) ~ 60–70 gives ξ ≈ 13–17 (≈15). - Axion radiation spectrum: The instantaneous emission index is q = 1.02 ± 0.03 (fiducial), consistent with a scale-invariant spectrum (q=1) within ~5%. A fit allowing q(t)=q_0+q_1 log(m_r/H) yields q_1 = −0.04 ± 0.08 (no evidence for growth), contrasting with earlier claims (q_1 ≈ 0.053 ± 0.005). Systematic variations in k_IR (30–100), k_UV (m_r/8–m_r/32), time step (Δlog up to log 2), and screening method leave this conclusion unchanged; low string-core resolution biases q high. - Moments of the spectrum: For q=1.06 (max at 1σ), β = 840 ± 70 and δ = 113 ± 7; for exact scale invariance q=1, δ = δ1 log(m_r/H) with δ1 = 6.2 ± 0.4; for q=0.98 at log=70, δ ≈ 820 ± 50. Nonlinearities during the QCD epoch are marginal: ⟨(α/f_a)^2⟩ ≲ 1.1 for log ≤ 70, implying number-changing effects are at most modest (~15%). - Effective string tension: Measured growth consistent with μ ∝ log(m_r/H), with fitted μ_1 = 3.7 ± 0.5 compared to μ_1 = π, indicating no significant numerical energy leakage. - Dark matter abundance: Using the measured q range (0.98–1.06) and conservative treatment of possible additional radiation during string collapse, the PQ scale consistent with Ω_DM is f_a ≈ (3.1×10^10 – 1.4×10^11) GeV, implying m_a ≈ (40–180) μeV. For q=1 exactly, m_a = 65 ± 6 μeV. Axions from domain wall/misalignment dynamics during the QCD epoch provide a sub-dominant contribution (~0.017 h^{-2} (f_a/4.3×10^{10} GeV)^{1.17}). - Robustness and convergence: AMR results agree with over-resolved static-grid runs on smaller volumes; results are insensitive to initial-state high-k cutoff (k_max=5 vs 25) and analysis choices. Low core resolution (resolving l by ~1 cell) biases spectra toward larger q at late times. - Experimental implication: Preferred mass range lies above current ADMX/HAYSTAC coverage but is accessible to future ADMX runs, MADMAX, and plasma haloscopes.

The AMR simulations demonstrate that axion strings radiate with an approximately scale-invariant spectrum over a wide dynamical range, leading to a robust prediction for the axion mass in the post-inflation PQ scenario. The linear growth of ξ with log(m_r/H) aligns with expectations from the logarithmically increasing string tension and contradicts models predicting constant ξ at late times. The absence of evidence for a time-growing q addresses a key discrepancy with some static-grid studies and underscores the importance of resolving string cores adequately: insufficient core resolution can bias the inferred spectrum. These findings imply that axion radiation from strings prior to the QCD phase transition likely dominates the axion number density and thus the dark matter abundance, sharpening predictions for the axion mass window. The results directly inform experimental strategies, suggesting a focus on m_a ≈ 40–180 μeV, with a central value near 65 μeV if the spectrum is exactly scale invariant. The methodological advances with AMR pave the way for more precise simulations through the QCD epoch, potentially reducing the dominant uncertainty from q and refining dark matter mass predictions.

This work presents the largest, highest-resolution simulations of axion string networks to date using adaptive mesh refinement, enabling precise characterization of network evolution and axion radiation. Key contributions include: (i) establishing linear growth of ξ with log(m_r/H) and validating the expected logarithmic growth of string tension; (ii) measuring an axion emission spectrum index q ≈ 1 with no evidence for time evolution; and (iii) deriving a sharpened axion mass prediction m_a ≈ 40–180 μeV (≈65 ± 6 μeV for q=1) for post-inflation PQ breaking. The AMR framework is well-suited to extend simulations into the QCD epoch, study domain wall collapse and potential shifts in the emission spectrum, and investigate related signals such as gravitational waves from global string networks. Future work should focus on increasing statistics to reduce the uncertainty on q, simulating the mass turn-on and network collapse in full detail, exploring scenarios with different domain wall numbers, and performing cross-method validation to further solidify the mass prediction and guide experimental searches.

- Extrapolation: Results are extrapolated from simulated log(m_r/H) up to ~9 to the QCD epoch at log ~ 60–70; while supported by theoretical arguments and consistency checks, uncertainties from long extrapolation remain. - QCD epoch not simulated: The simulations do not include the axion mass turn-on and full string-domain wall collapse dynamics; contributions from this epoch are conservatively bounded but not explicitly computed. - Finite-volume/time and resolution: Although volume and late-time resolution are carefully managed (AMR maintaining ≥3–4 cells across string cores), ultra-late-time core resolution and potential boundary/systematic effects may still introduce small biases. - Model assumptions: Assumes N_DW=1 (unstable walls), m_r ~ f_a (though TeV-scale m_r scenario argued to be consistent), radiation-dominated cosmology, and specific treatments of screening and spectral fitting; while tested for robustness, residual systematics may persist. - Data availability: Full simulation data (>50 TB) are not publicly hosted, limiting external reanalysis; uncertainties in some fitted quantities are modeled in a data-driven way due to difficulty in direct error estimation.