Theory

Axion dark matter

The QCD axion was originally proposed to resolve the strong-\(CP\) problem in quantum chromodynamics. More general axion-like particles (ALPs) arise in various extensions of the Standard Model. Axions and ALPs are well-motivated dark-matter candidates; for brevity this document uses “axion” for both.

In the Standard Halo Model (SHM) the local axion dark-matter field behaves as a classical, coherently oscillating pseudoscalar:

\[a(\mathbf{r}, t) = a_0 \cos\!\bigl(2\pi\nu_a t - \mathbf{p}_a \cdot \mathbf{r}/\hbar + \phi\bigr),\]

where \(\phi\) is a random phase, \(\mathbf{p}_a\) is the axion momentum, and the field amplitude is set by the local dark-matter energy density \(\rho_\mathrm{DM}\):

\[a_0 = \sqrt{\frac{2\hbar^3 \rho_\mathrm{DM}}{m_a^2 c}}.\]

The axion Compton frequency \(\nu_a = m_a c^2 / h\) is the fundamental frequency of oscillation. The field is coherent on a timescale

\[\tau_a \sim \frac{Q_a}{\pi\,\nu_a}, \qquad Q_a = \left(\frac{c}{v_\mathrm{lab}}\right)^2 \approx 10^6,\]

where \(v_\mathrm{lab} \approx 233\,\mathrm{km\,s^{-1}}\) is the speed of the laboratory in the galactic rest frame.

Axion–spin interaction

Axion field gradients couple to fermionic spins through a pseudoscalar interaction. The Hamiltonian is

\[\mathcal{H}_\mathrm{int} = -\,g_{aNN}\,\nabla a \cdot \mathbf{S},\]

where \(g_{aNN}\) is the axion–nucleon coupling strength (GeV\(^{-1}\)) and \(\mathbf{S}\) is the spin operator. This has the same form as the Zeeman interaction

\[\mathcal{H} = -\hbar\gamma\,\mathbf{B} \cdot \mathbf{S},\]

so the axion field gradient acts as an effective pseudomagnetic field

\[\mathbf{B}_\mathrm{axion} = \frac{g_{aNN}}{\gamma}\,\nabla a(\mathbf{r}, t).\]

This field drives resonant transitions when \(\nu_a\) equals the nuclear Larmor frequency \(\nu_L = \gamma B_0 / (2\pi)\).

The rms Rabi frequency in the perpendicular-gradient coupling geometry (grad_perp) is

\[\Omega_\mathrm{rms} = \frac{1}{2}\,g_{aNN}\sqrt{2\hbar c\,\rho_\mathrm{DM}}\;v_\mathrm{lab}.\]

Bloch equations

The nuclear spin magnetization \(\mathbf{M}\) evolves under the total effective field \(\mathbf{B}_\mathrm{eff}\) (bias field plus pseudomagnetic field) according to

\[\frac{dM_x}{dt} = \gamma\left(\mathbf{M}\times\mathbf{B}_\mathrm{eff}\right)_x - \frac{M_x}{T_2},\]

\[\frac{dM_y}{dt} = \gamma\left(\mathbf{M}\times\mathbf{B}_\mathrm{eff}\right)_y - \frac{M_y}{T_2},\]

\[\frac{dM_z}{dt} = \gamma\left(\mathbf{M}\times\mathbf{B}_\mathrm{eff}\right)_z - \frac{M_z - M_0}{T_1},\]

where \(T_1\) is the longitudinal (spin–lattice) relaxation time, \(T_2\) is the intrinsic transverse relaxation time, and \(M_0\) is the equilibrium magnetization.

Rotating coordinate frame

The bias field \(B_0\) is typically much stronger than the axion pseudomagnetic field. To remove the fast Larmor precession at \(\nu_L = \gamma B_0 / (2\pi)\) and enable efficient numerical integration of the slower axion-induced dynamics, axionbloch transforms to the rotating coordinate frame (RCF) that co-rotates with \(B_0\).

Field inhomogeneity and \(T_2^*\)

In a real magnet the bias field is not perfectly uniform. The ensemble of nuclear spins in a sample experiences a spread of Larmor frequencies, causing additional dephasing characterised by \(T_2^* < T_2\). axionbloch models this by sampling \(N_\mathrm{pt}\) discrete values from a Lorentzian field distribution with a user-specified fractional FWHM (e.g. 2 ppm), solving the Bloch equations independently for each spin packet, and averaging the resulting signals.

Numerical integration (RK4)

The Bloch equations are integrated with the fourth-order Runge–Kutta (RK4) method. Given a time step \(\Delta t\) and current magnetization \(\mathbf{M}^n\):

\[\mathbf{M}^{n+1} = \mathbf{M}^n + \frac{\Delta t}{6} \left(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4\right),\]

where

\[\mathbf{k}_1 = \left.\frac{d\mathbf{M}}{dt}\right|_{\mathbf{M}^n}, \qquad \mathbf{k}_2 = \left.\frac{d\mathbf{M}}{dt}\right|_{\mathbf{M}^n+\frac{1}{2}\mathbf{k}_1\Delta t},\]

\[\mathbf{k}_3 = \left.\frac{d\mathbf{M}}{dt}\right|_{\mathbf{M}^n+\frac{1}{2}\mathbf{k}_2\Delta t}, \qquad \mathbf{k}_4 = \left.\frac{d\mathbf{M}}{dt}\right|_{\mathbf{M}^n+\mathbf{k}_3\Delta t}.\]

The local truncation error is \(\mathcal{O}(\Delta t^5)\) and the accumulated error over \(N = 1/\Delta t\) steps is \(\mathcal{O}(\Delta t^4)\). The time step is chosen to be at least one order of magnitude smaller than all characteristic timescales of the system (relaxation times, axion period). The RK4 kernel is implemented in C++ and exposed to Python via pybind11 for computational efficiency.

Axion power spectral density lineshapes

The analytical PSD of the SHM axion field (Gramolin et al.) is supported for three coupling geometries:

Case	Description
non-grad	Non-gradient coupling
grad_par	Gradient coupling, sensitive axis ∥ \(\mathbf{v}_\mathrm{lab}\)
grad_perp	Gradient coupling, sensitive axis ⊥ \(\mathbf{v}_\mathrm{lab}\)

All lineshapes are one-sided (zero below \(\nu_a\)) and normalized so that \(\int S(\nu)\,d\nu = 1\).

Gravitationally bound axion halo

When axions are gravitationally bound to a compact body (e.g. Earth), they occupy discrete quantum states described by the time-independent Schrödinger equation (TISE):

\[\hat{H}\,\psi = E\,\psi, \qquad \hat{H} = -\frac{\hbar^2}{2m_a}\nabla^2 + m_a\,\Phi(\mathbf{r}),\]

where \(\Phi(\mathbf{r})\) is the gravitational potential. Separating radial and angular parts via \(\psi_{nlm} = R_{nl}(r)\,Y_l^m(\theta,\phi)\), the radial equation for \(u(r) = rR_{nl}(r)\) becomes

\[-\frac{\hbar^2}{2m_a}u'' + V_\mathrm{eff}(r)\,u = E\,u, \qquad V_\mathrm{eff}(r) = m_a\,\Phi(r) + \frac{\hbar^2 l(l+1)}{2m_a r^2}.\]

The package discretises the Hamiltonian on a uniform 1-D radial grid with \(N\) points spanning \(\pm L/2\) using a three-point finite-difference stencil and diagonalizes the resulting dense matrix with scipy.linalg.eigh. Eigenstates are labelled by the spectroscopic convention \(n = n_r + l + 1\) (1s, 2s, 2p, 3s, …).

Earth gravitational potential

EarthBoundAxionHalo pre-configures the solver with Earth’s gravitational potential from the Preliminary Earth Model (PEM) density profile. The interior potential is obtained by integrating the spherical-shell mass distribution; outside the surface it reduces to the point-mass approximation \(\Phi(r) = -GM_\oplus/r\). The reference point is the Earth’s centre (\(\Phi = 0\) at \(r = 0\)).

Axion field gradient at a ground station

For a superposition of eigenstates the total wavefunction is

\[\Psi(\mathbf{r}) = \sum_{n_r, l} c_{n_r l}\,R_{n_r l}(r)\,Y_l^0(\theta,\phi),\]

and its spherical-coordinate gradient is

\[\nabla\Psi = \partial_r\Psi\;\hat{r} + \frac{1}{r}\partial_\theta\Psi\;\hat{\theta} + \frac{1}{r\sin\theta}\partial_\phi\Psi\;\hat{\phi}.\]

The package evaluates each component on a 3-D \((r,\theta,\phi)\) mesh, then interpolates along the radial line pointing toward the experimental station to obtain \(\nabla\Psi(r_\mathrm{station})\).