Prospects for relic neutrino detection using nuclear spin experiments

TL;DR

This paper explores the potential of using nuclear spin experiments, modeled with open quantum systems, to detect the cosmic neutrino background and forecasts sensitivities of future experiments like CASPEr.

Contribution

It extends previous studies by modeling neutrino effects with a Lindblad master equation including experimental imperfections and provides sensitivity forecasts for future quantum sensing experiments.

Findings

CASPEr could constrain C$ u$B overdensity to $ imes 10^{13}$

Optimized setups could reach $ imes 10^{11}$ sensitivity

Quantum sensing offers a promising avenue for fundamental physics

Abstract

Direct detection of the cosmic neutrino background (C$\nu$B) remains one of the most formidable experimental challenges in modern physics. In this work, we extend recent studies of C$\nu$B-induced coherent transitions in polarised nuclear spin ensembles. Adopting an open quantum system framework, we model coherent neutrino effects in large spin ensembles using a Lindblad master equation that also incorporates realistic experimental imperfections such as local dephasing and imperfect polarisation. We solve the Lindblad equation numerically by way of a fast and computationally inexpensive method that can be extended to an arbitrarily large number of spins. Using our numerical solutions, we forecast the sensitivities of future experiments such as CASPEr to the local C$\nu$B overdensity parameter $\delta_\nu$. Our findings indicate that a CASPEr-like experiment, though primarily aimed at axion dark matter search, could also constrain the C$\nu$B overdensity to $\delta_\nu \sim 10^{13}$ in configurations achievable by currently planned experimental efforts, and down to $\delta_\nu \sim 10^{11}$ in the most optimised scenario. While C$\nu$B detection remains out of reach in the foreseeable future, our results highlight the potential of using quantum sensing to probe fundamental physics.

Equations

$$\mathcal{L}$$

$${\mathcal{L}}_{\mathrm{int}}$$

$${\mathcal{L}}_{\mathrm{int}}=-\sum _{ij}\frac{G_F}{\sqrt{2}}\bar{\Psi }{\gamma }^{\mu }\left(g_V^{i,j}-g_A^{i,j}{\gamma }^5\right)\Psi {\bar{\nu }}_i{\gamma }_{\mu }(1-{\gamma }^5){\nu }_j,$$

$$H^{(1/2)}$$

$$={\omega }_0{J}_z+{\Sigma }_z{J}_z+{\Sigma }_{-}{J}_{-}+{\Sigma }_{+}{J}_{+}.$$

$${\Sigma }_z$$

$${\Sigma }_{\pm }$$

$${\gamma }_{\pm }=$$

$$\times \left[(1-N_i(p_{i\pm }^0))N_i(p^0)+(1-{\bar{N}}_i(p_{i\pm }^0)){\bar{N}}_i(p^0)\right],$$

$${\gamma }_{-}^{\mathrm{SM}}\simeq 1.1\times 10^{-48}\mathrm{Hz},\frac{{\gamma }_{+}^{\mathrm{SM}}}{{\gamma }_{-}^{\mathrm{SM}}}\simeq 0.994,$$

$$H=\sum _{\alpha }\left({\omega }_0{J}_z^{\alpha }+{\Sigma }_z^{\alpha }{J}_z^{\alpha }+{\Sigma }_{-}^{\alpha }{J}_{-}^{\alpha }+{\Sigma }_{+}^{\alpha }{J}_{+}^{\alpha }\right),$$

$${\Gamma }_{\pm }=$$

$$\times \left[(1-N_i(p_{i\pm }^0))N_i(p^0)+(1-{\bar{N}}_i(p_{i\pm }^0)){\bar{N}}_i(p^0)\right]$$

$$\times {\int }_{-1}^1\mathrm{d}\mu F_{\pm }({\vec{p}}_{i\pm }-\vec{p}),$$

$$F_{\pm }(\vec{q})\equiv ⟨{\mathrm{i}}_{\mathrm{spin}}|\sum _{\alpha ,\beta }{J}_{\mp }^{\beta }{J}_{\pm }^{\alpha }{e}^{-\mathrm{i}\vec{q}\cdot ({\vec{x}}_{\alpha }-{\vec{x}}_{\beta })}|{\mathrm{i}}_{\mathrm{spin}}⟩$$

$$F_{\pm }^P(\vec{q})$$

$${\Gamma }_{\pm }\approx \frac{N^2}{16p_{\nu }^2{R}^2}{\gamma }_{\pm },$$

$${\Gamma }_{\pm }^{\mathrm{SM}}\simeq 1.2\times 10^{-6}\text{Hz}{\left(\frac{R}{1\text{cm}}\right)}^4,$$

$$H=H_S+H_B+H_{\mathrm{int}},$$

$$\frac{\mathrm{d}{\rho }_S}{\mathrm{d}t}$$

$${\mathcal{D}}_{\mathcal{O}}[{\rho }_S(t)]$$

$$\frac{\mathrm{d}{\rho }_S^{(1/2)}}{\mathrm{d}t}=$$

$$+{\gamma }_{-}(J_{-}{\rho }_S^{(1/2)}{J}_{+}-\frac{1}{2}\{J_{+}{J}_{-},{\rho }_S^{(1/2)}\})$$

$$+{\gamma }_{+}(J_{+}{\rho }_S^{(1/2)}{J}_{-}-\frac{1}{2}\{J_{-}{J}_{+},{\rho }_S^{(1/2)}\})$$

$$+{\gamma }_z(J_z{\rho }_S^{(1/2)}{J}_z-\frac{1}{2}\{J_z^2,{\rho }_S^{(1/2)}\}),$$

$$\frac{\mathrm{d}{\rho }_S}{\mathrm{d}t}=$$

$$+\sum _{\alpha \beta }{\gamma }_{+}^{\alpha \beta }(J_{+}^{\alpha }{\rho }_S{J}_{-}^{\beta }-\frac{1}{2}\{J_{-}^{\alpha }{J}_{+}^{\beta },{\rho }_S\})$$

$$+\sum _{\alpha \beta }{\gamma }_z^{\alpha \beta }(J_z^{\alpha }{\rho }_S{J}_z^{\beta }-\frac{1}{2}\{J_z^{\alpha }{J}_z^{\beta },{\rho }_S\}),$$

$${\gamma }_{\pm }^{\alpha \beta }=$$

$$\times {\delta }^{(1)}({\omega }_0\pm p^{\prime 0}\mp p^0)$$

$$\times \left[(1-N_k(p^{\prime 0}))N_j(p^0)+(1-{\bar{N}}_k(p^{\prime 0})){\bar{N}}_j(p^0)\right].$$

$$\frac{\mathrm{d}{\rho }_S}{\mathrm{d}t}=$$

$$+{\gamma }_{+}({\mathcal{J}}_{+}{\rho }_S{\mathcal{J}}_{-}-\frac{1}{2}\{{\mathcal{J}}_{-}{\mathcal{J}}_{+},{\rho }_S\})$$

$$+{\gamma }_z({\mathcal{J}}_z{\rho }_S{\mathcal{J}}_z-\frac{1}{2}\{{\mathcal{J}}_z^2,{\rho }_S\}),$$

$$\frac{\mathrm{d}{\rho }_S}{\mathrm{d}t}=$$

$$+\sum _{\alpha }{\gamma }_{+}(J_{+}^{\alpha }{\rho }_S{J}_{-}^{\alpha }-\frac{1}{2}\{J_{-}^{\alpha }{J}_{+}^{\alpha },{\rho }_S\})$$

$$+\sum _{\alpha }{\gamma }_z(J_z^{\alpha }{\rho }_S{J}_z^{\alpha }-\frac{1}{2}\{(J_z^{\alpha }{)}^2,{\rho }_S\}).$$

$$n_D=\sum _{j_{\min }}^{N/2}(2j+1{)}^2=\mathcal{O}(N^3)$$

$$|\theta ,\phi {⟩}_{\mathrm{CSS}}$$

Full text

††preprint: CPPC-2025-06

I Introduction

Standard hot big bang cosmology predicts a cosmic neutrino background (C$\nu$B), a relic from the early universe analogous to the cosmic microwave background (CMB), that decoupled from the primordial plasma when it cooled to temperatures of $T\sim 1$ MeV ParticleDataGroup:2024cfk . At decoupling these neutrinos had a nearly thermal spectrum of the relativistic Fermi-Dirac form $n_{F}=[\exp(p/T_{\nu})+1]^{-1}$, which is preserved by cosmological redshifting induced by the spatial expansion of the universe. Together with the minimum neutrino mass values implied by flavour oscillation experiments, we expect the C$\nu$B today to be largely non-relativistic, with a temperature of $T_{\nu,0}\simeq 1.7\times 10^{-4}~{\rm eV}\simeq 1.95$ K, slightly cooler than the $2.73$ K CMB photons. Its mean number density is predicted to be $n_{\nu}=n_{\bar{\nu}}\simeq 56~{\rm cm}^{-3}$ per flavour, but can be enhanced locally, i.e., in the solar neighbourhood, by factors of ${\cal O}(1-10)$, depending on the neutrino mass, through gravitational clustering Singh:2002de ; Ringwald:2004np ; Brandbyge:2010ge ; Villaescusa-Navarro:2011loy ; LoVerde:2013lta ; deSalas:2017wtt ; Zhang:2017ljh ; Mertsch:2019qjv ; Holm:2023rml ; Zimmer:2023jbb ; Zimmer:2024max ; Holm:2024zpr ; Worku:2024kwv .

Nothwithstanding its abundance, the C$\nu$B has so far evaded direct detection in the laboratory owing to its extremely low energies and weak interactions. To date, the only evidence we have of its existence is indirect and based entirely on its gravitational effects in the early universe and subsequent impact on precision cosmological observables. Specifically, measurements of the primordial light element abundances (e.g., Pitrou:2018cgg ) and the CMB anisotropies (e.g., Planck:2018vyg ) consistently prefer an expansion rate—at minutes and 400,000 years post-big bang, respectively—compatible at the 10% level with the standard prediction of the C$\nu$B energy density.111In cosmology parlance, the effect of the C$\nu$B energy density on the early universal expansion rate is parameterised by the effective number of neutrinos $N_{\rm eff}$. The inferred value of $N_{\rm eff}$ from primordial light element abundances and the CMB are $N_{\rm eff}=2.88\pm 0.27$ (68% C.I.) Pitrou:2018cgg and $N_{\rm eff}=2.99\pm 0.17$ (68% C.I.) Planck:2018vyg , respectively, consistent with the energy density in three light neutrino flavours. Anisotropic stress constraints from the CMB likewise support the weakly-interacting nature of the C$\nu$B Oldengott:2017fhy ; Chen:2022idm .

In the face of ever-improving cosmological/gravity-based constraints, it nonetheless remains of great interest to detect or at least constrain the C$\nu$B in a more direct way, using methods based on particle scattering. In this connection, a number of possible probes and constraints have been discussed. Below is a non-exhaustive list.

• •

At leading order in the Fermi constant $G_{F}$, the so-called Stodolsky effect Stodolsky:1974aq refers to a spin-dependent shift in the energy of a standard-model (SM) fermion sitting in a $CP$ asymmetric bath of neutrinos and antineutrinos. While cosmological neutrino-antineutrino asymmetries of $\eta_{\nu}\equiv(n_{\nu}-n_{\bar{\nu}})/n_{\gamma}\lesssim{\cal O}(10^{-2})$ are not excluded by observations (e.g., Oldengott:2017tzj ), it is also not well motivated within the SM, where ${\cal O}(10^{-10})$ particle-antiparticle asymmetries are generically expected in light of the observed matter-antimatter asymmetry in the universe. We remark however on a recent claim of a large $\eta_{\nu}\sim 10^{-2}-0.25$ being preferred by extragalactic 4He observations Burns:2022hkq . Note also discussions of the survival of such asymmetries in the local universe Arvanitaki:2022oby ; Ruchayskiy:2022eog .

• •

At next-to-leading order ${\cal O}(G_{F}^{2})$, coherent elastic scattering of a C$\nu$B wind with a macroscopic target can produce a mechanical force that causes the target to accelerate Shvartsman:1982sn ; Langacker:1982ih ; Smith:1983jj ; Duda:2001hd ; Domcke:2017aqj ; Shergold:2021evs . This force does not require the presence of a neutrino-antineutrino asymmetry, but works best for Dirac neutrinos. For Majorana neutrinos, the absence of vector currents causes the interaction cross section to be suppressed by a factor $v^{2}\sim 10^{-6}$, where $v$ is the Earth-C$\nu$B relative speed. Current torsion balance setups cannot measure this acceleration unless the C$\nu$B number density is enhanced by a factor $\delta_{\nu}\sim 10^{12}$ or more Shergold:2021evs . Advancement on this front consists in using new techniques, e.g., interferometry Domcke:2017aqj , to enhance the sensitivity to the acceleration.

• •

Another ${\cal O}(G_{F}^{2})$ effect, neutrino capture by $\beta$-decaying nuclei remains unique to C$\nu$B detection Weinberg:1962zza ; Cocco:2007za ; Long:2014zva ; Akita:2020jbo and goes hand-in-hand with direct neutrino mass searches based upon $\beta$-decay end-point spectrum measurements KATRIN:2022kkv ; KATRIN:2024cdt ; Project8:2022wqh . The flagship tritium $\beta$-decay neutrino mass experiment KATRIN currently limits the local C$\nu$B overdensity to $\delta_{\nu}\lesssim 10^{11}$ KATRIN:2022kkv . The proposed PTOLEMY experiment aims to constrain this parameter to $\delta_{\nu}\sim{\cal O}(1)$ PTOLEMY:2018jst .

• •

An accelerator-based variant of the above has also been discussed in Ref. Bauer:2021uyj , which proposes to use accelerated beams of ions in ion storage rings to capture C$\nu$B neutrinos at nuclear $\beta$-resonances to enhance the capture cross section.

• •

Inelastic ${\cal O}(G_{F}^{2})$ scattering of the C$\nu$B on a target system can also induce transitions between different energy states. These include atomic transitions Yoshimura:2014hfa ; Huang:2025yqu , as well as nuclear spin flips Bauer:2025hdq ; Arvanitaki:2024taq . Of the latter, Hydrogen spin flip induced by absorption of a relic neutrino-antineutrino pair may be sensitive to $\delta_{\nu}\sim 10^{9}$ for neutrino masses $m_{\nu}\sim 10^{-3}$ eV Bauer:2025hdq .

• •

On astrophysical scales, probes of the C$\nu$B range from small-scale (i.e., sub-1 AU) perturbations to planetary orbits (current constraint $\delta_{\nu}\lesssim 5.34\times 10^{10}$ Aliberti:2024udm ) and the cooling of old neutron stars from C$\nu$B capture (potential sensitivity $\delta_{\nu}\sim 10^{9}$ Das:2024thc ), to large-scale ($\gtrsim{\cal O}(10)$ Mpc) up-scattering of the C$\nu$B by cosmic rays Ciscar-Monsalvatje:2024tvm ; DeMarchi:2024zer ; Herrera:2024upj (current constraint $\delta_{\nu}\lesssim 10^{4}$ Herrera:2024upj ) and features in the energy spectra of ultra-high energy neutrinos due to C$\nu$B scattering (current constraint $\delta_{\nu}\lesssim 10^{14}$ Franklin:2024amy ). Also proposed are resonant absorption dips in the cosmogenic neutrino energy spectrum when these neutrinos interact with the C$\nu$B on a meson resonance Brdar:2022kpu , a variant of the original $Z$-resonance absorption dips for ultra-high energy neutrinos Weiler:1982qy .

We are particularly interested in the proposal of Ref. Arvanitaki:2024taq , which exploits spin ensembles such as in nuclear magnetic resonance (NMR) experiments to detect coherent effects induced by the C$\nu$B through inelastic neutrino-spin interactions. Specifically, Ref. Arvanitaki:2024taq argued that the momentum transfer $q$ in such interactions can be macroscopic in wavelength terms, causing the collective C$\nu$B-induced spin flip rates across the spin ensemble to scale as $N^{2}$, where $N$ is the number of spins, if the spin sample’s physical size $R$ is smaller than or comparable to $q^{-1}$ (this situation is to be compared with an $\propto N$ scaling in an incoherent, $qR\gg 1$ setting). Such macroscopic coherence of the target system can significantly enhance sensitivity to weak interactions and may conceivably yield a constraint on the C$\nu$B overdensity parameter comparable to the current KATRIN bound, $\delta_{\nu}\lesssim 10^{11}$ KATRIN:2022kkv , for specific initial spin configurations Arvanitaki:2024taq . Added to its appeal is that planned NMR experiments to search for dark matter axions and axion-like particles, e.g., CASPEr Budker:2013hfa ; JacksonKimball:2017elr , are in principle sensitive to such coherent signals. Besides opening up a new avenue for ultra-low-threshold searches for fundamental physics, the potential for these experiments to simultaneously search for and/or constrain both the axion and C$\nu$B parameter space is a most tantalising prospect.

In this work, we expand on the study of Ref. Arvanitaki:2024taq as follows: (i) We explore realistic experimental effects such as local dephasing and imperfect initial polarisation of the spin ensemble, both of which can degrade the coherent C$\nu$B signature. (ii) Modelling these effects using an open quantum system framework, we solve the Lindblad master equation numerically, both for the full density matrix of $N$ spins, and using a computationally efficient second-order approximation that tracks only a limited set of observables. Using these solutions, we are able to confirm for large $N$ the validity of the perturbative solutions put forward in Ref. Arvanitaki:2024taq . (iii) We forecast the sensitivity of realistic future experiments to the C$\nu$B overdensity parameter $\delta_{\nu}$, especially the reach of the CASPEr experiment Budker:2013hfa ; JacksonKimball:2017elr ; Walter:2025ldb , by including a more realistic modelling of the instrumental noise, initial polarisation of the spin ensemble, and sample size.

In addition to the sample size, which restricts the degree of C$\nu$B coherence across the spin ensemble, we find the initial spin polarisation to be a particularly important limiting factor—and also the primary reason why the current generation of NMR axion searches cannot place competitive constraints on $\delta_{\nu}$. In the most optimistic case (no instrumental noise, $R\sim 10$ cm, 100% polarisation), a sensitivity of $\delta_{\nu}\sim 10^{11}$ may be achievable using 129Xe nuclei. In more realistic settings (instrumental noise, $R\sim 1$ cm, 25% polarisation), however, the sensitivity to the C$\nu$B deteriorates to $\delta_{\nu}\sim 10^{13}$.

The paper is structured as follows. In Sec. II we introduce the scattering rates of the C$\nu$B with a single two-level system, such as a nuclear spin with $I=1/2$, and with an ensemble of $N$ spins, highlighting a regime where a coherent $\sim N^{2}$ enhancement can be achieved. In Sec. III, we review the Lindblad master equation commonly used to model open quantum systems, and describe the dynamics in the Dicke basis. Sections IV and V discuss, respectively, the emergence of $N^{2}$ effects in collective observables and the limitations arising from incoherent local effects that eventually degrade the signal. We show in Sec. VI how the Lindblad master equation combining coherent and incoherent local effects can be efficiently solved using a numerical approximation. We use these numerical solutions to forecast the sensitivity of currently-planned experiments to the C$\nu$B overdensity parameter $\delta_{\nu}$, as detailed in Sec. VII. Section VIII contains our conclusions. The computational details of this work are reported in eight appendices.

II Neutrino-spin interaction

Consider a spin-$1/2$ fermion—such as an electron or a nucleus—interacting with neutrinos via the weak interaction. We assume coherence over nuclear distances, allowing us to treat the fermion as a point-like particle and neglect its internal structure. At low energies, this interaction can be described by an effective 4-fermi interaction specified by the Lagrangian

[TABLE]

where $\Psi$ denotes the fermion field of mass $m$ and electric charge $q$, $A_{\mu}=(A_{0},\vec{A})$ is the electromagnetic 4-potential, $G_{F}$ the Fermi constant, $\nu_{f}$ the neutrino field of flavour $f=e,\mu,\tau$, and $g_{V}^{f}$ and $g_{A}^{f}$ are the vector and axial-vector couplings of the fermion to $\nu_{f}$.

The couplings $g_{V}^{f}$ and $g_{A}^{f}$ generally depend on the neutrino flavour as well as the nature of the fermion. Here, we focus exclusively on neutrino–nucleus interaction through the neutral-current channel. The corresponding couplings are therefore independent of flavour. The values of $g_{V}$ and $g_{A}$ are nucleus-dependent. In what follows, we consider the axial coupling $g_{A}$, which is only sensitive to the spin of the unpaired nucleon (see, e.g., bohr1998nuclear ).

It is convenient to express the neutrino fields in the mass basis. Performing the basis transformation $\nu_{f}=\sum_{i}U_{fi}\nu_{i}$, where $\nu_{i}$ are the mass eigenstates (for $i=1,2,3$) and $U_{fi}$ is the Pontecorvo-Maki-Nakagawa-Sakata matrix, the interaction Lagrangian becomes

[TABLE]

where the effective couplings in the mass basis are given by $g^{i,j}_{V(A)}=\sum_{f}g_{V(A)}^{f}U_{fi}^{*}U_{fj}$. In the flavour-independent case under consideration, the effective couplings are diagonal, i.e., non-zero only for $i=j$, and equal for all mass eigenstates.222If instead of nuclear spins we consider electron spins, then the neutrino-spin interaction can also proceed via the charged-current channel. In that case, the couplings $g_{V(A)}^{i,j}$ generally contain non-zero off-diagonal terms, allowing the neutrino to change mass eigenstates via the interaction. This case has been considered in Ref. Arvanitaki:2024taq .

In a typical NMR setup, a sample of (non-relativistic) nuclear spins is placed in a constant magnetic bias field $\vec{B}=\nabla\times\vec{A}$, which we choose to point in the negative-$z$ direction, i.e., $\vec{B}=-B\hat{z}$. We can therefore take the non-relativistic limit and derive from the relativistic Lagrangian (1) the corresponding Hamiltonian for a single spin,

[TABLE]

Here, $H_{S}^{(1/2)}=\omega_{0}J_{z}$ describes the interaction of the spin with the magnetic field, where $\omega_{0}=\gamma B$ is the energy splitting between adjacent spin states—also called the Larmor frequency, and $\gamma=g\,\mu_{N}$ is the gyromagnetic ratio specific to the nucleus under consideration, defined via the nuclear magneton $\mu_{N}=e/2m_{p}$ and the nucleus-specific $g$-factor. The $H_{\rm int}^{(1/2)}$ terms describe interaction of the spin with the C$\nu$B, where the spin raising and lowering operators are defined as $J_{\pm}\equiv J_{x}\pm{\rm i}J_{y}$. In our flavour-independent case the neutrino terms $\Sigma_{\pm,z}$ read

[TABLE]

and they arise exclusively from the axial part of the weak interaction. We refer to the reader to Appendix A for the derivation of the low-energy Hamiltonian, including the flavour-dependent case.

Physically, the term $\Sigma_{z}J_{z}$ corresponds to elastic scattering of neutrinos off a spin, which alters the phase evolution of the spin but does not change its energy. In contrast, the $\Sigma_{\pm}J_{\pm}$ terms describe inelastic processes, in which the neutrino exchanges energy with the spin, either exciting it ($\Sigma_{+}$) or de-exciting it ($\Sigma_{-}$). The single-spin excitation and de-excitation rates $\gamma_{\pm}$, defined as the transition rate of a spin from down to up ($+$) and vice versa ($-$), can be estimated using Fermi’s golden rule. Accounting only for inelastic neutrino-spin scattering, the rates read

[TABLE]

where $p_{i\pm}=[(\sqrt{|\vec{p}|^{2}+m_{i}^{2}}\mp\omega_{0})^{2}-m_{i}^{2}]^{1/2}$, $p^{0}_{i\pm}=\sqrt{|\vec{p}|^{2}+m_{i}^{2}}\mp\omega_{0}$, $m_{i}$ is the mass of the $i$th mass eigenstate, and the lower integration limits are given by $p^{\rm low}_{i+}=\sqrt{(\omega_{0}+m_{i})^{2}-m_{i}^{2}}$ and $p^{\rm low}_{i-}=0$.333All neutrino momenta throughout this work are 3-momenta. We do however use a superscript “0”, e.g., $p^{0}\equiv\sqrt{|\vec{p}|^{2}+m^{2}}$ to denote energy. This expression applies to both Dirac and Majorana neutrinos, but with the understanding that the neutrino and antineutrino occupation numbers, $N_{i}$ and $\bar{N}_{i}$, in the Dirac case are a factor two smaller than in the Majorana case.444Although at decoupling relativistic Dirac neutrinos (antineutrinos) are exclusively in the left- (right-) helicity state, as they redshift to non-relativistic velocities at late times, the generic expectation is they should relax into an equal mixture of left- and right-helicity states via gravitational interactions. The detailed calculation of $\gamma_{\pm}$ from the low-energy Hamiltonian (3) can be found in Appendix B. Assuming a standard, $CP$-symmetric C$\nu$B with $N_{i}=\bar{N}_{i}\sim n_{F}$, and a spin sample consisting of ${}^{129}\mathrm{Xe}$ nuclei with $g_{A}=1/2$, the rates (5) evaluate numerically to

[TABLE]

for a typical NMR energy splitting of $\omega_{0}=1\times 10^{-8}$ eV and a normally-ordered neutrino mass sum of $\sum m_{\nu}=0.15$ eV.555We adopt mass values consistent with global neutrino oscillation fits Esteban:2020cvm ; deSalas:2020pgw ; Capozzi:2017ipn . For normal ordering, the masses are fixed by the total mass sum and the measured mass-squared splittings, $\Delta m_{21}^{2}\simeq 7.4\times 10^{-5}~\mathrm{eV}^{2}$ and $\Delta m_{31}^{2}\simeq 2.5\times 10^{-3}~\mathrm{eV}^{2}$, giving $m_{1}\simeq 0.0419~\mathrm{eV}$, $m_{2}\simeq 0.0428~\mathrm{eV}$, and $m_{3}\simeq 0.0652~\mathrm{eV}$ for $\sum m_{\nu}=0.15$ eV.

In the case of $N$ spins, the low-energy Hamiltonian (3) generalises to

[TABLE]

where $J_{\pm,z}^{\alpha}$ are understood to operate only on the spin at spatial position $\vec{x}_{\alpha}$, and the neutrino terms $\Sigma^{\alpha}_{\pm,z}$ are the same as those given in Eq. (4) but now with the neutrino fields specified by $\nu_{i}(t,\vec{x}_{\alpha})$. The total excitation and de-excitation rates of the ensemble, defined here as the rate at which the total energy of the ensemble increases or decreases by one unit of $\omega_{0}$, can be written as666Throughout this work, we always denote integrated (ensemble-level) rates by an upper-case $\Gamma$, while single-spin interaction rates—relevant for the Lindblad dynamics discussed in Sec. III—are represented by a lower-case $\gamma$.

[TABLE]

where $\mu\equiv\hat{p}_{i\pm}\cdot\hat{p}$, is the polar angle between the momentum vectors $\vec{p}$ and $\vec{p}\,^{\prime}$, and

[TABLE]

are the excitation and de-excitation form factors.

Depending on the initial spin configuration and the momentum transfer $\vec{q}\equiv\vec{p}_{i\pm}-\vec{p}$ from the neutrinos to the ensemble, the form factors $F_{\pm}(\vec{q})$ can evaluate to a range of values. For example, a system prepared in the ground state $|G\rangle\equiv|\downarrow\rangle^{\otimes N}$ in which all spins align initially with the magnetic field has $F^{G}_{+}(\vec{q})=N$ and $F^{G}_{-}(\vec{q})=0$, and hence a non-vanishing total excitation rate $\Gamma_{+}=\gamma_{+}N$. In contrast, for an initial spin state aligned with the equatorial plane, $|P\rangle=\prod_{\alpha}\frac{1}{\sqrt{2}}(|\uparrow\rangle+|\downarrow\rangle)$, the form factors evaluate to

[TABLE]

with $j_{1}(x)$ the spherical Bessel function of order one. Clearly, $F^{P}_{\pm}(\vec{q})$ can range from ${\cal O}(N)$ in the fully incoherent regime, to ${\cal O}(N^{2})$ in the fully coherent regime, discussed in the following subsection.

II.1 Coherent regime

An inspection of the form factor (10) reveals that there are essentially two physical scales that determine the spatial coherence of the neutrino-spin interactions. The first is the characteristic size of the spin ensemble $R$, which sets the spatial extent over which the spins are distributed. The second is the typical momentum transfer induced by the relic neutrinos $q$.

Since the incoming neutrinos follow a momentum distribution $N_{i}(p^{0})$, the momentum transfer $\vec{q}=\vec{p}_{i\pm}(\vec{p})-\vec{p}$ is also dictated by the distribution. The coherence properties must therefore be understood in terms of the typical momentum transfer ${}\vec{q}_{\rm typ}$ extracted from this distribution. When the condition $q_{\rm typ}R\ll 1$ is satisfied, we see in Eq. (10) that $j_{1}(q_{\rm typ}R)\approx q_{\rm typ}R/3$, leading to $F_{\pm}\sim{\cal O}(N^{2})$. This condition can always be fulfilled if the typical neutrino momentum $p_{\nu}$ itself is small compared with the inverse system size, i.e., $p_{\nu}R\ll 1$. For relic neutrinos with a typical incoming momentum of $p_{\nu}\sim 3.15\,T_{\nu,0}\simeq 5.3\times 10^{-4}~\mathrm{eV}\sim 27~\mathrm{cm}^{-1}$, the condition $p_{\nu}R\ll 1$ therefore requires a very small $R\ll 0.04~\mathrm{cm}$. However, even for the larger-size spin samples to be considered in this work, e.g., $R\sim{\cal O}(1)$ cm, where formally $p_{\nu}R\gg 1$, it is still possible to fulfil the broader condition $q_{\rm typ}R\ll 1$ and hence maintain coherence across the sample in some limits. Specifically, the typical C$\nu$B kinetic energy is much larger than the energy splittings $\omega_{0}\sim 10^{-8}$ eV considered in this work. In this limit, the typical momentum transfer is $q_{\rm typ}\approx 2p_{\nu}\sin(\theta/2)$, where $\theta$ is the angle between the incoming and outgoing neutrino. Thus, we see that coherence can be maintained for those scattering angles satisfying $\theta\ll(p_{\nu}R)^{-1}$.

In terms of the total excitation and de-excitation rates (8), the $p_{\nu}R\gg 1$ limit is effectively equivalent to demanding that only outgoing momentum directions falling in the range $\cos((p_{\nu}R)^{-1})<\mu<1$ contribute to the $N^{2}$ enhancement. Thus, evaluating the $\mu$-integral immediately yields

[TABLE]

which we note is suppressed by a factor $\sim(2p_{\nu}R)^{-2}$ relative to the naïve coherent rate. The typical momentum transfer in this regime is $q_{\rm typ}\sim\omega_{0}m_{\nu}/p_{\nu}$. This will be the case of interest in our study and was first introduced in Ref. Arvanitaki:2024taq as an analogy to the Rayleigh-Gans regime 1925AnP…381…29G . See Appendix B for details.

As a concrete example, consider a spin ensemble composed of liquid 129Xe characterised with a spin density of $n_{s}\sim N/R^{3}\sim 10^{22}~\text{cm}^{-3}$. In this case, assuming an energy splitting $\omega_{0}=10^{-8}~\text{eV}$ and a standard C$\nu$B with three normally-ordered masses summing to $\sum m_{\nu}=0.15~\text{eV}$, the total excitation rate of the spin ensemble evaluates to

[TABLE]

where the $R^{4}$ scaling represents the aforementioned forward limit of the $N^{2}$ enhancement Arvanitaki:2024taq . The dependence of the total interaction rate ($\Gamma_{\rm tot}=\Gamma_{+}+\Gamma_{-}$) on the system parameters, including $\omega_{0}$ and the sample size $R$, is illustrated in the left panel of Fig. 1, where regions exhibiting enhanced coherence can be identified. The right panel of the same figure shows the dissipation ratio $\Gamma_{+}/\Gamma_{-}\equiv\gamma_{+}/\gamma_{-}$ across the relevant parameter space, and highlights a slight dominance of neutrino-induced de-excitation. This dominance arises because de-excitation is always kinematically allowed as it has no energy threshold: any neutrino regardless of its momentum can induce the transition; see Eq. (5). Excitation is on the other hand not always possible. The fractional difference is approximately $\gamma_{+}/\gamma_{-}-1\sim-2m_{\nu}\omega_{0}/p_{\nu}^{2}$ Arvanitaki:2024taq for $\omega_{0}\ll p_{\nu}^{2}/(2m_{\nu})\ll m_{\nu}$. See also Appendix B.

III An open quantum system

Although our full system comprises a spin ensemble $S$ and a relic neutrino bath $B$ that interact with each other, in practice we are interested only in the dynamics of the former. To this end, we adopt the framework of open quantum systems, which assumes unitary evolution of the density matrix of the full system, $\rho=\rho_{S}\otimes\rho_{B}$, under a total Hamiltonian

[TABLE]

where $H_{S}$ and $H_{B}$ denote, respectively, the free Hamiltonian of the spin ensemble and of the bath, and $H_{\rm int}$ describes their interaction. Tracing out the bath degrees of freedom, the evolution of the reduced density matrix $\rho_{S}$ is described by the Lindblad master equation Lindblad:1975ef ; Gorini:1975nb ; Manzano:2020yyw ,

[TABLE]

where the so-called Lamb-shift Hamiltonian $H_{LS}$ is a small hermitian correction to $H_{S}$ arising from $H_{\rm int}$, and

[TABLE]

is a super-operator that acts on the reduced density matrix and describes the dissipative effects of the system-bath interaction. Dissipation can occur through multiple channels, each specified by its own operator ${\cal O}_{k}$ and the associated interaction rate $\gamma_{k}$.

The Lindblad equation (14) serves well under the condition of weak coupling of the system to a large bath with no memory and when the Larmor frequency $\omega_{0}$ is much larger than the interaction rates. These conditions enter the derivation of Eq. (14) at various stages via the so-called Born, Markovian, and rotating-wave/secular approximation, which are discussed in more detail in Appendix C. See also, e.g., Ref. Manzano:2020yyw for a general introduction to the Lindblad formalism.

III.1 The master equation

In the case of a single spin-$1/2$ interacting with a neutrino bath, the Hamiltonians of the spin system and its interaction with the neutrinos, $H_{S}^{(1/2)}$ and $H_{I}^{(1/2)}$, are given in Eq. (3). The corresponding Lindblad equation for the spin density matrix in the interaction picture takes the form

[TABLE]

where we have highlighted $\rho^{(1/2)}_{S}$ as the reduced density matrix for a single two-level system. We refer the reader to Appendix C for an ab initio derivation of this equation.

Inspecting Eq. (16), we see firstly an energy shift due to the $\Delta$, which is a combination of the ${\cal O}(G_{F})$, Stodolsky effect that is present only when the neutrino-antineutrino asymmetry is non-zero, and ${\cal O}(G^{2}_{F})$ terms from the Hermitian part of the Lindblad equation. The $J_{\pm}$ terms and the corresponding coefficients $\gamma_{\pm}$ represent the non-Hermitian interaction of the spin with the neutrinos, causing the spin to either excite or de-excite via inelastic scattering, analogously to absorption and emission in quantum optics.777While we draw parallels between our system and quantum optics and introduced terminology commonly used in the latter subject, we always use “excitation” and “de-excitation” to refer to the effects of the C$\nu$B. The $J_{z}$/$\gamma_{z}$ term is, on the other hand, energy conserving and represents the effect of dephasing. The computations of the coefficients $\gamma_{\pm,z}$ can be found in Appendix D. Their general forms including both neutrino scattering and pair emission/absorption processes can be found in Eqs. (133), (134), and (135), and, for $\gamma_{\pm}$, match the heuristic single-spin excitation and de-excitation rates (5) under the same settings.

The single-spin equation (16) can be generalised to $N$ spins, where the Hilbert space of the spin ensemble is now $\rho_{S}=\otimes^{N}\rho^{(1/2)}$ and Eq. (7) gives the corresponding Hamiltonian. As we show in Appendix E, omitting the Hermitian part, the Lindblad equation for the $N$ two-level systems reads

[TABLE]

where $\alpha$ and $\beta$ run over all the spins of the system, and the angular operators $J_{\pm,z}^{\alpha}$ act only on the spin $\alpha$. The coefficients $\gamma^{\alpha\beta}_{\pm,z}$ again encode information of the spin-bath interaction, but now correlate the interaction at spatial position $\vec{x}_{\alpha}$ to that at $\vec{x}_{\beta}$. For instance, assuming only inelastic neutrino-spin scattering, the $\gamma_{\pm}^{\alpha\beta}$ coefficients read

[TABLE]

Depending on the extent to which correlations can be maintained across spatial distances $\Delta\vec{x}_{\alpha\beta}\equiv\vec{x}_{\alpha}-\vec{x}_{\beta}$ of the spin ensemble, qualitatively very different dynamics can arise. Two contrasting limits are the fully-coherent and the fully-incoherent regimes.

Fully-coherent regime

Full coherence is achieved when $e^{-{\rm i}\Delta\vec{x}_{\alpha\beta}.(\vec{p}-\vec{p}\,^{\prime})}\approx 1$ across the whole spin ensemble, or, equivalently, when the condition $qR\ll 1$ is satisfied for momentum transfer $q$ and a sample size $R$. In this limit, the interaction coefficients are effectively the same for all spin pairs. The Lindblad equation then reduces to

[TABLE]

where we have introduced the collective spin operators ${\cal J}_{\pm,z}=\sum_{\alpha}J_{\pm,z}^{\alpha}$. This collective Lindblad equation is structurally identical to that for a single spin and the coefficients $\gamma_{\pm,z}$ are the same single-spin interaction rates of Eqs. (133), (134), and (135). In the language of quantum optics, these terms describe, in the order of appearance, collective emission, collective absorption, and collective dephasing (see, e.g., Ref. PhysRevA.98.063815 ). Since the system Hamiltonian is $H_{S}=\omega_{0}{\cal J}_{z}$, the collective dephasing term commutes with $H_{S}$ and has no dynamical effect on the observables of interest; we will therefore neglect it in the following.

Incoherent regime

In the opposite, $qR\gg 1$ limit, the phase differences across the system are large, leading to a complete loss of coherence. This defines the incoherent regime, in which neutrino-spin interactions are fully local. The master equation in this case factorises into $N$ independent Lindblad equations, i.e.,

[TABLE]

Here, each spin evolves independently under local interactions with the bath, and the rates $\gamma_{\pm,z}$ are again the single-spin rates. This local Lindblad form is general for any uncorrelated ensemble-bath interaction and applies to a wide class of physical systems for which the Lindblad approximation holds. In Sec. V, we will use Eq. (20) to model other (non-neutrino) environmental effects that eventually drive the spin ensemble towards thermal equilibrium in a realistic NMR experiment.

III.2 Dicke basis and initial states

Suppose all $N$ spins are identically prepared and indistinguishable, and experience the same interactions with the environment. It is then convenient to use the framework of Dicke states PhysRev.93.99 , which are eigenstates of the collective operators ${\cal J}^{2}$ and ${\cal J}_{z}$, and form a basis for the permutationally invariant subspace of the full Hilbert space. As permutational invariance is preserved by the dynamics of the scenarios of interest, adoption of the Dicke basis allows us to drastically reduce the computational complexity of the problem, both analytically and numerically: instead of working in the full $2^{N}$-dimensional Hilbert space, we can restrict the dynamics to a much smaller subspace, spanned by the Dicke states labelled by the total spin and its projection, $|j,m\rangle$, where $j=N/2,N/2-1,\ldots,j_{\min}$, with $j_{\min}=0$ if $N$ even or $j_{\min}=1/2$ if $N$ odd, and $|m|\leq j$. The full density matrix therefore contains

[TABLE]

elements. In the presence of purely-coherent interactions, where the dynamics are confined to the fully symmetric subspace ($j=N/2$), this reduces further to $n_{D}=N+1$.

The left panel of Fig. 2 illustrates the intuition provided by the Dicke representation—commonly visualised as the Dicke triangle—which organises the Hilbert space according to total spin $j$, including both the symmetric subspace ($j=N/2$) and the non-symmetric sectors ($j<N/2$) that do not exhibit cooperative or superradiant behaviours. When the coherence condition is satisfied, the neutrino-spin interaction is represented by collective, “vertical” transitions along the Dicke ladder in steps of $|N/2,m\rangle\to|N/2,m\pm 1\rangle$. This behavior contrasts with that of incoherent local processes such as dephasing, which induce only “horizontal” motion along the $j$-axis of the Dicke triangle. For more details on the Dicke basis and the effects of local interactions, we refer the reader to Ref. PhysRevA.98.063815 and references therein. See also Ref. Gross:1982dkt ; Masson:2021eyk for an exposition on Dicke superradiance.

We are interested in three types of permutation-symmetric initial states, whose density matrices in the symmetric Dicke basis and Bloch sphere visualisations are shown in the centre and right panels of Fig. 2.

• •

The ground state $|G\rangle\equiv|\downarrow\rangle^{\otimes N}=|N/2,-N/2\rangle$ represents the case in which all $N$ spins are in the low-energy state of the splitting (pointing down in our notation), such that the polarisation is unity in the Bloch sphere. From an experimental perspective, such a state is called hyperpolarised, as the polarisation far exceeds the thermal equilibrium expectation (see Sec. V). Nonetheless, hyperpolarisation can be achieved experimentally for, e.g., 129Xe-enriched ensembles, via spin exchange optical pumping, as planned in the upcoming stages of the CASPEr experiment Walter:2025ldb .

• •

The fully excited state $|E\rangle\equiv|\uparrow\rangle^{\otimes N}=|N/2,N/2\rangle$ has all $N$ spins in the high-energy state of the splitting (pointing up). Starting from a hyperpolarised state like $|G\rangle$, the $|E\rangle$ state can be produced by applying a $\pi$ pulse, i.e., applying a weak magnetic field oscillating at the Larmor frequency, until the polarisation is inverted.

• •

The equatorial coherent spin state (CSS), or “product state”, is defined by $|P\rangle\equiv|\theta=\pi/2,\varphi=0\rangle_{\rm CSS}$, where

[TABLE]

In the Dicke basis it can be written as $|P\rangle=\sum_{m}c_{m}|j,m\rangle$, with $c_{m}=2^{-j}\binom{2j}{j+m}$. This state serves as the standard benchmark in quantum sensing because of its resemblance to the superradiant state $|S\rangle\equiv|N/2,0\rangle$ PhysRev.93.99 . It can be prepared by applying a $\pi/2$ pulse to the ground state.

The preparation of a hyperpolarised ground state $|G\rangle$ is critical for all three initial configurations. However, achieving a pure state with near-unity polarisation approximating $|G\rangle$ may not always be possible under realistic experimental settings. A more faithful description of the spin ensemble’s initial configuration should therefore be given in terms of a density matrix, shown here in the Dicke basis in a generic form,888The complete Dicke basis should in principle include a degeneracy label for each $|j,m\rangle$ state, i.e., $|j,m,\alpha\rangle$, with $1\leq\alpha\leq d_{j}^{N}$, where $d_{j}^{N}$ can be found in, e.g., Eq. (5) of Ref. Shammah_2017 , in order to recover the full Hilbert space. For simplicity, however, we have omitted this label, since it does not affect our arguments.

[TABLE]

The initial polarisation of the spin ensemble prior to the application of any pulse can then be trivially specified via the expectation value of the collective operator ${\cal J}_{z}$, i.e.,

[TABLE]

which always evaluates to $|p|<1$ if the configuration is not purely the ground ($m=-N/2$) or purely the excited ($m=N/2$) state.

We caution however that the definition (24) of $p$ alone does not completely specify the effect of partial polarisation on the spin ensemble. To see this, consider also the expectation values

[TABLE]

where the upper bounds correspond to the expectations for the ground ($j=N/2,m=-N/2$) and the excited ($j=N/2,m=N/2$) states. Clearly, for $\langle{\cal J}_{z}^{2}\rangle$ to be suppressed relative to the optimal value we must have admixtures of $|m|<N/2$ states (i.e., states along the Dicke ladder in Fig. 2), while for $\langle{\cal J}^{2}\rangle$ to become suboptimal the ensemble must contain $j<N/2$ states (i.e., states along the $j$-axis in the Dicke triangle). Crucially, neither $\langle{\cal J}_{z}^{2}\rangle$ nor $\langle{\cal J}^{2}\rangle$ can be expressed in general in terms of $p$ as defined in Eq. (24).

We shall return to the issue of polarisation in Sec. VI.1, where we numerically evaluate in the Lindblad equation and connect our solutions with the hyperpolarisation protocols and level of polarisation achievable in current and future experimental setups.

IV Coherent observables

Given an initial state, the Lindblad equation (17) can be solved to find the evolution of the spin ensemble. Several observables are of interest. The most common first-order observable is the collective magnetisation, whose expectation value can be constructed from $\langle{\cal J}_{z}\rangle={\rm tr}(\rho_{S}{\cal J}_{z})$, and is related to the ensemble’s energy. At second order, relevant observables include $\langle{\cal J}^{2}_{z}\rangle$ and $\langle{\cal J}^{2}_{x}\rangle$, where the latter represents the squared transverse spin component.

Exact analytical solutions for these observables generally do not exist. Numerical solutions, on the other hand, can become prohibitive for large values of $N$. We therefore consider first approximate analytical solutions in different limits and under different settings in order to map out the limiting behaviours of the system, before attempting a numerical treatment. We discuss first in this section approximate analytical solutions in the presence of coherent neutrino effects only. Section V considers the effects of non-neutrino “noise” from local interactions. We present numerical solutions including both neutrino-induced coherent and local noise effects in Sec. VI.

Before we delve into the solutions, we first define two timescales that characterise the evolution of the spin ensemble in the presence of coherent neutrino interactions.

• •

The superradiant timescale,

[TABLE]

characterises the timescale for a single collective spin-flip, i.e., the transition $|j,m\rangle\to|j,m\pm 1\rangle$. The numerical estimate follows from setting $\gamma_{\pm}=\gamma_{\pm}^{\rm SM}$, where $\gamma_{\pm}^{\rm SM}$ is given in Eq. (6).

• •

The dissipative timescale,

[TABLE]

is the dynamical timescale of the system, defined as the duration over which an ${\cal O}(1)$ change in a collective observable such as $\langle{\cal J}_{z}\rangle$ becomes apparent. This timescale also marks the end of the coherent superradiant regime, and limits the validity of the perturbative solution discussed below. The numerical estimate assumes again $\gamma_{\pm}=\gamma_{\pm}^{\rm SM}$ from Eq. (6).

IV.1 Perturbative solution

Considering the weakness of the neutrino-spin interactions, Ref. Arvanitaki:2024taq proposed a perturbative/iterative solution to Eq. (19) around the initial state. Assuming only coherent neutrino interactions are present and to linear order in $t$, the solution reads

[TABLE]

recalling the notation of ${\cal D}_{\cal O}[\rho]$ from Eq. (15), and we have again omitted the collective dephasing term.

In this perturbative limit, the evolution of $\langle{\cal J}_{z}\rangle$ for an ensemble initially prepared in the $|P\rangle$ state is given by

[TABLE]

with $\gamma_{\rm net}\equiv\gamma_{-}-\gamma_{+}$. As already noted in Ref. Arvanitaki:2024taq , $\langle{\cal J}_{z}\rangle$ benefits from an $N^{2}$ enhancement when the ensemble is initialised in the $|P\rangle$ state. However, this enhancement is partially offset by a suppression arising from the near cancellation between the $\gamma_{+}$ and $\gamma_{-}$ coefficients—typically $\gamma_{\rm net}\sim 2(m_{\nu}\omega_{0}/p_{\nu}^{2})\gamma_{\pm}\sim{\cal O}(10^{-2}-10^{-3})\gamma_{\pm}$, and can be even more strongly suppressed for smaller values of $\omega_{0}$ and $m_{\nu}$, as illustrated in Fig. 1.

In a similar vein, the second-order observables evolve in the perturbative limit according to

[TABLE]

where $\gamma_{\rm tot}=\gamma_{-}+\gamma_{+}$. Note that in the case of $\langle{\cal J}^{2}_{x}\rangle_{|G\rangle}$, we have assumed the ensemble to be prepared in the ground state: this observable can be advantageous, as it captures the same neutrino-induced signal as $\langle{\cal J}_{z}^{2}\rangle_{|P\rangle}$ (for the $m_{\nu}$ and $\omega_{0}$ values considered here, $2\gamma_{+}\simeq\gamma_{\rm tot}$) but without the need to rotate the spins to the equatorial plane.999As pointed out in Ref. Boyers:2025qgc , pulse imperfections can lead to errors in the manipulation of the collective spin. It is also a quantity directly accessed by experiments like CASPEr that measure the transverse magnetisation. The perturbative solutions (29)-(31) should remain valid as long as the state of the system does not deviate too much from the initial state. For the states considered in this work, we therefore expect the perturbative solutions to hold up to $t\lesssim t_{\rm diss}$, even in the case of the superradiant state.101010In the case of the superradiant state, although perturbation theory formally breaks down at times $t\gtrsim t_{\rm SR}$ Arvanitaki:2024taq , the coherent effect persists up to the dissipation timescale $t_{\rm diss}$, i.e., until the system has evolved significantly away from the initial state $|N/2,0\rangle$. This suggests that, even though the perturbative solution no longer accurately describes the full density matrix, it can still capture the expectation values of the relevant observables within $t_{\rm diss}$.

Inspecting the solutions (29)-(31), we note first of all that the well-known momentum-transfer condition $qR\lesssim 1$ is not, by itself, sufficient to guarantee a coherent effect. One must also select an observable suited to the given initial state. For a first-order observable, measuring $\langle{\cal J}_{z}\rangle$ is advantageous when starting from the $|P\rangle$ state, as the quantum states before and after scattering are not fully orthogonal. This requirement was already emphasised in Ref. Arvanitaki:2024taq . Nevertheless, coherence can also build up dynamically through the action of collective neutrino-spin interactions, even if the initial state has no transverse coherence at $t=0$, such as in the ground state $|G\rangle$. This is evident from Eqs. (30) and (31), where we have shown that measuring $\langle{\cal J}^{2}_{x}\rangle_{|G\rangle}$ produces the same coherent effect as measuring $\langle{\cal J}^{2}_{z}\rangle_{|P\rangle}$.

Secondly, while Eqs. (29)-(31) clearly demonstrate that both the first- and second-order observables benefit from an $N^{2}$ enhancement, to assess the degree of observability we need also to consider the intrinsic quantum noise of the states as quantified by the variance of the observable. For the first-order observable $\langle{\cal J}_{z}\rangle_{|P\rangle}$, the variance is

[TABLE]

where we have supplied at the second equality the perturbative expression. The first, $N/4$ term in the second line is commonly known as the spin projection noise (SPN), and respresents (uncorrelated) noise in the initial state. Figure 3 displays the perturbative solution for $\sigma_{{\cal J}_{z}}$ alongside its numerical solution (see Sec. VI.1) and illustrates that deviations from the SPN value of $N/4$ emerge around the dissipative timescale $t_{\rm diss}$. Therefore, at $t\lesssim t_{\rm diss}$, the signal and noise scale as $\langle{\cal J}_{z}\rangle\sim\gamma_{\rm net}N^{2}$ and $\sigma_{\mathcal{J}_{z}}\sim\sqrt{N}$, respectively, leading to a highly favourable signal-to-noise ratio scaling of $\sim\gamma_{\rm net}N^{3/2}$.

On the other hand, at $t\lesssim t_{\rm diss}$, the second-order observables $\langle{\cal J}_{z}^{2}\rangle_{|P\rangle}$ and $\langle{\cal J}_{x}^{2}\rangle_{|G\rangle}$ have an intrinsic variance of $\sigma^{2}_{{\cal J}^{2}}\sim N^{2}$, such that the signal-to-noise ratio scales as $\sim\gamma_{\rm tot}N$. Thus, while higher-order observables may appear at first glance to contain more information, they do not necessarily imply an improvement in sensitivity to the C$\nu$B signal unless the intrinsic noise of the observable can be reduced, for instance, through squeezing Arvanitaki:2024taq . Thus, we expect the $\langle J_{z}\rangle_{|P\rangle}$ observable to yield the best constraint on the C$\nu$B overdensity because of its favourable scaling with the number of spins.

IV.2 Steady-state solution

At $t\gg t_{\rm diss}$, we expect the ensemble to reach a steady state characterised by saturation of the observable. This steady state is not captured by the perturbative solution (29)–(31), but can be found by solving the Lindblad equation (19) in the limit $\dot{\rho}_{S}=0$. Writing the density matrix in the Dicke basis as $\rho_{S}=\sum_{mm^{\prime}}\rho^{mm^{\prime}}_{S}|j,m\rangle\langle j,m^{\prime}|$ and assuming the stationary solution to be unique (i.e., $\rho^{mm^{\prime}}_{S}=0$ for $m\neq m^{\prime}$), a recursive relation for the diagonal elements $\rho^{m,m}_{S}\equiv\rho^{m}_{S}$ can be obtained,

[TABLE]

for $m\in[-N/2-1,N/2+1]$, and $\rho_{S}^{-N/2-1}$=$\rho_{S}^{N/2+1}=0$. Equation (33) is effectively a linear equation for the entries of the diagonal with a tridiagonal matrix, which can be solved iteratively for $\rho^{m}_{t\to\infty}$, from which we can compute the steady-state values of the observables via

[TABLE]

Note that, since we consider only coherent neutrino interactions, the final state of the spin ensemble implied by Eqs. (33) and (34) is not a thermal equilibrium state one would generically expect the system to tend to: for this to happen local interactions are required, to be discussed in Sec. V.

Figure 4 shows the late-time values of the normalised collective observable $\sigma_{x}\equiv(\langle{\cal J}^{2}_{x}\rangle_{|G\rangle})^{1/2}$, as a function of the number of spins $N$. Although not shown here, we also find that the steady-state value of $\langle{\cal J}_{z}\rangle_{|P\rangle}$ approaches the ground-state expectation $-N/2$ around $t\sim(N\gamma_{\rm net})^{-1}$ in the large $N$ limit, following a similar scaling behaviour in $N$ as $\sigma_{x}$. For all observables convergence is reached for $N\gtrsim\mathcal{O}(10^{3})$, but at larger values of $N$ in those cases where a near cancellation between excitation and de-excitation, i.e., $\gamma_{+}\approx\gamma_{-}$, exists.

V Local noise effects

Our analysis so far has focused exclusively on the coherent interaction between the spin ensemble and the C$\nu$B. Realistic experimental setups, however, are inevitably subject to other, non-neutrino environmental interactions that relax the system, i.e., drive it to a state of thermal equilibrium, by locally exciting or de-exciting the spins, as well as locally dephasing the spin precessions. Amongst these, local dephasing plays a particularly important role.

This form of noise originates from random, uncorrelated fluctuations acting independently on each spin, due to, e.g., spin-spin interactions and magnetic field fluctuations, and, phenomenologically, can be described by the operator $\sqrt{\gamma_{\phi}^{\rm loc}}J_{z}^{\alpha}$, where $\alpha$ labels the individual spin. As illustrated in Fig. 2, the key feature of local dephasing is that it acts incoherently and forces the system to move into lower-$j$ sectors where collective effects such as superradiance are suppressed or entirely lost. The relevant timescale is $t_{\phi}\sim 1/\gamma_{\phi}^{\rm loc}$, which is easily much shorter than the dissipative timescale $t_{\rm diss}$ from Eq. (27). Consequently, local dephasing can wash out coherent signatures induced by the neutrino-spin interaction, effectively setting a limit on the observable impact of the C$\nu$B.

In a purely phenomenological approach, let us assume that relaxation is dominated by local (single-particle) effects, i.e., we do not model spin-spin interactions. This simplification allows us to “repurpose” the Lindblad equation (20) in the incoherent regime to model the relaxation characteristics of a realistic experimental sample. Identifying $\gamma_{z}$ in Eq. (20) with the local dephasing rate $\gamma^{\rm loc}_{\phi}$ and the rates $\gamma_{\pm}$ with the local emission and absortion rates $\gamma^{\rm loc}_{\pm}$, we can tune $\gamma_{\pm,\phi}^{\rm loc}$ to match the experimental longitudinal and transverse relaxation times, $T_{1}$ and $T_{2}$, and the final equilibrium state.

Explicitly, starting from Eq. (20), we take the trace ${\rm tr}[\rho_{S}{\cal J}_{i}]$, with $i=x,y,z$, on both sides of the equation, in order to obtain an evolution equation for the collective magnetisation $\langle{\cal J}_{i}\rangle$ in the $i$th direction,

[TABLE]

where $J_{\pm,z,i}$ are single-spin operators, and $\rho_{S}^{(1/2)}$ is a single-spin density matrix. Solving this equation analytically, we find for the longitudinal magnetisation

[TABLE]

with $\langle{\cal J}_{z}\rangle(0)=0,-N/2$ for the $|P\rangle$ and the $|G\rangle$ state, respectively. The longitudinal magnetisation tends to an equilibrium value of

[TABLE]

at a rate that can be matched to the longitudinal relaxation time $T_{1}$ via

[TABLE]

Similarly, starting from the $|P\rangle$ state, we find the magnetisation in the transverse direction to be

[TABLE]

yielding a match

[TABLE]

to the transverse relaxation time $T_{2}$.

Observe in Eq. (37) that the equilibrium magnetisation $\langle{\cal J}_{z}\rangle^{\text{eq}}$ depends only on the ratio $\gamma^{\text{loc}}_{+}/\gamma^{\text{loc}}_{-}$. This is a consequence of the detailed balance condition $\gamma^{\text{loc}}_{+}/\gamma^{\text{loc}}_{-}=e^{-\beta\omega}$, which ensures that the ensemble approaches the well-known thermal steady state

[TABLE]

as long as all interaction channels couple the system to a thermal reservoir at the same temperature $T=1/\beta$. Equation (41) also highlights the inherent difficulty in achieving a strong initial polarisation: even under cryogenic temperatures and extremely large magnetic fields, the polarisation remains small. Reinstating all physical constants, the thermal polarisation takes the form

[TABLE]

where the numerical estimate assumes the gyromagnetic ratio $\gamma$ of 129Xe nuclei. At ${\cal O}(10^{-3})$, this estimate demonstrates that thermal polarisation alone does not suffice to reach the strongly polarised regime, necessitating the use of advanced techniques to hyperpolarise the spin ensemble.

As established in Eqs. (38) and (40), the local relaxation times, $T_{1}$ and $T_{2}$, depend on the absolute values of the local rates $\gamma^{\text{loc}}_{\pm,\phi}$. In typical solid-state NMR systems, the longitudinal relaxation timescale is much larger than the transverse relaxation timescale, i.e., $T_{1}\gg T_{2}$. This implies $\gamma^{\rm loc}_{\pm}\ll\gamma^{\rm loc}_{\phi}$, and the transverse relaxation time $T_{2}$ is well-approximated by the effective dephasing time set by $1/\gamma^{\rm loc}_{\phi}$.111111In standard NMR terminology, the experimentally-observed dephasing time is denoted $T_{2}^{*}$ and includes contributions from magnetic field inhomogeneities and chemical shifts. In contrast, $T_{2}$ refers to the intrinsic dephasing time due to dipolar spin-spin interactions, and typically satisfies $T_{2}\gg T_{2}^{*}$. For our purposes, we do not distinguish between $T_{2}$ and $T_{2}^{*}$, but adopt the convention $T_{2}=T_{2}^{*}=\gamma_{\phi}^{-1}$, which we always use to refer to the actual (measured) dephasing time. Note however that in a generic (i.e., not specific to solid state) NMR setup, $T_{1}$ can be as small as $T_{2}$. In what follows, we model the $\gamma^{\rm loc}_{\pm}$ and $\gamma_{\phi}^{\rm loc}$ rates according to the ratio $T_{1}/T_{2}$, with the additional stipulation that the ratio $\gamma^{\rm loc}_{+}/\gamma^{\rm loc}_{-}$ should reproduce the equilibrium magnetisation achievable at cryogenic temperatures.

VI Numerical solutions

Having discussed both coherent neutrino and incoherent local effects, we are now in a position to explore a more complete, numerical solution of a master equation that accounts for both effects. The strong separation of timescales established above, $T_{2}\ll T_{1}$, implies that coherent effects such as those induced by neutrino-spin interactions are primarily limited by local dephasing. However, because our observable of interest, $\langle{\cal J}_{z}\rangle_{|P\rangle}$, experiences longitudinal relaxation, we must also account for the effects of $T_{1}$, particularly on the correlations of the signal between different times. Therefore, a minimum master equation that should capture the essential behaviours of the system would read

[TABLE]

where ${\cal L}$ is the Lindblad super-operator of this system, and we remind the reader that $\gamma_{\pm}$ are the C$\nu$B-induced excitation and de-excitation rates, while $\gamma_{\pm}^{\rm loc}$ and $\gamma_{\phi}^{\rm loc}$ are the local absorption/emission and dephasing rates.

Given Eq. (43) and assuming the C$\nu$B-induced excitation and de-excitation rates are tunable by an overdensity factor $\delta_{\nu}$, i.e.,121212The overdensity parameter $\delta_{\nu}$ is strictly speaking not an independent parameter, in the sense that other characteristics of the C$\nu$B must also change when $\delta_{\nu}$ is varied. For example, enhancement of $\delta_{\nu}$ due to local gravitational clustering generally results in a larger typical momentum Ringwald:2004np . For excessively large values of $\delta_{\nu}$ (e.g., $\sim 10^{11}$, such as the current KATRIN $\delta_{\nu}$ limit KATRIN:2022kkv ), physically consistent modelling can in fact lead to a signal that is completely different from the $\delta_{\nu}\sim{\cal O}(1)$ case Bondarenko:2023ukx . We nonetheless adopt the definition (44), firstly in order to facilitate comparison with constraints and forecasted sensitivities in the literature that have been derived under similar assumptions. Secondly, for extremely large values of $\delta_{\nu}$, there is also no model-independent way to consistently relate $\delta_{\nu}$ to other C$\nu$B properties. Consequently, the interpretation of $\delta_{\nu}$ must be limited to a mere characterisation of how away an experimental technique is from achieving detection, rather than a necessarily physically meaningful measure of the local C$\nu$B overdensity.

[TABLE]

we can now identify three distinct dynamical regimes defined by the competition between collective neutrino-induced interactions and local dephasing, in the limit where $T_{1}\gg T_{2}$.

• •

Coherence-dominated regime: When $\gamma_{\pm}\gg\gamma^{\rm loc}_{\phi}\approx T_{2}^{-1}$, coherent effects dominate even at the single-spin level, keeping the dynamics well confined to the symmetric subspace. The system exhibits collective behaviours with enhancement scaling as $\sim N^{2}$, providing optimal conditions for signal-to-noise amplification. This regime is depicted in blue in Fig. 5 and is unattainable in realistic experiments under realistic dephasing, unless $\delta_{\nu}$ is very large.

• •

Intermediate regime: In the case $\gamma_{\pm}\ll\gamma^{\rm loc}_{\phi}$, but $N\gamma_{\pm}\gg\gamma^{\rm loc}_{\phi}$ is still maintained, individual spins are dominated by local dephasing. Yet collective effects are stronger in the global dynamics. Thus, collective effects remain observable, even while symmetry is partially broken by local noise. This regime is depicted in Fig. 5 in green, highlighting a complete suppression of the coherent signal only at $t\sim\gamma^{-1}_{\pm}$. The intermediate regime is also the most relevant experimentally for detecting weak collective signatures under realistic dephasing, although an enhanced $\delta_{\nu}\gg{\cal O}(1)$ will still be necessary.

• •

Dephasing-dominated regime: When $N\gamma_{\pm}\ll\gamma^{\rm loc}_{\phi}$, incoherent processes dominate over the collective interaction. In this regime, the system rapidly departs from the symmetric subspace, the signal-to-noise ratio saturates well below unity, and sensitivity to the neutrino–spin interaction is effectively lost. This region, shown in orange in Fig. 5, corresponds to the case of a standard C$\nu$B, where $\delta_{\nu}\sim{\cal O}(1)$ is the generic expectation.

In the following, we first solve Eq. (43) in Sec. VI.1 numerically by way of an efficient method that tracks explicitly the observables of interest: we call it the “second-order approximation” because of a truncation that needs to be applied at second order. This is to be compared to a complete numerical solution for the full density matrix $\rho_{S}$, whose details can be found in Appendix F. In Sec. VI.2 we apply the same method to investigate the evolution of the two-time correlation functions under Eq. (43), which will be used later in Sec. VII to form the covariance matrix that quantifies the theoretical uncertainties in the observables.

VI.1 Second-order approximation

While it is possible to numerically solve the Lindblad master equation (43) at the level of the density matrix for a system of $N$ spins, it can quickly become impractical or even intractable for values of $N$ larger than about 100 (see Appendix F). Rather than dealing with the full density matrix, then, we shift the focus to tracking just a handful physically meaningful quantities.

Specifically, as was done in, e.g., Refs. Shammah_2017 ; carmichael1 , for the case of Dicke superradiance and also touched upon in Sec. V, we construct differential equations from the master equation (43) directly for the expectation values of the relevant observables, namely, the collective magnetisation $\langle{\cal J}_{z}\rangle$, $\langle{\cal J}_{z}^{2}\rangle$, and $\langle{\cal J}_{x}^{2}\rangle$. This is akin to tracking the moments of a continuous system, which is widely used in many different fields. However, differential equations thus constructed generally do not form a closed system, in the sense that the equation of motion for, e.g., $\langle{\cal J}_{z}\rangle$ depends on a higher moment $\langle{\cal J}_{z}^{2}\rangle$, whose own equation of motion depends in turn on $\langle{\cal J}_{z}^{3}\rangle$, and so on, giving rise an infinite hierarchy of coupled equations. To solve the system, therefore, one needs to truncate the hierarchy at a chosen order Shammah_2017 .

Here, we opt to truncate the hierarchy at second-order. For the observables $\langle{\cal J}_{z}\rangle$ and $\langle{\cal J}_{z}^{2}\rangle$, this means we adopt the closure conditions $\langle{\cal J}_{z}^{3}\rangle\approx\langle{\cal J}_{z}\rangle\langle{\cal J}_{z}^{2}\rangle$ and $\langle{\cal J}_{z}{\cal J}^{2}\rangle\approx\langle{\cal J}_{z}\rangle\langle{\cal J}^{2}\rangle$ to arrive at Shammah_2017

[TABLE]

Note that our closure conditions are by no means unique and other choices exist. For example, closure can be achieved alternatively by demanding that all cumulants (i.e., connected pieces of the moments) vanish at third order and above, e.g., $\langle{\cal J}_{z}^{3}\rangle=3\langle{\cal J}_{z}^{2}\rangle\langle{\cal J}_{z}\rangle+2\langle{\cal J}_{z}\rangle^{3}+\langle{\cal J}_{z}^{3}\rangle_{c}$, where $\langle{\cal J}_{z}^{3}\rangle_{c}$ is the connected piece to be set to zero. We observe however that all closure-dependent terms in equation (45) are proportional to $\gamma_{\rm net}$ and vanish in the limit $\gamma_{\rm net}\to 0$. Indeed, given the small $\gamma_{\rm net}\ll\gamma_{\pm}$ in our scenario, we have verified that adopting a different closure condition (e.g., the cumulant truncation scheme mentioned above) has no discernible impact on our numerical solutions within in the time range of interest (i.e., $t\lesssim T_{2}$) even for extremely large $N$ values.

Initial conditions including partial polarisation

As discussed in Sec. III.2, in experimental settings spin ensembles are prepared in a hyperpolarised configuration that approximates the ground state and, if need be, rotated to the desired orientation by applying an appropriate pulse. Common hyperpolarisation protocols such as spin-exchange optical pumping achieve close-to-unit polarisation by way of spin excitation and de-excitation processes that can be modelled as local interactions. See Appendix G for details. Thus, in order to establish the initial conditions of a partially-polarised spin ensemble prepared in the aforementioned manner, we can simply take Eq. (45) assuming only the local rates $\gamma_{\pm}^{\rm loc}$ are non-zero, re-interpret these rates as the excitation and de-excitation rates $\gamma_{\pm}^{\rm hyper}$ of the hyperpolarisation technique in use, and compute the corresponding steady-state solutions for $\langle{\cal J}_{z}\rangle$ and $\langle{\cal J}_{z}^{2}\rangle$.

These solutions then serve as initial conditions for the partially-polarised “ground state”, which read

[TABLE]

where we have used $\langle{\cal J}^{2}\rangle=\langle{\cal J}_{x}^{2}\rangle+\langle{\cal J}_{y}^{2}\rangle+\langle{\cal J}_{z}^{2}\rangle=N/2+\langle{\cal J}_{z}^{2}\rangle$, and

[TABLE]

is the initial polarisation of the ensemble expressed in terms of the rates of the hyperpolarisation method. See also Eq. (37). If on the other hand we wish to initialise in a product state-like configuration, a $\pi/2$ rotation of the coordinates should be applied to Eq. (46) to get

[TABLE]

while $\langle{\cal J}^{2}\rangle$ remains the same under rotation. Then, together, Eqs. (45) and (46) or (48) form a closed set of equations that completely determine $\langle{\cal J}_{z}\rangle,\langle{\cal J}^{2}_{z}\rangle$, and $\langle{\cal J}^{2}\rangle$.

To obtain approximate numerical solutions for the transverse observable $\langle{\cal J}^{2}_{x}\rangle$, we supplement the differential equations (45) with an additional equation of motion for the non-hermitian variables $\langle{\cal J}_{\pm}\rangle\equiv\langle{\cal J}_{x}\rangle\pm{\rm i}\,\langle{\cal J}_{y}\rangle$, i.e.,

[TABLE]

where we have applied the closure conditions $\langle{\cal J}_{+}{\cal J}_{z}\rangle\approx\langle{\cal J}_{+}\rangle\langle{\cal J}_{z}\rangle$ and $\langle{\cal J}_{z}{\cal J}_{-}\rangle\approx\langle{\cal J}_{z}\rangle\langle{\cal J}_{-}\rangle$. Here, we are interested in a spin ensemble initially prepared in the ground state $|G\rangle$, such that the initial conditions are $\langle{\cal J}_{\pm}(0)\rangle_{|G\rangle}=0$. Then, a solution for $\langle{\cal J}^{2}_{x}\rangle$ can be constructed from

[TABLE]

where $\langle{\cal J}_{z}^{2}\rangle$, $\langle{\cal J}^{2}\rangle$, and $\langle{\cal J}_{\pm}\rangle$ are solutions of Eqs. (45) and (49).

Figures 6-8 show the numerical solutions obtained from the above equations of motion under different settings, which we discuss in more detail below.

The ideal case: no local noise and perfect polarisation

We consider first the case in which there is no local noise and the polarisation of the initial state is perfect, i.e., $p=1$. Figure 6 shows the solutions for the collective observable $\langle{\cal J}_{z}\rangle$ under these settings for several choices of $N$. These are contrasted with numerical solutions of the full density matrix for those choices of $N$ values where a full numerical solution is tractable. Details of the full density matrix solution can be found in Appendix F.

Remarkably, where the full density matrix solution can be found, the second-order solution is in excellent agreement with it at all times. This agreement follows from the small $\gamma_{\rm net}\ll\gamma_{\pm}$ of our system, where, as discussed above, the hierarchy of equations (45) becomes exact and independent of the truncation scheme in the limit $\gamma_{\rm net}\to 0$, and indicates that the second-order method is indeed a reliable way to obtain predictions for the collective observables. Importantly, the method remains valid for arbitrarily large $N$ values—$N=10^{22}$ in this figure—at no additional cost in terms of computational runtime or memory relative to the small $N$ case. In contrast, the computational cost of numerically solving the full density matrix scales as $\propto N^{3}$, and becomes prohibitively expensive even on modern computers when $N\gtrsim 10^{4}$. For completeness, we also show in the same figure the perturbative solution (29). Evidently, our numerical solutions are well approximated by the perturbative solution at $t\lesssim(N\gamma_{\rm net})^{-1}$, especially for large $N$ values.

With local noise and partial polarisation

Having established that the second-order method is reliable, we now use it to investigate more realistic cases in which local noise is present and the polarisation of the initial state may only be partial.

Figure 7 illustrates the impact of non-zero dephasing on the observable $\langle{\cal J}_{z}\rangle$ for a large system size $N$ (red curve). As shown, dephasing suppresses the growth of $\langle{\cal J}_{z}\rangle$ at late times. In this example, the suppression begins around the dissipative timescale $t_{\rm diss}$ due to our specific choice of local dephasing rate, $\gamma_{\phi}^{\rm loc}=N\gamma_{-}$. Different choices for this rate would lead to alternative behaviours, as previously demonstrated in Fig. 5. While dephasing does not directly alter the dynamics of $\langle{\cal J}_{z}\rangle$, it inhibits the coherent dissipation mediated by neutrino interactions. This dissipation then drives the system toward its own equilibrium value, given by $\langle{\cal J}_{z}\rangle=-(N/2)\gamma_{\rm net}/\gamma_{\rm tot}$ in the limit $T_{1}\to\infty$.

Finite $T_{1}$ effects are introduced via the local emission and absorption rates $\gamma_{\pm}^{\rm loc}$ in the master equation. We select a ratio $\gamma_{+}^{\rm loc}/\gamma_{-}^{\rm loc}=0.995$, corresponding to an equilibrium polarisation of approximately $p_{\rm eq}\sim 10^{-3}$—a value achievable at cryogenic temperatures. The strengths of the local relaxation processes are chosen such that $T_{1}$ is ten times longer than $T_{2}$, and we allow this ratio to vary in the forecasting analysis in Sec. VII. The blue and green curves in Fig. 7 show the effects of local relaxation and dephasing, with and without superradiant neutrino interactions, respectively.

Figure 8 shows the variance $\sigma_{{\cal J}_{z}}^{2}$ of ${\cal J}_{z}$ and the transverse observable $\langle J_{x}^{2}\rangle_{|G\rangle}$. Here, in addition to demonstrating the effects of local dephasing, we also show how partial polarisation affects the time evolution of the observables. Clearly, in the presence of dephasing, the variance $\sigma_{{\cal J}_{z}}^{2}$ shows minimal growth unless $p>0.1$. In the case of $\langle J_{x}^{2}\rangle_{|G\rangle}$, the polarisation requirement is even stronger—$p>0.3$—although, as discussed in Sec. IV, we do not expect this observable to be competitive against $\langle J_{z}\rangle_{|P\rangle}$. This result therefore highlights the importance of polarisation for the purpose of C$\nu$B detection. For instance, a liquid methanol-based setup Walter:2025ldb , with a polarisation of order $p\sim 10^{-5}$, is currently insufficient for NMR experiments to target C$\nu$B effects. More generally, thermal polarisation levels achievable at typical experimental conditions—${\cal O}(10^{-3})$ for 129Xe from Eq. (42)—remain too low to yield detectable signals. Advanced polarisation techniques are thus essential to this venture.

VI.2 Two-time correlation functions

We envisage a future experiment that measures an observable, e.g., $\langle{\cal J}_{z}\rangle_{|P\rangle}$, at a series of time points. To estimate the sensitivity to the C$\nu$B overdensity parameter of such an experiment, besides the observable’s variance $\sigma^{2}(t)$ at each time $t$, we will also need its covariance $\sigma^{2}(t,t+\tau)$ at two different times. This covariance can be constructed from the two-time correlation function, generally written as $\langle{\cal O}_{1}(t){\cal O}_{2}(t+\tau)\rangle$, where ${\cal O}_{1}$ and ${\cal O}_{2}$ are two operators that act only on the spins.

Following the quantum regression theorem QRT , the two-time correlation function can be computed as

[TABLE]

where the modified density matrix $\rho_{S;{\cal O}_{1}}$ is defined at $t$ via $\rho_{S;{\cal O}_{1}}(t)=\rho_{S}(t){\cal O}_{1}(0)$. The evolution of $\rho_{S;{\cal O}_{1}}(t+\tau)$ from $t$ to $t+\tau$ for $\tau\geq 0$ is governed by the same Lindblad master equation (43) used to compute our observables but for a change of the time variable $t\to\tau$, i.e.,

[TABLE]

where ${\cal L}$ is the Lindblad super-operator of Eq. (43). A more detailed discussion about this equation of motion and the construction of two-time correlation functions can be found in Appendix H.

Since our main observable is $\langle{\cal J}_{z}\rangle_{|P\rangle}$, we are particularly interested in the two-time correlation function $\langle{\cal J}_{z}(t){\cal J}_{z}(t+\tau)\rangle$. Following the second-order method of Sec. VI.1—or, equivalently, using the correspondence between Eqs. (159) and (160)) in Appendix H—a hierarchy of differential equations can be constructed from Eq. (52) to describe its time evolution, namely,

[TABLE]

where we have truncated the hierarchy by approximating $\langle{\cal J}_{z}(t){\cal J}_{z}^{3}(t+\tau)\rangle\approx\langle{\cal J}_{z}(t){\cal J}_{z}^{2}(t+\tau)\rangle\langle{\cal J}_{z}(t+\tau)\rangle$ and $\langle{\cal J}_{z}(t){\cal J}_{z}(t+\tau){\cal J}^{2}(t+\tau)\rangle\approx\langle{\cal J}_{z}(t){\cal J}^{2}(t+\tau)\rangle\langle{\cal J}_{z}(t+\tau)\rangle$, with $\langle{\cal J}_{z}(t+\tau)\rangle\equiv{\rm tr}[\rho_{S}(t+\tau){\cal J}_{z}]$. We also remind the reader that the time variable here is $\tau$, and the function $\langle{\cal J}(t)\rangle$ which appears in all three equations above does not evolve with $\tau$.

Equation (53) can be solved from any time point $t$ in the duration of an experiment, given the “initial conditions” at $\tau=0$, $\langle{\cal J}_{z}(t)\rangle$, $\langle{\cal J}_{z}^{2}(t)\rangle$, $\langle{\cal J}^{2}(t){\cal J}_{z}(t)\rangle$, and $\langle{\cal J}_{z}^{3}(t)\rangle$. For a spin ensemble originally prepared in the $|P\rangle$ state, these initial conditions amount at $t=0$ to

[TABLE]

For $t>0$, on the other hand, $\langle{\cal J}_{z}(t)\rangle$ and $\langle{\cal J}_{z}^{2}(t)\rangle$ can be obtained directly from the solutions to Eq. (45), while we use the same solutions to approximate $\langle{\cal J}^{2}(t){\cal J}_{z}(t)\rangle\approx\langle{\cal J}^{2}(t)\rangle\langle{\cal J}_{z}(t)\rangle$ and $\langle{\cal J}_{z}^{3}(t)\rangle\approx\langle{\cal J}_{z}(t)\rangle\langle{\cal J}_{z}^{2}(t)\rangle$.

We are particularly interested in the covariance of the observable $\langle{\cal J}_{z}\rangle_{|P\rangle}$ at different times, i.e, the connected piece of the two-time correlation function, defined as

[TABLE]

which we will use in Sec. VII to characterise the theoretical uncertainties in the observable. This connected piece can be obtained from the solutions of Eqs. (45) and (53) for every pair of $(t,t+\tau)$. It is however also possible—and instructive—to construct an equation of motion directly for $\sigma^{2}_{{\cal J}_{z}}(t,t+\tau)$ from Eqs. (45) and (53), which reads

[TABLE]

Where the $\propto\gamma_{\rm net}$ terms are suppressed relative to the $\propto(\gamma_{\rm tot}+\gamma_{-}^{\rm loc}+\gamma_{+}^{\rm loc})$ term, this equation of motion is solved approximately by

[TABLE]

where we have used $\gamma_{\pm}\ll\gamma_{\pm}^{\rm loc}$ at the second equality, and $\sigma^{2}_{{\cal J}_{z}}(t)\equiv\langle{\cal J}_{z}^{2}(t)\rangle-\langle{\cal J}_{z}(t)\rangle^{2}$ is the usual variance at time $t$ also defined earlier in Eq. (32).

Figure 9 shows the time evolution of the two-time correlation function $\sigma^{2}_{{\cal J}_{z}}(t,t+\tau)$, obtained from solving Eq. (53), for $N=10^{3}$ and $t=0$. In the case $\gamma_{-}^{\rm loc}=N\gamma_{-}$ (dark blue solid line), we see an exponential decay at a characteristic time $\tau\sim 1/(2N\gamma_{-})\approx 1/(\gamma_{-}^{\rm loc}+\gamma_{+}^{\rm loc})$, while in the absence of local emission/absorption the correlation decays at $\tau\sim 1/(2\gamma_{-})\approx 1/\gamma_{\rm tot}$ due to neutrino dissipation alone, confirming the analytical estimate (57). Observe also that local dephasing effects do not influence the decay in the longitudinal correlation. This is primarily a consequence of the small $\gamma_{\rm net}$ in our scenario, which inhibits the effects of $\gamma_{\phi}^{\rm loc}$ on the evolution of $\sigma_{{\cal J}_{z}}^{2}$. Lastly, while Fig. 9 presents the numerical solutions for $t=0$, the two-time correlation for other choices of $t$ exhibits the same decay behaviour, as suggested by Eq. (57).

VII Future constraints on the C$\nu$B

To assess the sensitivity of future experiments to collective effects of the C$\nu$B, we consider an idealised NMR experiment using polarised 129Xe nuclei, as envisioned in future stages of the CASPEr experiment Walter:2025ldb . 129Xe is a noble gas isotope with a nuclear spin of $I=1/2$, making it naturally suited to modelling as a spin-$1/2$ system. In liquid form Xenon has a density of $\rho=2.94$ g/cm3, which gives

[TABLE]

as the number density of Xenon nuclei given an atomic weight of $A_{r}({\rm Xe})=131.29$ Prohaska . Assuming 100% isotopic enrichment of 129Xe can be achieved, we use the same number to approximate the 129Xe spin density $n_{s}$.

The gyromagnetic ratio of 129Xe is approximately $\gamma_{\text{Xe}}\simeq 11.8~\text{MHz/T}$, which gives a maximum energy splitting of $\omega_{0}\simeq 5.8\times 10^{-7}~\text{eV}$ for a maximum magnetic field strength of $B_{\rm max}=12$ T. For a fixed neutrino mass, the C$\nu$B-induced net interaction rate $\gamma_{\rm net}\propto 1-\gamma_{+}/\gamma_{-}\sim 2m_{\nu}\omega_{0}/p_{\nu}^{2}$ generally increases with $\omega_{0}$, as also shown in Fig. 1. However, $\omega_{0}$ cannot be made arbitrarily large: for the typical momentum transfer $q_{\rm typ}\sim\omega_{0}m_{\nu}/p_{\nu}$, the coherence condition must eventually break down when $q_{\rm typ}$ reaches $\sim R^{-1}$. Then, setting the coherence limit at $q_{\rm typ}R<1$ immediately yields an upper bound on the energy splitting $\omega_{0}<p_{\nu}/m_{\nu}R$. Thus, in our modelling, the optimal energy splitting is determined by the smaller of the maximum value achievable in the laboratory and the coherence limit, namely,

[TABLE]

where we have set $m_{\nu}=m_{1}$. Assuming the normal hierarchy, this choice of $m_{\nu}$ gives the largest possible value of $\omega^{\star}_{0}$, leading to the best sensitivity to $\delta\nu$.

As already mentioned in Sec. III.2, the preparation of the ground state, with all spins aligned along the magnetic field direction (i.e., the $-z$-direction in this work), can be achieved via spin-exchange optical pumping. This technique is commonly used to hyperpolarise noble gases like Xenon PhysRevLett.87.067601 ; RevModPhys.69.629 and can reach polarisations up to 90% doi:10.1021/jp501493k . Notably, hyperpolarised 129Xe exhibits exceptionally long transverse relaxation times, with values estimated to be up to $T_{2}\sim 1000$ s PhysRevLett.87.067601 ; Walter:2025ldb . After the nuclear spin ensemble has attained the desired polarisation, a coherent spin (product) state can be prepared by applying a resonant $\pi/2$ radio-frequency pulse at the Larmor frequency. This effectively rotates the collective spin vector into the equatorial plane (see Fig. 2), thereby allowing us to use the $\langle{\cal J}_{z}\rangle_{|P\rangle}$ observable to constrain collective neutrino effects. Under certain conditions, the longitudinal relaxation time $T_{1}$ can also be long, e.g., $T_{1}\gtrsim$ hours for solid-state 129Xe PhysRevB.94.094309 . In the following we take $T_{1}/T_{2}=10$, which corresponds to a value $T_{1}\simeq 2.8$ hours for $T_{2}=1000$ s, modelled according to Eqs. (38) and (40).

VII.1 $\chi^{2}$ analysis

To estimate the sensitivity to the C$\nu$B overdensity parameter $\delta_{\nu}$, we simulate a mock experimental data set ${\cal\hat{J}}_{z}(t_{i})$ of a fiducial model that does not contain the coherent neutrino signal,131313A precessing transverse magnetisation can in principle induce a current in the detection coil and generate a feedback magnetic field that accelerates the decay of the magnetisation. This effect is commonly known as radiation damping and may exhibit a collective behaviour scaling. In this work, we adopt a fiducial model that does not include such collective background effects, which can be justified under certain conditions Galanis:2025amc . A full modeling of the experimental environment lies beyond the scope of this study. comprising $N_{\text{shots}}$ measurements at $n_{t}$ discrete time points sampled from $t_{1}=1/f_{s}$ to $t_{n_{t}}=T_{2}$. The sampling frequency is fixed to $f_{s}=14.3$ kHz, in line with the most recent CASPEr setup Walter:2025ldb . Then, assuming the likelihood function (i.e., probability of data given a theory) to be gaussian, we can define a $\chi^{2}$-statistics via

[TABLE]

where $C$ is the covariance matrix to be defined below, and

[TABLE]

is theoretical prediction minus the mock data at time $t_{i}$.

The theoretical prediction $\langle{\cal J}_{z}(t_{i};\delta_{\nu})\rangle_{|P\rangle}$ includes neutrino effects tuned by the overdensity parameter $\delta_{\nu}$. To ensure consistency, we take $\langle{\cal J}_{z}(t_{i};\delta_{\nu})\rangle_{|P\rangle}$ from the numerical solutions of Sec. VI, with the time variable rescaled by a factor $(2p_{\nu}R)^{2}$ to account for partial coherence in the $2p_{\nu}R>1$ regime (note that the local rates $\gamma_{\pm,\phi}^{\rm loc}$ must also be simultaneously rescaled by the same factor in order to preserve their effects). For the mock data ${\cal\hat{J}}_{z}(t_{i})$, on the other hand, we assume for simplicity that they coincide with the expectation value $\langle{\cal\hat{J}}_{z}(t_{i})\rangle_{|P\rangle}$ under the fiducial model (i.e., without neutrino effects), as given by Eq. (36) with $\gamma_{-}^{\rm loc}\simeq 0.026\gamma_{\phi}^{\rm loc}$ and $\gamma_{+}^{\rm loc}/\gamma_{-}^{\rm loc}=0.995$ (i.e., blue line in Fig. 7).

The covariance matrix $C_{ij}=C_{ji}$ is given by

[TABLE]

and includes contributions from the total instrumental noise $\sigma_{\rm SQUID}^{2}$, as well as the theoretical covariance of the model including neutrinos,

[TABLE]

where $\sigma^{2}_{{\cal J}_{z}}(t_{i};\delta_{\nu})$ is the theoretical variance of ${\cal J}_{z}$, and we have used the approximation (57) for the two-time correlation functions $\sigma^{2}_{{\cal J}_{z}}(t_{i},t_{j};\delta_{\nu})$ in the off-diagonal ($i\neq j$) entries. Again, we take $\sigma^{2}_{{\cal J}_{z}}(t_{i};\delta_{\nu})$ from the numerical solutions of Sec. VI, with a rescaling factor $(2p_{\nu}R)^{2}$ applied to the time scale and the local rates $\gamma_{\pm,\phi}^{\rm loc}$.

The total instrumental noise, which we refer to as “SQUID noise”, includes intrinsic SQUID contributions, thermal (Johnson–Nyquist) noise from the input circuit, as well as other noise sources originating from, e.g., the readout strategy. These noises are assumed to be random, uncorrelated between independent measurements, and uncorrelated with the signal. We model their total variance as

[TABLE]

so that the ratio

[TABLE]

quantifies the ratio of instrumental noise relative to the spin projection noise. The SQUID and SPN power spectral densities, $\mathcal{P}_{\rm SQUID}(f)$ and $\mathcal{P}_{\rm SPIN}(f)$, are specific to the experimental setup. For the CASPEr experiment Walter:2025ldb , we estimate that $\mathcal{F}_{\rm SQUID}\sim 85$.141414Using the specifications given in Ref. Walter:2025ldb , we estimate that at $f=1.2$ MHz (corresponding to $B=0.1$ T for 129Xe), the SQUID spectral density is approximately independent of frequency and takes on the value of $\mathcal{P}_{\rm SQUID}\sim 10^{-12}~\Phi_{0}^{2}/\mathrm{Hz}$. The SPN power is on the other hand $\mathcal{P}_{\rm SPN}\sim 1.2\times 10^{-14}~\Phi_{0}^{2}/\mathrm{Hz}$. Combining these numbers per Eq. (65) yields $\mathcal{F}_{\rm SQUID}\sim 85$. In the following analysis, however, we shall retain $\mathcal{F}_{\rm SQUID}$ as a free parameter.

VII.2 Forecasted constraints

By construction, the $\chi^{2}$ of Eq. (60) is automatically zero when $\delta_{\nu}=0$. Then, to determine an upper limit on $\delta_{\nu}$ we simply need to find the smallest value of $\delta_{\nu}$ that satisfies

[TABLE]

where $\Delta\chi^{2}=2.71~(3.84)$ corresponds to a 90% (95%) confidence level for one degree of freedom.

We consider first the best possible constraint on $\delta_{\nu}$. This can be achieved by assuming perfect polarisation $p=1$ and an optimal splitting energy $\omega^{\star}_{0}$ as defined in Eq. (59). In this idealised scenario, we find that a 10 cm-radius spherical sample of 129Xe nuclei would be able to constrain the C$\nu$B overdensity at

[TABLE]

for the normally-ordered $\sum m_{\nu}=0.15$ eV benchmark scenario. The limit assumes a coherence time of $T_{2}=10^{3}$ s, $T_{1}=10~T_{2}$, and $N_{\rm shots}=10^{3}$ along the lines of Ref. Arvanitaki:2024taq , and is comparable to the estimated sensitivity quoted therein (which assumes hydrogen spins) within an order of magnitude. At face value, such an optimistic sensitivity appears uncompetitive against the current KATRIN limit on the neutrino overdensity of $\delta_{\nu}\lesssim 10^{11}$ for the benchmark scenario (but may be marginally better for a small range of neutrino masses, depending on the ratio $T_{1}/T_{2}$; see Fig. 10). Note however that this conclusion may change if the system achieves at least an $\mathcal{O}(10)$ value in squeezing Galanis:2025amc .

Realistic NMR experiments face even more limitations: imperfect initial polarisation, smaller sample volumes, and instrumental noise. Assuming an experimental setup likely to be realised by CASPEr—instrumental-noise-limited and using a smaller sample volume of $V\sim R^{3}\sim{\rm cm}^{3}$—the sensitivity to $\delta_{\nu}$ reduces significantly to

[TABLE]

at 90% confidence level, again for the normally-ordered $\sum m_{\nu}=0.15$ eV benchmark.

Figure 10 shows the dependence of the projected $\delta_{\nu}$ sensitivities on the lightest neutrino mass $m_{1}$, assuming normal ordering. From $m_{1}\simeq 0.05$ eV up to the KATRIN neutrino mass limit $<0.45$ eV KATRIN:2024cdt , future setups following the specifications of CASPEr will be able to set limits of $\delta_{\nu}\lesssim 10^{12}-10^{14}$, depending on the SQUID noise floor and provided that the initial sample polarisation reaches ${\cal O}(1)$ values. Because the sensitivity to $\delta_{\nu}$ increases with neutrino mass for $m_{1}\gtrsim 0.05$ eV, our forecasted CASPEr constraints are generally better than the benchmark value (68) for larger neutrino masses.

For $m_{1}\lesssim 0.05$ eV, the mass dependence in the rates is lost, and the sensitivities to $\delta_{\nu}$ asymptote to $\delta_{\nu}\sim 10^{13}$. On the other hand, while the $m_{1}$ dependence also drops out in the optimistic ($R=10$ cm) constraint at very small masses, we still observe a somewhat favourable scaling in the range $m_{1}\in(2\times 10^{-3},3\times 10^{-2})$ eV. This arises because, for a larger sample size $R=10$ cm, the optimal energy splitting can be tuned to higher values and only reaches the ceiling imposed by the maximum magnetic field $B_{\rm max}$ at around $m_{1}\sim 10^{-3}$ eV. The effect of this ceiling can be clearly discerned in Fig. 10, where the optimistic sensitivity plateaus at a value of $\delta_{\nu}\sim 5\times 10^{10}$ for $m_{1}\lesssim 10^{-3}$ eV.

We display in Fig. 11 the parametric dependence of the experimental sensitivity to $\delta_{\nu}$ on the system size $R$ and the transverse coherence time $T_{2}$ in the SQUID-noise-dominate regime (i.e., ${\cal F}_{\rm SQUID}\gg 1$). As can be seen, when the condition $p_{\nu}R<1$ is satisfied, the system is in the fully coherent regime, where the interaction rates $\Gamma_{\pm}$ scale as $\propto R^{6}$. Taking into account the scaling of the total uncertainty, $\sigma\propto N^{1/2}\propto R^{3/2}$, this leads to a scaling of $R^{-9/2}$ for the $\delta_{\nu}$ sensitivity. This behaviour applies to both a fixed energy splitting and the optimal case, since the optimal splitting here is effectively constant and determined by the maximum magnetic field achievable in the experiment. As the radius increases, coherence is eventually lost around $R\sim p_{\nu}^{-1}$ except in the forward direction, resulting in a reduced interaction rate that scales as $\Gamma_{\pm}\propto R^{4}$ and hence a $\delta_{\nu}$ constraint scaling of $R^{-5/2}$. This scaling remains the same even at $p_{\nu}R\gg 1$ in the case of a fixed energy splitting.

In the case of optimal splitting, however, a change in the scaling behaviour occurs around $R\sim 0.5$ cm for 129Xe spins. This transition arises because, at these radii, values of the magnetic field required to tune $\gamma B_{\rm max}$ to the optimal value $p_{\nu}/(m_{1}R)$ are achievable experimentally. As such, the energy splitting $\omega_{0}$ becomes dependent on $R$, which leads the net rate $\gamma_{\rm net}\approx(2m_{1}\omega_{0}/p_{\nu}^{2})\gamma_{\pm}$ itself to directly scale with $R^{-1}$ at the optimal $\omega_{0}$. Combined with the loss of coherence at $p_{\nu}R\gg 1$, the collective net interaction rate here scales as $\propto R^{3}$. Consequently, the sensitivity to $\delta_{\nu}$ scales as $R^{-3/2}$ in this regime. Note that, had we used a different target material with a different gyromagnetic ratio, the transition in the scaling behaviour would occur at a different value of $R$.

In the right panel of Fig. 11, we present the $T_{2}$ scalings. In the standard quantum-limited scenario, where SQUID noise can be suppressed, the sensitivity to $\delta_{\nu}$ scales as $T_{2}^{-1}$. For a fixed total integration time $t$, this is equivalent to the standard sensitivity scaling of $\sqrt{tT_{2}}$. In contrast, when uncorrelated SQUID noise dominates the uncertainties, we observe a more favourable scaling of $\delta_{\nu}\propto T_{2}^{-3/2}$. In the case of $T_{1}$, we find a scaling of $\delta_{\nu}\propto(T_{1}/T_{2})^{-1/2}$. Scaling with other parameters can be straightforwardly deduced from Eq. (60). Combining all information, we find that the sensitivity to $\delta_{\nu}$ in the SQUID-noise-dominated regime scales as

[TABLE]

where $r$ denotes the $R$ dependence from Fig. 11, which can be parameterised as

[TABLE]

This shows, as expected, that the sensitivity to $\delta_{\nu}$ improves for larger sample sizes and higher spin densities, as well as longer coherence and measurement times. The sensitivity is, however, expected to degrade quadratically with imperfect polarisation and is further limited by instrumental noise.

Finally, it is worth mentioning that one could also constrain $\delta_{\nu}$ using the second-order observable $\langle{\cal J}_{x}^{2}\rangle_{|G\rangle}$. While this observable has the advantage of probing a C$\nu$B-induced signal without requiring spin rotation into the equatorial plane, its signal-to-noise ratio scales as $\sim\gamma_{\rm tot}N$, as pointed out in Sec. IV. As a result, it provides a weaker sensitivity to $\delta_{\nu}$ compared to $\langle{\cal J}_{z}\rangle_{|P\rangle}$ for spin ensembles with $N\gtrsim(\gamma_{\rm tot}/\gamma_{\rm net})^{2}\simeq 10^{6}$. Indeed, even in the optimisitc $R=10$ cm scenario, the sensitivity of $\langle{\cal J}_{x}^{2}\rangle_{|G\rangle}$ to the CMB overdensity is only $\delta_{\nu}\gtrsim 10^{20}$, about ten orders of magnitude weaker than what can be achieved with $\langle{\cal J}_{z}\rangle_{|P\rangle}$.

VIII Conclusions

We have investigated in this work the dynamics of large spin ensembles interacting with the cosmic neutrino background, with a particular focus on identifying the conditions under which coherent collective effects—scaling as $N^{2}$, where $N$ is the number of spins—can amplify the signal from an otherwise extremely weak neutrino-spin coupling. Such collective neutrino-spin interactions are expected to contribute to dissipating a spin ensemble and can potentially be exploited as a competitive way to set limits on the local C$\nu$B overdensity, as first put forward in Ref. Arvanitaki:2024taq .

After briefly reviewing the interaction rates of the C$\nu$B with a system of $N$ two-level systems (e.g., spin-$1/2$ particles) and the Lindblad master equation in the open quantum system framework, we reproduced the perturbative solutions presented in Ref. Arvanitaki:2024taq , valid when the spin ensemble has not deviated too much from its initial state. We then extend the analysis by considering steady-state analytical solutions as well as numerical simulations of the Lindblad master equation for different initial spin configurations. Notably, our numerical simulations account for realistic noise sources that can degrade the collective C$\nu$B signal, such as local dephasing as well as imperfect initial spin polarisation. We make use of the Dicke basis—illustrated in Fig. 2—to facilitate this analysis, which provides an intuitive representation of how collective and local interactions influence the spin ensemble’s evolution under the Lindblad master equation.

We have furthermore presented a fast and computationally inexpensive method to derive approximate numerical solution from the Lindblad equation. Specifically, the method solves directly for the observables and correlation functions of interest while capturing both coherent neutrino effects and local interactions, and is able to reproduce the full density matrix solution to ${\cal O}(0.01)\%$ accuracy. This fast method enables us to use statistical inference methods to forecast the sensitivity of future NMR experiments to the C$\nu$B overdensity, while incorporating the complete quantum evolution of the spin ensemble.

Considering the projected experimental specifications of the forthcoming CASPEr experiment Budker:2013hfa ; JacksonKimball:2017elr ; Walter:2025ldb and nuclear-spin-based experiments more generally, we compute future experimental sensitivities to the C$\nu$B overdensity parameter $\delta_{\nu}$ using a $\chi^{2}$ analysis. We find that achievable near-future setups along the lines of CASPEr could constrain the C$\nu$B overdensity parameter at the level of $\delta_{\nu}\sim 10^{13}$ for neutrino masses currently within the KATRIN neutrino mass bound KATRIN:2024cdt , with larger masses yielding a more favourable sensitivity. Notably, in addition to the scaling with the sample volume and the obvious dependence on instrumental noise, we find that the initial polarisation $p$ of the spin ensemble plays a prominent role in our projected sensitivities: scaling as $\delta_{\nu}\propto p^{-2}$, a competitive constraint requires at least $p\sim 0.25$. This requirement also means that thermal polarisation levels achievable at typical experimental conditions—$p\sim 10^{-3}$ for 129Xe at $T\sim 4$ K and $B\sim 10$ T—are not suited to C$\nu$B detection. Advanced polarisition techniques are needed to reach the desired $p\gtrsim 0.25$.

In optimistic experimental scenarios—with a sample size up to $R=10$ cm, and 100% polarisation—the constraint may improve to $\delta_{\nu}\sim 10^{11}$ for a small range of neutrino mass values. This number is largely consistent with the estimate of Ref. Arvanitaki:2024taq , and also at face value comparable to the current KATRIN limit on the C$\nu$B overdensity KATRIN:2022kkv .

Our results indicate that C$\nu$B direct detection is likely beyond the reach of current NMR technology, even in the most optimistic scenario. Indeed, even meeting the “realistic” projected $\delta_{\nu}$ constraints represent significant experimental challenges. We emphasise however much effort is currently being devoted to overcoming these challenges—in the case of CASPEr, in the context of dark matter axion searches. Should the desired experimental configurations materialise, we have shown that, by means of a pulsed-NMR measurement such as performed recently in Walter:2025ldb , the potential exists for such axion searches to simultaneously set potentially competitive constraints on the C$\nu$B for free. Constraining the C$\nu$B with nuclear spin experiments therefore remains an attractive possibility.

Acknowledgements.

We thank Caterina Braggio, Giovanni Carugno, Jan Hamann, Alex Sushkov, and Oleg Sushkov for enlightening discussions, and Marios Galanis for insightful comments on the manuscript. YG, GP, and Y3W are supported by the Australian Research Council’s Discovery Project (DP240103130) scheme.

Data Availability Statement

The data that support the findings of this article are openly available at https://github.com/gpierobon/OpenNu.

Appendix A Low energy non-relativistic spin Hamiltonian

We outline the derivation of the non-relativistic Hamiltonian for the neutrino-fermion spin interaction. We begin with the Lagrangian for a fermion $\Psi$ in the presence of an external electromagnetic 4-potential $A_{\mu}=(A_{0},\vec{A})$ and a 4-fermi interaction with a neutrino $\nu_{i}$ in the mass basis,

[TABLE]

where $i,j$ label the neutrino mass eigenstate, and $q$ is the electric charge of the fermion. The corresponding equation of motion is

[TABLE]

or, equivalently,

[TABLE]

where we have used $\vec{p}=-{\rm i}\vec{\nabla}$, and $p^{0}={\rm i}\partial_{t}$, and defined $\Sigma^{i,j}_{\mu}\equiv\bar{\nu}_{i}\gamma_{\mu}(1-\gamma^{5})\nu_{j}$.

We write the fermion spinor in 2-component notation,

[TABLE]

such that in the non-relativistic limit, we can make the ansatz

[TABLE]

which yields the solution $\chi_{0}\approx[\vec{\sigma}\cdot(\vec{p}-q\vec{A})/(2m)]\phi_{0}$, where $m$ is the mass of the fermion. Substituting this solution back into the equation of motion (73) in 2-component form, we find a Schrödinger-Pauli-like equation with additional neutrino-fermion interaction terms,

[TABLE]

where $\vec{B}=\vec{\nabla}\times\vec{A}$ is the external magnetic field, and $p^{0}=E$ is the energy.

In the form (76), the total Hamiltonian can be immediately identified with the terms in parentheses on the left-hand side of the equation. We concentrate exclusively on the spin interaction of the fermion with the C$\nu$B. Assuming the fermion (spin) is at rest, this interaction is described by the spin Hamiltonian,

[TABLE]

where we have used $\vec{J}=\vec{\sigma}/{2}$, and the magnetic moment is defined as $\vec{\mu}=(q/m)\vec{J}$. Although this result has been obtained assuming a charged fermion, an analogous spin-magnetic field interaction is also present for neutral particles like the neutron or nuclei. Here, the interaction of a magnetic field with the magnetic moment due to the particle’s internal structure can be described by an effective Lagrangian of the form $\bar{\psi}\sigma^{\mu\nu}\psi F_{\mu\nu}$. In the non-relativistic limit, this interaction term also reduces to a Hamiltonian of the form $H=-\vec{\mu}\cdot\vec{B}$. For nuclei, the convention is to define $\vec{\mu}=g\mu_{N}\vec{J}$, where $g$ is the $g$-factor of the nucleus under consideration, and $\mu_{N}=e/2m_{p}$ is the nuclear magneton.

Supposing without any loss of generality that the magnetic field points along the negative $z$-axis, i.e., $\vec{B}=-B\hat{z}$, the spin Hamiltonian can be recast into a more suggestive form,

[TABLE]

where the Larmor frequency $\omega_{0}=\gamma B$ is also the energy splitting between two adjacent spin states, with $\gamma$ the gyromagnetic ratio, and we have defined $J_{\pm}\equiv J_{1}\pm{\rm i}J_{2}$, and

[TABLE]

Then, splitting the Hamiltonian into a spin ensemble and an interaction term as

[TABLE]

we identify $H_{S}\equiv\omega_{0}J_{z}$, $H_{I}\equiv\Sigma_{z}J_{z}+\Sigma_{-}J_{-}+\Sigma_{+}J_{+}$, and

[TABLE]

gives the interaction Hamiltonian in the interaction picture.

In the case of a system with $N$ spins, the Hamiltonians of the spin ensemble and the interaction can be straightforwardly generalised from Eqs. (80) and (81) to

[TABLE]

where the operators $J^{\alpha}_{\pm,z}$ act only on the spin $\alpha$ at the spatial position $\vec{x}_{\alpha}$, and $\Sigma^{\alpha}_{\pm,z}$ are defined as per Eq. (79) but with the understanding that the neutrino wavefunctions must now be specified as $\nu_{i}(t,\vec{x}_{\alpha})$.

Appendix B Excitation and de-excitation rates

B.1 Single spin

Consider first a single spin immersed in the cosmic neutrino bath. Following Fermi’s golden rule the excitation and de-excitation rates are given by

[TABLE]

where the superscripts “$+$” and “$-$” refer respectively to excitation and de-excitation, $|T^{+}_{\rm fi}|^{2}$ is the transition probability from the initial state $|{\rm i}\rangle=|\downarrow\rangle\otimes|{\rm i}_{\nu}\rangle$ to the final state $|{\rm f}\rangle=|\uparrow\rangle\otimes|{\rm f}_{\nu}\rangle$, while $|T^{-}_{\rm fi}|^{2}$ is its counterpart for the transition $|{\rm i}\rangle=|\uparrow\rangle\otimes|{\rm i}_{\nu}\rangle$ to $|{\rm f}\rangle=|\downarrow\rangle\otimes|{\rm f}_{\nu}\rangle$. The quantity $\rho(E_{\rm f})$ is the density of states at the final energy $E_{\rm f}$, which in our case can be written as,

[TABLE]

where $\vec{p}$ and $\vec{p}\,^{\prime}$ are the 3-momenta of the incoming and outgoing neutrinos respectively, $\vec{n}$ is the momentum of the final-state fermion (assuming it is initially stationary), $P_{{\rm i}_{\nu}}$ and $P_{{\rm f}_{\nu}}$ represent respectively the thermal weight of the initial-state neutrino and any Pauli blocking factor associated with the final state, $j$ and $k$ are mass indices of the initial- and final-state neutrinos, while $a$ and $c$ index their helicities. Since the recoil of the fermion is expected to be tiny, the angle of the outgoing neutrino is essentially unrestricted. We therefore proceed to use the Dirac delta in Eq. (84) to eliminate the integral over the nuclear momentum, to arrive at

[TABLE]

where it is now understood that the relation between $|\vec{p}|$ and $|\vec{p}\,^{\prime}|$ is determined by energy conservation alone.

At first order in perturbation theory (i.e., the Born approximation), the transition amplitude is given by

[TABLE]

Focussing on excitation, consider an initial-state neutrino specified by a wavefunction

[TABLE]

where $a$ is a helicity index. Suppose it scatters off a spin inelastically, thereby transferring energy $\omega_{0}$ to it, and then exits the interaction in a final-state plane wave given by $\nu_{k}(t,\vec{x})=u_{k,c}(\vec{p}\,^{\prime})e^{-{\rm i}(p^{\prime 0}t-\vec{p}^{\prime}\cdot\vec{x})}$. This process is represented by the leftmost diagram in Fig. 12, and its transition probability can be written as

[TABLE]

where we have used $|\langle\uparrow|J_{+}|\downarrow\rangle|^{2}=1$ at the second equality, and the factors of $1/2p^{0}$ originates from the normalisation of the $u$-spinors to $2p^{0}$. The thermal weight factors in Eq. (84) are then $P_{{\rm i}_{\nu}}=N_{j,a}(p^{0})$ and $P_{{\rm f}_{\nu}}=1-N_{k,a}(p^{\prime 0})$, where $N_{i,a}(p^{0})$ is the neutrino occupation number, while the Dirac delta $\delta^{(1)}(E_{\rm f}-E_{\rm i})$ becomes $\delta^{(1)}(\omega_{0}+p^{\prime 0}-p^{0})$, enforcing the change in energy between the initial- and final-state neutrinos to be equal to $\omega_{0}$, the energy splitting due to the external magnetic field. A similar calculation can also be performed in the case of an antineutrino exciting a spin, as represented by the middle diagram in Fig. 12.

Putting it altogether, we then find the excitation rate of one spin to be

[TABLE]

where we have defined the phase space measure ${\rm d}\tilde{p}\equiv{\rm d}^{3}\vec{p}/((2\pi)^{3}2p^{0})$. Ignoring the relative motion between the spin and the C$\nu$B, and assuming the neutrino occupation numbers to be independent of helicity, the expression simplifies to

[TABLE]

which, using the Dirac delta, can be further reduced to

[TABLE]

with $p_{k+}=[(\sqrt{|\vec{p}|^{2}+m_{j}^{2}}-\omega_{0})^{2}-m_{k}^{2}]^{1/2}$, and $p^{0}_{k+}=\sqrt{m_{k}^{2}+p^{2}_{k+}}$.

Computation of the de-excitation rate from inelastic neutrino-spin scattering proceeds in a similar fashion, and yields

[TABLE]

Again, ignoring any neutrino-spin relative motion and assuming helicity-independent neutrino occupation numbers, we find

[TABLE]

with $p_{k-}=[(\sqrt{|\vec{p}|^{2}+m_{j}^{2}}+\omega_{0})^{2}-m_{k}^{2}]^{1/2}$, and $p^{0}_{k-}=\sqrt{m_{k}^{2}+p^{2}_{k-}}$. The lower integration holds if it evaluates to a real number. Otherwise, it needs to be set to zero.

Suppose for simplicity $m_{j}=m_{k}=m$ and $\omega_{0}\ll m$. The lower integration limits of the $\gamma_{+}$ and $\gamma_{-}$ integrals, Eqs. (91) and (93), are $\sqrt{(\omega_{0}+m)^{2}-m^{2}}\approx\sqrt{2m\omega_{0}}$ and zero, respectively. We can take the former to zero too if $\omega_{0}\ll|\vec{p}|^{2}/(2m)\ll m$ is satisfied by the bulk of the neutrino distribution. The same condition also allows us to approximate $p_{k\pm}\approx\sqrt{|\vec{p}|^{2}\mp 2m\omega_{0}}\approx|\vec{p}|\mp m\omega_{0}/|\vec{p}|$. Then, to leading order, we find

[TABLE]

and

[TABLE]

where we have omitted for brevity the Pauli-blocking terms in both expressions. Thus, we see immediately that $\gamma_{+}<\gamma_{-}$, and $\gamma_{\rm net}$ is suppressed by a factor $\sim 2m\omega_{0}/|\vec{p}|^{2}$ relative to $\gamma_{\pm}$.

B.2 $N$ spins

Consider now a system of $N$ spins described by the interaction Hamiltonian (82). We are interested in the excitation and de-excitation rates of the whole spin ensemble, which can be defined via

[TABLE]

where

[TABLE]

are the excitation and de-excitation amplitudes of the $N$-spin ensemble from an initial state $|{\rm i}_{\rm spin}\rangle$ to an unspecified final state $|{\rm f}_{\rm spin}\rangle$. The summation over ${\rm f}_{\rm spin}$ encompasses a complete basis of final states; the ladder operators will project out the final states that have one unit of energy $\omega_{0}$ more ($+$) or less ($-$) than the initial state.

Again, focussing first on spin excitation due to inelastic scattering of a plane-wave neutrino from $\nu_{j}(t,\vec{x}_{\alpha})=u_{j,a}(\vec{p})e^{-{\rm i}(p^{0}t-\vec{p}\cdot\vec{x}_{\alpha})}$ to $\nu_{k}(t,\vec{x}_{\beta})=u_{k,c}(\vec{p}\,^{\prime})e^{-{\rm i}(p^{\prime 0}t-\vec{p}\,^{\prime}\cdot\vec{x}_{\beta})}$, we find

[TABLE]

A similar expression can be found for $|T^{-}_{\rm fi}|^{2}$. Then, following the same procedure as for a single spin, it is straightforward to arrive at the collective excitation and de-excitation rates

[TABLE]

where $\mu\equiv\hat{p}^{\prime}\cdot\hat{p}=\cos\theta_{\vec{p}\vec{p}\,^{\prime}}$, with $\theta_{\vec{p}\vec{p}\,^{\prime}}$ the angle between the incoming and outgoing momenta, and

[TABLE]

are form factors that depend on the initial state of the whole $N$ spin ensemble. Since the operator $\sum_{\alpha}J_{\pm}^{\alpha}$ already projects the initial state into the subspace of final states that has the correct energy requirements, we can use the completeness relation of all possible final states $\sum_{\rm f^{\prime}_{spin}}|{\rm f^{\prime}}_{\rm spin}\rangle\langle{\rm f^{\prime}}_{\rm spin}|=1$ to simplify the form factors to

[TABLE]

where it is understood that $F_{\pm}(\vec{p}\,^{\prime}-\vec{p})$ are to be evaluated for $\vec{p}\,^{\prime}=\vec{p}_{k+}$ or $\vec{p}\,^{\prime}=\vec{p}_{k-}$.

Suppose the spin ensemble is initially in the ground state, i.e., $|{\rm i}_{\rm spin}\rangle=|\downarrow_{1},\downarrow_{2},\cdots,\downarrow_{\gamma},\cdots,\downarrow_{N}\rangle$. We see immediately that the de-excitation form factor must always yield $F^{G}(\vec{p}\,^{\prime}-\vec{p})=0$, while the excitation form factor evaluates to

[TABLE]

irrespective of the momentum transfer. Thus, in this case there is no coherent effect from the neutrino-spin interaction. If, on the other hand, the spin ensemble is initially in the product state $|{\rm i}_{\rm spin}\rangle=|\rightarrow_{1},\rightarrow_{2},\cdots,\rightarrow_{\gamma},\cdots,\rightarrow_{N}\rangle$, where $|\rightarrow\rangle\equiv(|\uparrow\rangle+|\downarrow\rangle)/\sqrt{2}$, then the excitation and de-excitation form factors evaluate to

[TABLE]

which features possible coherent effects in the second term. When the spatial separation between spins $|\vec{x}_{\alpha}-\vec{x}_{\beta}|$ is much smaller than $V^{1/3}$, where $V$ is the volume of the sample, we can take the continuum limit and, for a constant spin density $n_{s}=N/V$, approximate the last expression as

[TABLE]

where $j_{1}(x)$ is the spherical Bessel function of order one, and we have assumed at the second equality the spin sample to be a sphere of radius $R$. Clearly, the second, $N^{2}$ term dominates if $|\vec{p}-\vec{p}\,^{\prime}|R\ll 1$, with $|j_{1}(x)|^{2}/x^{2}\to 1/9$ as $x\to 0$. In this limit, the system exhibits coherent behaviours.

Given the typical momentum of the neutrino bath $\vec{p}_{\nu}$, in order for the coherence condition to hold over the whole length $R$ of the spin sample, we must demand $|\vec{p}_{\nu}-\vec{p}\,^{\prime}(\vec{p}_{\nu})|R\ll 1$. Assuming for simplicity $m_{i}=m_{j}=m$, the typical momentum transfers during excitation and de-excitation are

[TABLE]

with $p_{\nu}\equiv|\vec{p}_{\nu}|$. If we want coherence over the whole object, we must demand that $|\vec{p}_{\nu}-\vec{p}\,^{\prime}(\vec{p}_{\nu})|^{2}\ll 1/R^{2}$ even for $\cos\theta_{\vec{p}\vec{p}\,^{\prime}}=-1$. For non-relativistic neutrinos, where $\sqrt{p_{\nu}^{2}+m^{2}}\approx m+p_{\nu}^{2}/(2m)$, this yields the condition

[TABLE]

or, equivalently,

[TABLE]

with $\epsilon_{\pm}\equiv p_{\nu}^{2}/(2m)\mp\omega_{0}$. This in turn leads to

[TABLE]

when expanded to linear order in $\epsilon_{\pm}$.

From Eq. (108), we see that excitation can only happen if $\omega_{0}<p_{\nu}^{2}/(2m)$. But, if excitation happens at all, the coherence condition is also automatically fulfilled for any incoming momentum satisfying $p_{\nu}R\ll 1$. On the other hand, while de-excitation has no kinematic threshold, coherence also does not come for free. Rather, in order to have coherent de-excitation, the condition $\omega_{0}\ll p_{\nu}^{2}/(2m)$ must be satisfied if $p_{\nu}R\ll 1$. Physically, this means that the energy splitting must be much smaller than the incoming neutrino kinetic energy. If this is not the case, the momentum transfer of the neutrino is large and coherence is lost. For a typical C$\nu$B momentum of $p_{\nu}\sim 5.3\times 10^{-4}~{\rm eV}\sim 27~{\rm cm}^{-1}$ and a neutrino mass of $m=0.05$ eV, this condition amounts to demanding that $\omega_{0}\ll 3\times 10^{-6}$ eV.

For typical spin samples of size $R\sim 1$ cm, the condition $p_{\nu}R\ll 1$ cannot be realised by the C$\nu$B across the whole sample volume. However, even in the case of $p_{\nu}R\gg 1$, coherence is still possible under limited conditions. Specifically, assuming $\omega_{0}\ll p_{\nu}^{2}/(2m)$, Eq. (105) evaluates to

[TABLE]

Demanding $|\vec{p}_{\nu}-\vec{p}\,^{\prime}(\vec{p}_{\nu})|^{2}R^{2}\ll 1$, we see that coherence can be achieved even in the case $p_{\nu}R\gg 1$ for small angles $\theta_{\vec{p}\vec{p}\,^{\prime}}\ll 1/(p_{\nu}R)$ between the incoming and outgoing neutrino momenta. Under this condition, we can approximate the coherent part of the form factors as

[TABLE]

In other words, the collective excitation and de-excitation rates $\Gamma_{\pm}$ in the $p_{\nu}R\gg 1$ regime are suppressed by a factor $(2p_{\nu}R)^{-2}$ relative to the fully coherent rates. Evaluating Eq. (105) for the momentum transfer squared $q^{2}\equiv|\vec{p}_{\nu}-\vec{p}\,^{\prime}(\vec{p}_{\nu})|^{2}$ in the forward direction (i.e., $\cos{\theta_{\vec{p}\vec{p}\,^{\prime}}}=1$) and expanding in $\epsilon_{\pm}$ as above, it is straightforward to show that the momentum transfer in this regime is $q\sim\omega_{0}m/p_{\nu}$.

Appendix C Derivation of the Lindblad equation

The Lindblad formalism assumes unitary dynamics of the total neutrino-spin system described by a Hamiltonian that can be decomposed into three parts,

[TABLE]

Here, $H_{S}$ is the free Hamiltonian of the spin ensemble, $H_{B}$ the free Hamiltonian of the neutrino bath, and $H_{I}$ specifies the interaction between the spin ensemble and the neutrino bath. The density matrix $\rho$ acts on the Hilbert space of the whole neutrino-spin system, and we assume it can be split into a direct product of the spin-system density matrix $\rho_{S}$ and the bath density matrix $\rho_{B}$, i.e.,

[TABLE]

Such a separation can be made in the limit of weak system-bath interactions and large baths (relative to the spin ensemble), where no significant correlations develop between the spin ensemble and the bath. This is the so-called Born approximation.

The dynamics of the total neutrino-spin system is governed by the Liouville-von-Neumann equation,

[TABLE]

where the density matrix $\rho$ is understood to be in the Schrödinger picture. In the interaction picture this equation can be equivalently written as

[TABLE]

where $H^{I}_{I}$ is given by Eq. (81), and has the formal solution

[TABLE]

Plugging (115) into (114) and doing a partial trace over the bath, we arrive at an equation of motion for the spin ensemble alone that reads

[TABLE]

where we have omitted the interaction label (superscript “$I$”) for notational simplicity.

Observe that, implicit in the interaction Hamiltonian $H_{I}$, the first term is formally $\mathcal{O}(G_{F})$, while the second term is $\mathcal{O}(G_{F}^{2})$, and repeatedly plugging Eq. (115) into (116) would yield in principle a perturbative series in $G_{F}$. Considering however that the interactions of interest are very weak, it suffices to iterate only once to capture the leading-order effects. Physically, the first, $\mathcal{O}(G_{F})$ term in Eq. (116) corresponds to the so-called Stodolsky effect Stodolsky:1974aq , which can be most easily seen in the Schrödinger picture where it manifests as a shift in the energy levels of the spin ensemble, provided that a non-zero neutrino-antineutrino asymmetry is present. We shall therefore not be concerned with this shift and will henceforth drop the $\mathcal{O}(G_{F})$ term.

To discern the role of the second, $\mathcal{O}(G_{F}^{2})$ term in Eq. (116), we apply the Markovian approximation, which is valid under two assumptions: (i) The bath’s correlation time is much shorter than the spin ensemble’s relaxation time. This ensures the bath quickly loses memory of the interaction, thereby preventing information from flowing back to the spin ensemble; and (ii) the bath does not change much over time. This is the case for the relatively weak interactions under consideration. Implementing the Markovian approximation consists formally in changing $\rho(s)\rightarrow\rho(t)$ on the right-hand side of Eq. (116), i.e.,

[TABLE]

such that any step taken by the spin ensemble at any time $t$ depends only on the state of the system at that time. Equation (117) is the so-called Redfield equation. Next, we define a new integration variable $s^{\prime}=t-s$ and take the integration limit to $\infty$ (or, equivalently in terms of the original integration variable, we take the lower limit to $s\to-\infty$). This approximation leads us to

[TABLE]

and corresponds to assuming that the correlations decay fast enough that the dynamics of the system remains the same irrespective of what happened at earlier times.

At this point we can simply evaluate Eq. (118) for the Hamiltonian (81). A brute-force evaluation will yield an intractably long expression. However, because the neutrino-spin interaction rate is expected to be much smaller than the Larmor frequency $\omega_{0}=g\mu_{N}B$ (i.e., many oscillations can occur over a time scale dictated by the interaction rate), we may adopt the so-called rotating-wave/secular approximation, which amounts to neglecting fast-oscillation terms of the form $e^{\pm{\rm i}\omega_{0}t}$ or $e^{\pm 2{\rm i}\omega_{0}t}$. Doing so yields

[TABLE]

where

[TABLE]

Note that $\Sigma_{\alpha}$ is an operator in the interaction picture and hence depends on time.

Rearranging Eq. (119) to separate the hermitian and non-hermitian terms yields the Lindblad equation for our system:

[TABLE]

with

[TABLE]

and

[TABLE]

and we have assumed that the bath is constant in time. Observe that the first, Hermitian term in Eq. (121) commutes with $H_{S}$. Its form therefore remains unchanged in Schrödinger picture and represents an energy shift in the free Hamiltonian that is suppressed by $G_{F}^{2}$ through the bath-dependent coefficients $G_{\pm,z}$. The remaining terms are non-hermitian and characterise dissipative effects in the spin ensemble due to the bath via the $\gamma_{\pm,z}$ coefficients. Note that if the $g$-factor was negative, e.g., in the case of electrons, the Larmor frequency $\omega_{0}$ would be positive, which would effectively exchange $\gamma_{+}$ and $\gamma_{-}$ in the Lindblad equation.

Appendix D Computation of $\gamma_{\pm},\gamma_{z}$

Our objective now is to compute $\gamma_{+}$, $\gamma_{-}$ and $\gamma_{z}$, which requires that we calculate correlators of $\Sigma_{\alpha}$ that depend on the neutrino fields. To this end, we assume that neutrinos are Dirac fermions, whose expansion in modes is given by

[TABLE]

where ${\rm d}\tilde{p}\equiv{\rm d}^{3}\vec{p}/((2\pi)^{3}2p^{0})$, $u_{s}(\vec{p})$ and $v_{s}(\vec{p})$ are the spinor solutions of the Dirac equation, and the creation/annihilation operators obey the anticommutation relations

[TABLE]

All other anticommutators of $a^{s}_{i,\vec{p}},a^{r{\dagger}}_{j,\vec{q}},b^{s}_{i,\vec{p}},b^{r{\dagger}}_{j,\vec{q}}$ are zero.

We start with $\gamma_{+}$ given in Eq. (122), which evaluates using (79) to

[TABLE]

where have introduced the notations $P_{L}\equiv(1-\gamma^{5})/2$ and $\Upsilon_{\pm}\equiv\gamma^{1}\pm{\rm i}\gamma^{2}$. Substituting Eq. (124) into this expression, the non-vanishing terms read

[TABLE]

To compute the thermal traces of the bath, we assume that each neutrino species assumes a thermal distribution, such that the background density matrix is given by

[TABLE]

where $T=1/\beta$ is the neutrino temperature, and the bath Hamilton can be expanded as

[TABLE]

Using this information and the anticommutation relations (125), the traces evaluate to

[TABLE]

where we have introduced the notation $\delta_{pq}\equiv(2\pi)^{3}2p^{0}\delta^{(3)}(\vec{p}-\vec{q})\delta_{rs}\delta_{ij}$, and $N_{i,a}(p^{0})$ ($\bar{N}_{i,a}(p^{0})$) represents the occupation number of the $i$th neutrino (antineutrino) mass eigenstate with helicity $a$.

Substituting the thermal traces into Eq. (127) and integrating over $s$, we find

[TABLE]

where we have assumed the neutrino occupation numbers to be independent of helicity. Compared with the single-spin excitation rate (89) computed in Appendix B, we recognise immediately that the first and third terms correspond to inelastic scattering of a neutrino and an antineutrino, respectively, off a spin, thereby transferring energy to it. The second term describes the absorption of a neutrino-antineutrino pair, which we did not include in the heuristic calculation of Appendix B, but can also excite a spin if kinematically allowed; the rightmost diagram in Fig 12 represents this process.

Neglecting the relative movement between the Earth and the C$\nu$B and using the properties of the spinors, Eq. (131) can be further simplified to

[TABLE]

Using the Dirac deltas to eliminate one integral, the final expression reads

[TABLE]

with $q_{k+}=[(\sqrt{|\vec{p}|^{2}+m_{i}^{2}}-\omega_{0})^{2}-m_{k}^{2}]^{1/2}$, and $q_{k+}^{0}=\sqrt{m_{k}^{2}+q_{k+}^{2}}$. Except for the appearance of the second integral, this expression is the same as the heuristic rate (91) estimated earlier. Observe also that the second integral contributes only if $m_{k}+m_{i}<\omega_{0}$. This is the kinematic requirement for neutrino-antineutrino pair absorption alluded to above. For typical values of the magnetic field considered in this work, $\omega_{0}\sim 10^{-9}\rm eV$. Thus this process plays a role only for extremely small neutrino masses $m_{i}\lesssim 10^{-9}\rm eV$.

The computations of $\gamma_{-}$ and $\gamma_{z}$ proceed in an analogous fashion. Contrary to $\gamma_{+}$, however, $\gamma_{-}$ describes de-excitation of a spin via inelastic scattering with a neutrino or an antineutrino, i.e., the spin loses energy. In addition, a spin (up) can de-excite (to a spin-down) by emitting a neutrino-antineutrino pair if kinematically permitted. The final expression for $\gamma_{-}$ reads

[TABLE]

with $q_{k-}=[(\sqrt{|\vec{p}|^{2}+m_{i}^{2}}+\omega_{0})^{2}-m_{k}^{2}]^{1/2}$, $q_{k-}^{0}=\sqrt{m_{k}^{2}+q_{-}^{2}}$, $q_{k^{\prime}}=[(\sqrt{|\vec{p}|^{2}+m_{i}^{2}}-\omega_{0})^{2}-m_{k}^{2}]^{1/2}$, and $q_{k}^{\prime 0}=\sqrt{m_{k}^{2}+q_{k}^{\prime 2}}$. Again, this expression can be contrasted with the heuristic rate (93), which does not contain the second integral that contributes only if $m_{k}+m_{i}<\omega_{0}$.

For $\gamma_{z}$ we find

[TABLE]

Contrary to $\gamma_{\pm}$, this term does not describe energy exchange processes. Rather, it describes the collective dephasing of the spin ensemble due to elastic scattering.

In the case of Majorana neutrinos, the main difference in the computation of $\gamma_{\pm,z}$ is that the expansion in modes of the neutrino fields is given by

[TABLE]

and the corresponding bath Hamiltonian is

[TABLE]

If we neglect pair absorption and emission (which, as discussed above, are normally not kinematically allowed in NMR experiments except for close-to-zero neutrino masses), the expressions for $\gamma_{\pm,z}$—Eqs. (133), (134), and (135)—are exactly identical for Dirac and Majorana neutrinos. The only difference between the two cases occurs in the value of the occupation number $N_{i,a}(p^{0})$. In the Dirac case, when ultra-relativistic neutrinos decouple from the early universe plasma, they are always in a left-helicity state. However, as the universe cools down and the neutrinos become non-relativistic, we would generically expect them to relax via gravitational interactions into an approximately equal number of left- and right-helicity states with an average occupation number (assuming no neutrino-antineutrino asymmetry) given by $N_{i,a}(p^{0})=\bar{N}_{i,a}(p^{0})\approx n_{F}(|\vec{p}|)/2$, where

[TABLE]

is the relativistic Fermi-Dirac distribution. On the other hand, Majorana neutrinos are their own antiparticles. Thus, relaxation into left- and right-helicity states merely turn an neutrino into an antineutrino and vice versa, leading to an occupation number of $N_{i,a}(p^{0})=\bar{N}_{i,a}(p^{0})\approx n_{F}(|\vec{p}|)$.

Appendix E Lindblad equation for $N$ spins

Consider now a system comprising $N$ $1/2$-spins, where the Hilbert space is spanned by $\rho_{S}=\otimes^{N}\rho_{1/2}$, and the Hamiltonian of neutrino-spin system is given in Eq. (82). The corresponding Lindblad equation can be generalised from Eq. (121). In the interaction picture (and omitting the hermitian terms) it reads

[TABLE]

with

[TABLE]

which are generalisations of the $\gamma_{\pm,z}$ coeffcients from Eq. (122).

Consider first the coefficient $\gamma^{\alpha\beta}_{+}$. Following the procedure outlined in Appendix D, we find

[TABLE]

After evaluating the thermal traces, the expression simplifies to

[TABLE]

where $\Delta\vec{x}_{\alpha\beta}\equiv\vec{x}_{\alpha}-\vec{x}_{\beta}$.

From Eq. (142) we see immediately that in the limit $|\vec{p}-\vec{q}|\ll|1/\Delta\vec{x}_{\alpha\beta}|$ and $|\vec{p}+\vec{q}|\ll|1/\Delta\vec{x}_{\alpha\beta}|$, the exponentials evaluate to unity such that $\gamma^{\alpha\beta}_{+}$ loses its dependence on $\alpha$ and $\beta$, and approximates $\gamma_{+}$ of Eq. (131) computed previously for the single-spin system. The same limiting behaviour is also seen in $\gamma_{-}^{\alpha\beta}$ and $\gamma_{z}^{\alpha\beta}$ (not shown here). Then, in this coherent regime, the Lindblad equation reads

[TABLE]

where we have defined ${\cal J}_{+}\equiv\sum_{\alpha}J^{\alpha}_{+}$, ${\cal J}_{-}\equiv\sum_{\alpha}J^{\alpha}_{-}$, and ${\cal J}_{z}\equiv\sum_{\alpha}J^{\alpha}_{z}$. If on the other hand $|\vec{p}-\vec{q}|\gg|1/\Delta\vec{x}_{\alpha\beta}|$, the cross (i.e., $\alpha\neq\beta$) terms in Eq. (142) average to zero and only the $\alpha=\beta$ terms survive. (The condition on $|\vec{p}+\vec{q}|$ versus $|1/\Delta\vec{x}_{\alpha\beta}|$ is not important, as neutrino-antineutrino pair absorption and emission are not kinematically allowed for neutrino masses $m_{i}\gtrsim 10^{-9}\rm eV$ given typical values of the magnetic field $B\sim 0.1$ T.) In this limit the Lindblad equation reads

[TABLE]

a form that is generally applicable for any “local” interaction that is not coherent over the length of the spin sample.

Appendix F Full numerical solutions

We demonstrate here how we find a full numerical solution of the Lindblad master equation (43) for the spin ensemble’s density matrix $\rho_{S}$. We show first the case of collective neutrino interactions only. This effectively restricts the dynamics to the fully symmetric subspace of total angular momentum $j=N/2$, where numerical solutions can scale with $N^{3}$ (see Ref. PhysRevA.98.063815 for a detailed discussion) and are hence tractable for relatively large values of $N$. We then extend the numerical treatment to include local interactions, where the dynamics can no longer be restricted to the $j=N/2$ manifold.

In the presence of coherent interactions alone, we set $\gamma_{\phi}^{\rm loc}=0$, and solve

[TABLE]

where we have introduced the Liouvillian super-operator at the second equality, expressed in the Dicke basis. We work in Liouville space by vectorising (i.e., column-stacking) the density matrix $\rho_{S}\to|\rho_{S}\rangle$ and using the Liouvillian super-operator

[TABLE]

where the action of the operators on $\rho_{S}$ from the left and the right becomes a Kronecker product, namely, $X\rho_{S}Y\to(Y^{\top}\otimes X)|\rho_{S}\rangle$.

To numerically integrate the master equation (145), we implement standard time-integration schemes such as the Euler method, as well as higher-order solvers such as the fourth-order Runge–Kutta (RK4) algorithm. To ensure numerical stability, especially for such explicit integration schemes, we compute the spectral radius of the Liouvillian, i.e., the largest eigenvalue of $L_{D}$, $\lambda_{\rm max}$, and set the time step $\Delta t$ such that $\Delta t\,\lambda_{\rm max}\ll 1$. This criterion ensures convergence and bounds the error at each step, particularly for stiff systems or rapidly decaying modes. In practice, we choose $\Delta t=h/\lambda_{\rm max}$ with $h\sim 0.1$ depending on the desired accuracy. We have additionally tested implicit integration methods and found no substantial improvement in convergence compared to explicit methods, provided the timestep factor $h$ is sufficiently small. Given their significantly lower computational cost, explicit schemes remain the optimal choice for efficiency in this context.

Our simulations are implemented in the publicly available code OpenNu,151515The source code is available at https://github.com/gpierobon/OpenNu. written in C++ and based on the Eigen161616Eigen is a high-performance C++ library for linear algebra: https://eigen.tuxfamily.org library for efficient sparse matrix operations. This implementation achieves high computational performance and a low memory footprint, which is essential for solving large-scale Lindblad dynamics. In particular, it allows us to simulate on a single modern CPU collective spin ensembles with Hilbert space dimensions up to the equivalent of $N\sim 10^{4}$. To avoid numerical instabilities due to very large or small values, we normalise time by the characteristic dissipation timescale $t_{\rm diss}$, achieved by rescaling the Liouvillian with the product of the number of spins and the de-excitation rate $\gamma_{-}$. As a result, the ratio $\gamma_{+}/\gamma_{-}$ remains as the sole control parameter governing the dynamics in this regime. The left panel of Fig. 13 displays our numerical solutions obtained using ${\tt OpenNu}$ for $N=100$ and different initial conditions.

To include local interactions and fully solve Eq. (43), we use the Permutational-Invariant Quantum Solver (PIQS) introduced in Ref. Shammah_2017 and implemented in the qutip library.171717The qutip source code is available at github.com/qutip/qutip PIQS also reduces the computational complexity of simulating open quantum dynamics from exponential to polynomial in $N$ by exploiting the permutational symmetry of two-level systems, under the assumption that they are identical and identically prepared. However, when the dynamics involve population outside the fully symmetric subspace, i.e., when $j<N/2$ states are included, the symmetric Dicke basis must be extended to account for the action of local operators such as $J^{\alpha}_{z}$ in Eq. (43). As a result, with respect to OpenNu the size of the Liouvillian super-operator increases by a factor of $N^{2}$, significantly raising the computational cost and limiting tractable simulations to system sizes of $N\sim{\cal O}(100)$. In the right panel of Fig. 13, we present numerical solutions computed with ${\tt qutip}$ for values up to $N=50$, demonstrating excellent agreement between the full density matrix solution and the fast second-order approximation.

Appendix G Hyperpolarisation

Hyperpolarised spin systems, which are out of thermal equilibrium, enhance signal intensity in NMR experiments Eills2023SpinHyperpolarization . A technique widely used in NMR applications is spin-exchange optical pumping (SEOP) RevModPhys.69.629 , particularly for noble gases such as 129Xe or 3He. SEOP can be broken down into two steps:

• •

Optical pumping of alkali metal atoms. A circularly polarised laser excites alkali atoms (typically Rb or Cs vapour) in a low magnetic field. A large electron spin polarisation $p_{\rm A}\approx 1$ in the alkali ground state develops rapidly via repeated excitation-relaxation cycles and from angular momentum conservation in the atomic transition. The dynamics of this process is governed by a rate equation Walker_2011 ,

[TABLE]

where $\Gamma_{\rm pump}$ is the pumping rate, and $\Gamma_{{\rm A,rel}}$ the total alkali spin relaxation rate. The steady-state $p_{A}$ is given by

[TABLE]

which shows that to obtain a large $p_{\rm A}$ requires $\Gamma_{\rm pump}\gg\Gamma_{\rm A,rel}$.

• •

Spin-exchange collisions. The polarised alkali atoms collide with the target noble gas, e.g., Xe, atoms in a gas cell, effectively transferring angular momentum to the noble gas nuclear spins (via a Fermi-contact interaction). The polarisation of the noble gas spins $p_{\rm Xe}$ can be computed from

[TABLE]

where $\gamma_{\rm SE}$ is the spin-exchange rate, and $\Gamma_{\rm rel}$ is the nuclear spin relaxation rate. Then,

[TABLE]

gives the steady-state value of $p_{\rm Xe}$.

SEOP typically operates in dilute gas phases, where nuclear spin-spin interactions are negligible and polarisation accumulates through random spin-transfer collisions that do not generate coherences between spins. The overall mechanism is analogous to the local interaction channels described in Sec. V. Matching $p^{\rm ss}_{\rm Xe}$ of Eq. (150) to the polarisation parameter in Eq. (47) then yields a mapping

[TABLE]

between the various rates to the hyperpolarisation excitation and de-excitation rates used in Sec. VI.1.

Appendix H Two-time correlation functions

We outline here the computation of the two-time correlation functions $\langle C(t)A(t+\tau)\rangle$ via the quantum regression theorem QRT , where $\tau\geq 0$, and the operators $A$ and $C$ act on the spin ensemble alone (i.e., not the bath). See also Ref. carmichael1 for a more comprehensive introduction. This correlation function is defined in the Heisenberg picture via

[TABLE]

where $\rho(0)$ is the density matrix of total system of spin and bath at time $t=0$, and the trace ${\rm tr}_{T}$ pertains to the total system. Applying time evolution to the operators, i.e., $A(t)=e^{{\rm i}Ht}A(0)e^{-{\rm i}Ht}$, $C(t)=e^{{\rm i}Ht}C(0)e^{-{\rm i}Ht}$, and $\rho(t)=e^{-{\rm i}Ht}\rho(0)e^{{\rm i}Ht}$ assuming a constant $H$, allows us to rewrite the expression in the Schrödinger picture,

[TABLE]

where the trace is now carried out over the state at time $t+\tau$, with $A(0)$ and $C(0)$ as time-independent Schrödinger-picture operators.

Assuming the total density matrix can be split into a spin ensemble and a bath piece, i.e., $\rho(t)=\rho_{S}(t)\otimes\rho_{B}(t)$, at all times allows us to “swap the order” of $\rho_{B}(t)$ and $C(0)$ to get

[TABLE]

where we have defined in the last line a modified spin density matrix $\rho_{S;C}(t)\equiv\rho_{S}(t)C(0)$. Writing out the total trace ${\rm tr}_{T}$ in terms of individual traces over the bath (subscript $B$) and spin ensemble (no subscript, for consistency with the notation of the main text), i.e., ${\rm tr}_{T}={\rm tr}_{B}\otimes{\rm tr}$, we arrive at

[TABLE]

where

[TABLE]

is simply the modified spin matrix $\rho_{S;C}$ forwarded in time in time from $t$ to $t+\tau$ under the total Hamiltonian $H$ with the bath degrees of freedom traced out.

Using $\tau$ as the time variable, the formal equation of motion for $\rho_{S;C}(\tau)$ is

[TABLE]

where ${\cal L}$ is some super-operator, which in our case would be the Lindbladian. The formal solution reads $\rho_{S;C}(t+\tau)=e^{{\cal L}\tau}\{\rho_{S;C}(t)\}$, such that the two-time correlation function can also be written as

[TABLE]

Note that while we have assumed a time-independent $H$ in the above derivation, the equation of motion (157) applies also to a time-dependent $H$. The ${\tt qutip}$ library contains functions that compute numerically two-time correlation functions for a given ${\cal L}$.

We are interested to derive equations of motion directly for the two-time correlation functions $\langle C(t)A_{\mu}(t+\tau)\rangle$, $A_{\mu}$ is a set of operators indexed by $\nu=1,2,3,\ldots$. Suppose $\langle A_{\mu}\rangle$ evolves under the Lindblad master equation as

[TABLE]

where $M_{\mu\alpha}$ is square matrix of constant coefficients. It then following straightforwardly from Eqs. (157) that

[TABLE]

In other words, the “same” set of equations of motion describe both the single-time expectation values and the two-time correlations.