Neutrino thermalization via randomization on a quantum processor

TL;DR

This work tackles neutrino flavor thermalization in dense astrophysical environments by modeling the system with an all-to-all Heisenberg-like Hamiltonian featuring random couplings. It introduces randomized quantum circuits to emulate non-local dynamics with depth that remains effectively independent of system size, and measures thermalization via the flavor variance $V(t)$, identifying a thermalization time $\tau_{th}$ that scales as $O(\sqrt{N})$ in large systems. The authors validate the approach through classical simulations (tensor networks, PSA, Pauli propagation) and IBM quantum hardware with error mitigation, observing consistent sqrt($N$) scaling and favorable agreement with semi-classical PSA predictions. This demonstrates the utility of near-term quantum devices as empirical validators of classical methods for complex many-body dynamics and opens pathways for applying random circuits to challenging problems beyond classical reach.

Abstract

The dynamical evolution of neutrino flavor in supernovae can be modeled by an all-to-all spin Hamiltonian with random couplings. Simulating such two-local Hamiltonian dynamics remains a major challenge, as methods with controllable accuracy require circuit depths that increase at least linearly with system size, exceeding the capabilities of current quantum devices. The eigenstate thermalization hypothesis predicts that these systems should thermalize, a behavior confirmed in small-scale classical simulations. In this work, we investigate flavor thermalization in much larger systems using random quantum circuits as an empirical tool to emulate the non-local dynamics, and demonstrate that the thermal behavior can be reproduced using a depth independent of the system size. By simulating dynamics of over one hundred qubits, we find that the thermalization time grows approximately as the square root of the system size, consistent with predictions from semi-classical methods. Beyond this specific result, our study illustrates that near-term quantum devices are useful tools to test and validate empirical classical methods. It also highlights a new application of random circuits in physics, providing insight into complex many-body dynamics that are classically intractable.

Figures (7)

• Figure 1: Trotterized dynamics. Scaling of the characteristic time as function of the system size for the single qubit observables (red), Loschmidt echo (purple), single qubit entropy (orange) and variance (blue), alongside square root and logarithm fits (dotted and dashed lines). The calculation are performed with a converged second-order Trotter product formula ($\delta t = 0.5$).
• Figure 2: Rescaling factor. The rescaling factor normalized with $\sqrt{N/l}$ shown as a function of the number of layers for different system sizes. The expectation value and standard deviation is computed across Hamiltonian realization and initial state $(k,s)$ at fixed depth $l$.
• Figure 3: Protocol A. The variance of a single-qubit observable, denoted $V(t)$, is shown as a function of time under protocol A for different number of (fix) layers, time step $\tau = t/L$ and system sizes: (a) $N=10$, (b) $N=16$, (c) $N=22$. The black solid line represents the exact dynamics computed using a converged product formula, while the colored line corresponds to the evolution with the randomized protocol. The shaded band around the latter indicates the interquartile range across disorder realizations.
• Figure 4: Protocol B. The variance of a single-qubit observable, denoted $V(t)$, is shown as a function of time under protocol B for different number of layers $L$, time step $\tau=t_{\text{end}}/L$ and system sizes: (a) $N=10$, (b) $N=16$, (c) $N=22$. The black solid line represents the exact dynamics computed using a converged product formula, while the colored line corresponds to the evolution under the random channel. The shaded band around the latter indicates the interquartile range across disorder realizations.
• Figure 5: MPS simulations. The variance $V(t)$ is shown as a function of the number of layers using a fixed time step $\delta t = 1$ for different bond dimension and system sizes $N=20$ (a) and $N=32$ (b). The insets shows the relative absolute error with the exact evolution, where for the larger system we use the largest bond dimension as a proxy.