2026 lattice simulation estimate for post-inflationary QCD-axion dark matter

Use natural units and the source's spatially flat FLRW model in conformal time: $\mathcal{X}=I_{\tau }\times {\mathbb{T}}_L^3$, with $I_{\tau }\subseteq \mathbb{R}$ an open interval, $L>0$, ${\mathbb{T}}_L^3=(\mathbb{R}\mathrm{/}L\mathbb{Z}{)}^3$, and metric $g_{\mu \nu }=a(\tau {)}^2{\eta }_{\mu \nu }$ for a specified $a\in C^2(I_{\tau },(0,\mathrm{\infty }))$ and mostly-minus Minkowski metric $\eta$. The minimal field content is the Standard Model plus one complex scalar $\mathrm{\Phi }\in C^2(\mathcal{X},\mathbb{C})$; the simulation evolves the discretized values ${\mathrm{\Phi }}_{\mathrm{\ell }}(\tau )\in \mathbb{C}$ for $\mathrm{\ell }\in (\mathbb{Z}\mathrm{/}\mathrm{12,288}\mathbb{Z}{)}^3$, with periodic boundary conditions. Let $T\in C^1(I_{\tau },(0,\mathrm{\infty }))$ be the radiation-era temperature. Let $U(1{)}_{\mathrm{P}\mathrm{Q}}=\{e^{i\alpha }:\alpha \in \mathbb{R}\mathrm{/}(2\pi \mathbb{Z})\}$ act globally, not as a gauge symmetry, by $\mathrm{\Phi }\mapsto e^{i\alpha }\mathrm{\Phi }$. Every other field $X_j$ has a declared integer PQ charge $q_j\in \mathbb{Z}$ and action $X_j\mapsto e^{i{q}_j\alpha }{X}_j$; the colored-field content and charges are model inputs and are not fixed uniquely by the simulation. With $\theta =\arg \mathrm{\Phi }(\mathrm{m}\mathrm{o}\mathrm{d}2\pi )$ where $\mathrm{\mid }\mathrm{\Phi }\mathrm{\mid }>0$, let $\chi \in C^0([0,\mathrm{\infty }),[0,\mathrm{\infty }))$ be the adopted temperature-dependent QCD topological susceptibility, including its zero-temperature input. Its convention is $\chi (T):={\mathrm{\partial }}_{\vartheta }^2\mathcal{F}(T,\vartheta ){\mathrm{\mid }}_{\vartheta =0}$, of mass dimension four, where $\vartheta$ is the dimensionless QCD vacuum angle and $\mathcal{F}$ is the QCD equilibrium free-energy density, assumed twice differentiable in $\vartheta$ near zero. For $T>0$, fix its angle-dependent part by $\mathcal{F}(T,\vartheta )-\mathcal{F}(T,0)=-{\lim }_{V\to \mathrm{\infty }}(T\mathrm{/}V)\log [Z_V(T,\vartheta )\mathrm{/}{Z}_V(T,0)]$, where $Z_V$ is the supplied finite-spatial-volume QCD partition function and the limit and angle derivatives are assumed to exist; at $T=0$, $\mathcal{F}$ denotes the vacuum-energy density. Continuity at zero is an adopted input, not a result of this simulation. The model potential is $V(\mathrm{\Phi },T)=\frac{\lambda }{4}(\mathrm{\mid }\mathrm{\Phi }{\mathrm{\mid }}^2-f_a^2{)}^2+\chi (T)\frac{\mathrm{\mid }\mathrm{\Phi }{\mathrm{\mid }}^2}{f_a^2}(1-\cos \theta )$, with $\lambda >0$ dimensionless and symmetry-breaking scale $f_a$ of mass dimension one supplied as model parameters; at $\mathrm{\Phi }=0$ the QCD term is defined by its continuous extension. The QCD color anomaly explicitly breaks the continuous PQ symmetry. Define $N_{\mathrm{D}\mathrm{W}}$ as the number of distinct minima of the QCD-induced angular potential in one fundamental $2\pi$ period, equivalently the absolute anomaly coefficient in this normalization. The simulation takes $N_{\mathrm{D}\mathrm{W}}=1$. A PQ string is a line defect where $\mathrm{\Phi }=0$ and the phase winds by $2\pi$ around a transverse loop. The QCD-induced potential creates angular minima; a domain wall is a surface across which $\theta$ interpolates between adjacent minima. With $N_{\mathrm{D}\mathrm{W}}=1$, one wall attaches to each string, and their coupled configuration is the string-wall system whose collapse is simulated.

For small angular displacements near $\theta =0$, define the temperature-dependent harmonic mass by $m_a(T{)}^2:=\chi (T)\mathrm{/}{f}_a^2$; let $H(\tau ):=a^{\mathrm{\prime }}(\tau )\mathrm{/}a(\tau {)}^2$ be the physical Hubble rate in conformal time. Let $H_0$ be the present expansion rate, $h:=H_0\mathrm{/}(100\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}\mathrm{M}\mathrm{p}{\mathrm{c}}^{-1})$, $G_N$ Newton's constant, ${\rho }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},0}:=3H_0^2\mathrm{/}(8\pi G_N)$ the present critical density, and ${\mathrm{\Omega }}_{c,0}:={\rho }_{c,0}\mathrm{/}{\rho }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},0}$ the present cold-dark-matter density parameter. Let ${\tau }_{\ast }$ be the source's time at which $m_a(T_{\ast })=H_{\ast }$, with $H_{\ast }=H({\tau }_{\ast })$. For $k>0$, use specifically the averaged estimator of Correia et al., §III Eqs. (21)–(22): $N_{\mathrm{a}\mathrm{x}}(k,\tau ):=[{\mathcal{P}}_0(k,\tau )\mathrm{/}{\mathrm{\Omega }}_k(\tau )+{\mathrm{\Omega }}_k(\tau ){\mathcal{P}}_s(k,\tau )\mathrm{/}{k}^2]\mathrm{/}[2a(\tau )]$, with ${\mathrm{\Omega }}_k(\tau ):=[k^2+a(\tau {)}^2{m}_a(T(\tau ){)}^2{]}^{1\mathrm{/}2}$. Here ${\mathcal{P}}_b=k^3⟨\mathrm{\mid }{\tilde{J}}_b{\mathrm{\mid }}^2{⟩}_{\widehat{k}}\mathrm{/}(2{\pi }^2{V}_{\mathrm{b}\mathrm{o}\mathrm{x}})$, $b\in \{0,s\}$, is power per logarithmic wavenumber in the source's box Fourier convention, ${\tilde{J}}_{\mu }(\mathbf{k})={\int }_{{\mathbb{T}}_L^3}{e}^{-i\mathbf{k}\cdot \mathbf{x}}{J}_{\mu }(\mathbf{x})d^{3}x$, $V_{\mathrm{b}\mathrm{o}\mathrm{x}}=L^3$, $\mathbf{k}\in (2\pi \mathrm{/}L){\mathbb{Z}}^3$. The currents are $J_{\mu }=\mathrm{Im}({\mathrm{\Phi }}^{\ast }{\mathrm{\partial }}_{\mu }\mathrm{\Phi })\mathrm{/}\mathrm{\mid }\mathrm{\Phi }\mathrm{\mid }$ off the string cores, and ${\tilde{J}}_s={\widehat{k}}_i{\tilde{J}}_i$; $⟨\cdot {⟩}_{\widehat{k}}$ denotes the angular-shell average (finite Fourier-shell binning on the lattice). Core values do not affect continuum integrals when the zero set has measure zero; no phase derivative at a core is asserted. This is the mean of the temporal and longitudinal estimators, not the temporal-only alternative of Appendix D. Write $n_{\mathrm{a}\mathrm{x}}(\tau )={\int }_0^{\mathrm{\infty }}{N}_{\mathrm{a}\mathrm{x}}(k,\tau )d\ln k$ for the continuum spectral notation, assuming integrability; the finite simulation uses its nonzero-mode shell bins and lattice cutoffs. Its particle-number interpretation assumes the source's late-time free-field/adiabatic regime, not an exact conserved count during defect collapse. The dimensionless comoving number-density parameter is ${\nu }_{\mathrm{a}\mathrm{x}}(\tau ):=[a(\tau )\mathrm{/}a({\tau }_{\ast }){]}^3{n}_{\mathrm{a}\mathrm{x}}(\tau )\mathrm{/}(f_a^2{H}_{\ast })$. The paper's ${\nu }_{\mathrm{a}\mathrm{x}\ast }$ is the asymptotic post-transition value of this parameter; it is estimated from finite-time simulation data by the stated scaling extrapolation, rather than measured at infinite time.

Correia et al. report classical-field simulations of this minimal post-inflationary model using ${\mathrm{12,288}}^3$ lattice sites. The initial momentum field is zero and the initial scalar field is a unit-variance Gaussian random field, with comoving correlation length selected from $l_{\mathrm{\Phi }}\in \{10,20,40,80\}$ at $\tau =50$ in simulation units. A diffusion stage runs to $\tau =70$, followed by core growth to $\tau =470$ and then physical evolution. Under their scaling model, five runs with different random seeds per parameter combination give the number-density estimate and its run-level standard error; the same five seeds are reused across parameter combinations, so those combinations cannot be treated as independent replicates; a power-law fit and extrapolation yield their final ${\nu }_{\mathrm{a}\mathrm{x}\ast }=9.48(23)$. Combining this with the selected QCD topological-susceptibility data yields the model inference $m_a=17.46(84)\mu \mathrm{e}\mathrm{V}$ when axions are assumed to supply all of the present cold-dark-matter density ${\rho }_{c,0}$. The conversion adopts Correia et al.’s Planck-attributed input ${\mathrm{\Omega }}_{c,0}{h}^2=0.1200(11)$, the CMB temperature $T_0=2.725(1)\mathrm{K}$, the zero-temperature susceptibility $\chi (0)=(75.5\pm 0.5{)}^4{\mathrm{M}\mathrm{e}\mathrm{V}}^4$, and the paper's selected journal-version susceptibility and thermodynamic tables. The parenthetical mass uncertainty is the paper's symmetrized propagated uncertainty for its selected inputs and scaling analysis; it does not cover the between-dataset susceptibility disagreement or alternative network extrapolations discussed by the authors. Planck 2018 VI, Table 2 and Eq. (23), instead give ${\mathrm{\Omega }}_{c,0}{h}^2=0.1200\pm 0.0012$ as a 68% interval for base-$\mathrm{\Lambda }$CDM with TT,TE,EE+lowE+lensing and the baseline Plik likelihood. The cited simulation adopts $0.1200(11)$; the inspected passage does not explain the difference in its uncertainty, so it is retained as the authors’ adopted conversion input rather than identified with that Planck interval. Likewise the inspected scaling-fit description gives least-squares fits, weighting and marginalization over a Gaussian scaling-density input, but does not specify a covariance treatment for the matched seeds across parameter combinations. The quoted error is therefore reported with that unresolved dependence limitation; no independent covariance-corrected uncertainty or revised mass is calculated here. This is a numerical-model inference, not a measured axion mass or a string-theory prediction.