Thermal axion production in QCD from the lattice

Technical Summary: Thermal Axion Production in QCD from the Lattice

Problem Statement
The production rate of thermal axions in the early Universe is governed by the QCD topological charge density spectral function, $\sigma(k, k)$, evaluated on the light-cone. While perturbative methods work at very high temperatures and chiral perturbation theory applies below the transition region, the crossover region ($T \sim 150\text{--}300$ MeV) is strongly coupled, rendering these techniques unreliable. Standard lattice QCD approaches face a fundamental obstacle: they compute Euclidean correlation functions, whereas the production rate requires the Minkowski spectral function. Previous attempts to bridge this gap relied on analytical continuation, a process prone to large systematic errors due to the ill-posed nature of the inversion. This paper addresses the challenge of computing weighted integrals of the axion production rate directly from Euclidean lattice data without analytical continuation.

Methodology
The authors propose a method inspired by recent work on photon production rates (Harvey Meyer et al.), which allows for the computation of weighted integrals of spectral functions using Euclidean correlators.

1. Theoretical Framework:

   • The axion production rate is related to the spectral function $\sigma(\omega, \omega)$.
   • The authors define a Euclidean correlator $H_E(\omega_n)$ involving imaginary spatial momentum, which relates to the spectral function via a kernel integral.
   • Direct evaluation of $H_E(\omega_n)$ is problematic because the kernel grows exponentially with distance, and the integral is divergent in the ultraviolet (UV) and potentially at high frequencies.
   • To resolve divergence issues, the authors construct linear combinations of Matsubara modes, specifically $H_E^{(0,1)}$ and $H_E^{(0,1,2)}$. These combinations introduce subtraction kernels that cancel the leading UV divergences and ensure convergence at high frequencies.

2. Lattice Implementation:

   • Setup: The study is performed in pure $SU(3)$ Yang-Mills theory (quenched approximation) at a temperature $T = 1.5 T_c$ on a lattice with $N_\tau = 16$ (aspect ratio 4).
   • Gradient Flow: To handle the multiplicative renormalization and short-distance noise of the topological charge density operator $q(x)$, the authors employ improved (Zeuthen) gradient flow. This smears the operator over a radius $r \sim \sqrt{8t_f}$, stabilizing the topological charge.
   • Tail Fitting: The Euclidean correlator $G(\omega_n, r)$ decays exponentially at large distances. However, the weight function $\sinh(\omega_n r)/(\omega_n r)$ grows exponentially. To compute the integral, the authors fit the "tail" of the correlator (large $r$) to a sum of exponentials ($e^{-mr}/r$) derived from the spectral representation of the theory in a compactified geometry. The integral is then computed by combining precise lattice data at small $r$ with the fitted function at large $r$.
   • Extrapolation: Results are obtained at finite flow times $t_f$. The authors argue that for $t_f \gtrsim 0.4a^2$, the smearing function is effectively $O(4)$ invariant, implying that the flow-time dependence of the integrated quantity $H_E$ should be negligible (scaling as $t_f^2$ or higher). They perform a constant fit to the data across different flow times to extrapolate to $t_f \to 0$.

Key Contributions

• Analytical Continuation-Free Method: The paper demonstrates a viable pathway to extract weighted integrals of the axion production rate directly from Euclidean lattice simulations, bypassing the need for unstable analytical continuation.
• Handling Exponential Noise: The authors develop a robust numerical strategy combining gradient flow with tail-fitting techniques to overcome the signal-to-noise ratio degradation caused by the exponentially growing weight function in the integral.
• Subtraction Schemes: The formulation of specific linear combinations of Matsubara modes ($H_E^{(0,1)}$ and $H_E^{(0,1,2)}$) that render the integrals finite and well-defined on the lattice.

Results

• Numerical Values: In pure-glue QCD at $T = 1.5 T_c$, the authors successfully extract:
  • $H_E^{(0,1)} = -0.0357(58) T^4$
  • $H_E^{(0,1,2)} = -0.0152(91) T^4$
    The negative sign is consistent with the anti-hermitian nature of the topological charge density operator.
• Comparison with Perturbation Theory: The results are compared with Next-to-Leading Order (NLO) perturbative calculations. The lattice results are consistent with NLO predictions (using $\alpha_s \approx 0.3$), whereas Leading Order (LO) results are significantly too small.
• Feasibility: The study confirms that extracting at least one, and possibly two, finite weighted integrals is feasible with modest numerical effort (approx. 60,000 CPU core-hours) on a single lattice spacing.

Significance and Claims
The paper claims to open a new avenue for determining thermal axion production rates in the strongly coupled QCD regime where other methods fail. By successfully computing these integrals in pure-glue theory, the authors establish the feasibility of the approach. They posit that this method can be extended to full QCD at temperatures around and slightly above the crossover ($T \in [200, 1000]$ MeV), a range critical for cosmological constraints on axion dark matter.

The authors remain modest regarding the immediate phenomenological impact, noting that their method yields integral moments rather than the differential production rate $d\Gamma/d^3k$ as a function of momentum. They acknowledge that translating these integral constraints into precise limits on axion mass and coupling requires further development. Furthermore, they highlight that the success of the method relies on the lightest flavorless $0^{-+}$ state being sufficiently heavy compared to the temperature scale to ensure the correlator decays faster than the growing weight function; this condition may become challenging at very high temperatures where the coupling becomes weak.