Quantum information and quantum simulation of neutrino physics

Neutrinos, owing to their feeble interactions with matter, are highly efficient at transporting energy, entropy, and lepton number in extreme astrophysical settings such as core-collapse supernovae and binary compact object mergers (binary neutron star or black hole - neutron star), as well as during certain epochs in the early universe (e.g., see Refs. Janka:2006fh ; Burrows:2020qrp ; Fuller:2022nbn ; Foucart:2022bth ; Kyutoku:2017voj ; Grohs:2015tfy ). As a result, they are expected to play a key role in influencing the dynamics and nucleosynthesis in these environments.

The flavor-specific interactions that govern the neutrino transport and determine the free neutron and proton abundances in these environments are charged current (anti-)neutrino captures that convert neutrons into protons, and vice versa. Given the typical temperatures and densities of these environments, neutrinos decouple with energies of $\mathcal{O}(1\text{--}10)\,$MeV, and therefore, the charged-current interactions of $\mu$ and $\tau$ flavor (anti-)neutrinos are energetically suppressed. Since these charged current processes govern the energy transport as well as the neutron-to-proton ratio (or equivalently electron fraction), which in turn determines nucleosynthesis yields (e.g.,  Surman:2003qt ; Martinez-Pinedo:2017ksl ; Kajino:2014bra ; Frohlich:2015spx ; Langanke:2019ggn ; Roberts:2016igt ; Steigman:2012ve ; Grohs:2015tfy ), the flavor-asymmetric nature of charged-current capture implies that a complete understanding of neutrino flavor evolution in these environments is crucial Qian:1993dg ; Yoshida:2006qz ; Duan:2010af ; Kajino:2012zz ; Wu:2014kaa ; Sasaki:2017jry ; Balantekin:2017bau ; Xiong:2019nvw ; Xiong:2020ntn ; George:2020veu .

In this paper we first summarize recent progress in our understanding of neutrino oscillations in extreme astrophysical environments, particularly the quantum many-body aspect of collective neutrino oscillation physics driven by $\nu$-$\nu$ interactions in high neutrino densities (Secs. 2–6). In Sec. 7, we show some new results of many-body neutrino calculations in a three-flavor setup. Finally, we offer concluding remarks in Sec. 8.

Neutrinos experience flavor oscillations thanks to a misalignment between the propagation (mass) eigenstates and the eigenstates of the weak interaction (flavor). This misalignment, present even in vacuum, can be modified in interesting ways when neutrinos pass through a medium. Where the medium consists of baryons and charged leptons, neutrino coherent forward-scattering with these particles can lead to suppressed or resonantly enhanced flavor conversion through the Mikheyev-Smirnov-Wolfenstein (MSW) mechanism Wolfenstein:1977ue ; Mikheyev:1985zog ; Mikheev:1986if . This mechanism is the widely accepted solution for explaining the observed solar neutrino deficit. Even more fascinating is the scenario where the neutrinos themselves form a significant component of the background, as can be the case in extreme astrophysical environments like core-collapse supernovae, binary compact object mergers, and the early universe. In this situation, neutrinos experience an additional forward-scattering potential on account of pairwise weak neutral current interactions with other neutrinos Notzold:1987ik ; Fuller:!987aa . In the low-energy limit, this can be represented by the four-Fermi effective operator

Unlike the interactions with baryons and charged leptons, which are diagonal in the flavor basis, the coherent $\nu$-$\nu$ scattering can give rise to both diagonal and off-diagonal potentials, each of which depend on the flavor composition of the neutrino background Pantaleone:1992xh ; Pantaleone:1992eq ; Samuel:1993 . In this sense, the background itself becomes dynamical and the flavor evolution histories of the interacting neutrinos thus become coupled to one another. This coupling introduces a non-linearity and a nontrivial geometrical complexity to the flavor evolution problem, giving rise to many kinds of collective oscillation modes in flavor space. See, for example, the reviews in Refs. Duan:2009cd ; Duan:2010 ; Chakraborty:2016yeg ; Tamborra:2020cul ; Richers:2022zug and references therein. In particular, a lot of the recent literature deals with the aptly named “fast flavor transformations”  which occur as a result of nontrivial angular distributions in the neutrino flavor field (reviewed in Refs. Chakraborty:2016yeg ; Tamborra:2020cul ; Richers:2022zug ). Flavor instabilities arising in this manner grow on a much faster timescale compared to the instabilities present in physical setups with a greater degree of angular symmetry, hence the moniker “fast”.

Another feature of this problem that has received attention in recent years is the quantum nature of these collective effects, as part of a broader focus on quantum simulations in high energy physics Bauer:2022hpo . The pairwise interactions between neutrinos that give rise to the coherent $\nu$-$\nu$ scattering potentials assume the form of a spin-spin coupling, when expressed using “isospin” operators defined in the Hilbert space of neutrino flavor states (described in detail in the following section). Depending on whether one considers two or three flavors, the isospin operators for each neutrino form a SU(2) or SU(3) algebra. The system of interacting neutrinos then constitutes a quantum many-body system, with the size of the associated Hilbert space scaling exponentially with the number of particles in the system ($2^{N}$ or $3^{N}$ for a system of $N$ neutrinos and anti-neutrinos in two or three flavors, respectively).

Owing to the immense computational difficulty of analyzing a many-body system with an appreciable number of neutrinos, it becomes necessary to invoke certain simplifying assumptions. One common approach involves postulating that the effect of multi-particle quantum correlations on the flavor evolution will remain negligible in the large-$N$ limit, as a result of a large number of random phases being added incoherently. With this assumption, one is allowed to construct a simplified Hamiltonian, where the two-particle interaction operator is replaced with an effective one-particle interaction with a “mean-field” representing all of the background particles Samuel:1993 ; Sigl:1993ctk ; Qian:1994wh . In this “Random phase approximation”, the effective dimensionality of the state-variable space is reduced to $n_{f}N$ (where $n_{f}$ is the number of flavors), thus making the problem much more tractable numerically.

However, this simplification naturally raises the question as to whether such an assumption can lead to the exclusion of any crucial physics from the problem. Several recent studies have attempted to answer this question using a variety of different physical setups and numerical methods. In what follows, we review some of the recent progress.

In much of the recent literature investigating many-body neutrino correlations, the neutrinos are modeled as interacting plane waves in a box of volume $V$ (which can be taken to be time-dependent to mimic the decrease of the neutrino number density with distance from the source). In this picture, the finite size of the neutrino wave-packet (and consequently, the finiteness of the interaction interval between any two neutrinos) is not taken into account. Therefore, this setup may not be fully suitable for analyzing certain features of this problem, such as effects of incoherent scattering on the neutrino flavor conversion Shalgar:2023ooi . Nevertheless, it can be a useful initial step for illustrating the feedback effect from non-trivial $\nu$-$\nu$ correlations on the flavor dynamics, in a simplified setting. In order to fully validate or invalidate the mean-field approximation, more detailed analyses including wave-packet effects may be needed in the future.

The Hamiltonian describing a system of interacting neutrinos can be written in terms of generators of $SU(n_{f})$ where $n_{f}$ are the number of neutrino flavors.

For instance, in a two-flavor model, the SU(2) generators, in terms of the Fermionic creation and annihilation operators, are defined as Balantekin_2006

in the mass basis. Alternatively, one may construct a flavor basis SU(2) algebra, with $a_{e}$ and $a_{x}$ replacing $a_{1}$ and $a_{2}$ (where $x$ represents an appropriate superposition of $\mu$ and $\tau$ flavors). In the spin-$1/2$ representation, one can write these operators in terms of Pauli matrices: i.e., $\vec{J}_{\mathbf{p}}=\vec{\sigma}_{\mathbf{p}}/2$, where $\vec{\sigma}_{\mathbf{p}}$ is a vector of Pauli matrices defined in the subspace of the neutrino with momentum $\mathbf{p}$. In what follows, we also refer to these SU($n_{f}$) generators as neutrino “isospin” operators. In the two-flavor context, isospin “up” and “down” can refer to $\ket{\nu_{1}}$ and $\ket{\nu_{2}}$, respectively, in the mass basis, or $\ket{\nu_{e}}$ and $\ket{\nu_{x}}$ in the flavor basis. In this two-flavor picture, a Hamiltonian consisting of terms that represent vacuum mixing as well as $\nu$-$\nu$ interactions¹¹1 Here, we exclude the term representing neutrino interactions with ordinary matter (e.g., baryons and charged leptons), since it has a structure that is conceptually similar to the vacuum oscillation term—i.e., consisting of individual neutrinos interacting with a background. In regimes where collective oscillation effects typically dominate, this matter-interaction term can be “rotated away” with a suitable change-of-basis transformation, resulting in a modified mixing angle and mass-squared splitting compared to the corresponding values in vacuum. can be written as

where $\vec{B}$ can be interpreted as a “background field” that points along the direction of the mass basis. It is equal to $(0,0,-1)$ in the mass-basis representation, or $(\sin 2\theta,0,-\cos 2\theta)$ in the flavor-basis representation. $\omega_{\mathbf{p}}={\delta m^{2}}/{(2|\mathbf{p}|)}$ are the vacuum oscillation frequencies for neutrinos with momenta $\mathbf{p}$, with $\delta m^{2}$ being the mass-squared difference between the eigenstates. $\mathbf{\widehat{p}}$ and $\mathbf{\widehat{q}}$ are the unit vectors along the momenta of interacting neutrino pairs, and $V$ is the quantization volume. One can define a $\nu$-$\nu$ coupling parameter $\mu\equiv\sqrt{2}G_{F}N/V$, where $N$ is the total number of interacting neutrinos. The strength of the $\nu$-$\nu$ interactions thus depends on the neutrino number density and the intersection angle between their trajectories of the interacting neutrinos. This dependence introduces an additional geometric complexity to the problem, aside from the complexity associated with the exponential scaling of the Hilbert space. Note that the Hamiltonian of Eq. (5) consists only of terms that either preserve or exchange the momenta of interacting neutrino pairs (forward and exchange terms), since these can be added up coherently in the mean-field limit. Recent work has argued that these terms should not enjoy special status in an exact many-body calculation, unlike in the mean field limit, and, as a result, consideration of a more generalized Hamiltonian with interaction terms besides forward and exchange scattering is warranted Johns:2023ewj . This addition could be an important avenue to pursue in future studies.

In the mean field approach, the interaction term in the Hamiltonian is replaced with an effective one-particle operator, using the following prescription:

This approximation in essence enforces that the wavefunctions of individual neutrinos remain uncorrelated through the course of the evolution (assuming that the neutrino ensemble starts out in an uncorrelated state), and the system as a whole thereby remains in a direct product state, greatly reducing the effective number of independent amplitudes from $n_{f}^{N}$ to $n_{f}N$.

To calculate quantum corrections beyond the mean-field coherent limit, one can systematically incorporate $n$-body density matrices $\rho_{1\ldots n}$ for $n\geq 1$, given by

into the coupled equations of motion for $N$ neutrinos, as follows Volpe:2013uxl :

where $H_{1\ldots n}$ is the Hamiltonian truncated for the first $n$ neutrinos in a given ordering and $V(i,j)$ is the two-body interaction potential for a pair of neutrinos $(i,j)$. This procedure follows the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy for density matrices: the mean-field equations can be recovered by restricting to $n=2$ and approximating $\rho_{12}\approx\rho_{1}\rho_{2}$ (i.e., requiring the two-body correlation function to be zero). In this description, investigating the importance of quantum corrections would involve incorporating the $n$-body density operators for progressively increasing values of $n$, while checking for convergence of results for physical observables.

Because of the exponential growth ($2^{n}$) of the size of the Hilbert space, classical computers are unable to exactly simulate systems of more than $\sim 30$ neutrinos. One possibility is to introduce compact representations of the wave function through tensor network methods Roggero:2021 ; Roggero:2021asb ; Cervia:2022 , and more specifically matrix product states Vidal:2003 ; 2011AnPhy.326…96S ; PAECKEL2019167998 . These methods allow for the computation of systems of hundreds of neutrinos in a two-beam setup with a time-independent $\nu$-$\nu$ interaction Roggero:2021 ; Roggero:2021asb . Alternatively, when considering very dense neutrino gases (if vacuum oscillations are ignored), methods based on generalized angular momentum representations, by analogy between two flavor oscillations and spin systems (see Eq.11), can reach up to $\mathcal{O}(10^{6})$ neutrinos and predict the thermodynamic limit in some cases Friedland:2003eh ; Friedland:2006ke ; Xiong:2021evk ; Martin:2021bri ; Roggero:2022 .

In the case of a time-dependent interaction strength, a more sophisticated tensor network method, namely, the time-dependent variational principle (TDVP) method has been utilized in Ref. Cervia:2022 . These techniques provided considerable computational benefit for an initial state with all neutrinos in the same flavor, allowing for evolution of a system with $\approx 50$ oscillation modes. This limit was a consequence of the entanglement among neutrinos being more localized in certain regions of the neutrino energy distribution. For systems with initial states being a mixture of $\nu_{e}$ and $\nu_{x}$ flavors, the entanglement is more de-localized, and therefore the comparative advantage gained through TDVP methods is less dramatic, although work remains in progress on this front. As shown in Fig. 2, for a system of eight neutrinos with an initial state $\ket{\psi}=\ket{\nu_{e}}^{\otimes 8}$ the survival probability for the neutrino in the highest frequency mode in the first mass eigenstate $P_{\nu_{1}}$, converges to the exact value obtained from RK4 for a very small bond dimension $D\sim 4$. Therefore, $D$ can be truncated to a significantly small number (for $N=8$ the maximum bond dimension is $2^{4}=16$) in the case of an initial state with all neutrinos in the same flavor. In the case of an initial state with half neutrinos in $\nu_{e}$ and $\nu_{x}$ flavor each, the $P_{\nu_{1}}$ does not converge to the exact value until $D\sim 15$. Furthermore, the results do not converge to those of RK4 exactly, and hence the truncation requires considering a smaller time step for more accurate results.

Non-vanishing connected correlations, being the essential characteristic of beyond mean-field behavior, deserve particular attention. Ursell functions Shlosman1986 provide a unified framework for $n$-connected correlations,

and have been studied in recent works Illa:2022zgu . It is worth pointing out that, despite only pairwise interactions in the system, higher order correlations dynamically develop as well.

For generic closed quantum systems, quantum simulation algorithms are promising tools to study quantum many-body evolution. Preliminary attempts Hall:2021rbv ; Yeter-Aydeniz:2021olz ; Illa:2022jqb ; Amitrano:2022yyn ; Illa:2022zgu to simulate collective neutrino oscillations on a quantum computer have already been taken, for small system sizes and short evolution times. In Ref. Hall:2021rbv , a system of four neutrinos was simulated using superconducting qubit hardware. In particular, the unitary evolution operator $U(t)=\exp(-\mathrm{i}Ht)$ was approximated via first-order Trotter-Suzuki decomposition, with error of $\mathcal{O}(t^{2})$.

Notably, since the interaction in this model is long-range, quantum devices with all-to-all connectivity are desirable. Nevertheless, on a quantum device having connectivity only among neighboring qubits, SWAP operations can still be used to implement this interaction Hall:2021rbv , though doing so requires more quantum gates in the simulation that may decohere the quantum state being simulated. Alternatively, hybrid quantum algorithms such as quantum Lanczos (QLanczos) could be used Yeter-Aydeniz:2021olz to approximately diagonalize the neutrino many-body Hamiltonian on a quantum computer. Further, real-time evolution using trotterization may allow for calculations of transition probabilities for interacting neutrinos. However, practical limitations of current quantum hardware prevent studies of larger systems in these earlier quantum simulations, i.e., limited number of unitary operations with low accuracy. More recently in Refs. Amitrano:2022yyn and Illa:2022zgu , trapped-ion quantum devices were utilized to perform the simulations eight and twelve neutrinos respectively, thanks to the all-to-all qubit connectivity in trapped-ion based architecture.

An extension of the frequently adopted two-flavor framework to three flavors in the mean-field approximation has revealed several unique phenomena, e.g., multiple spectral splits, in collective neutrino oscillations Fogli:2008fj ; Duan:2008prl ; Dasgupta:2008prd ; Dasgupta:2009prl ; Dasgupta:2010prd ; Friedland:2010prl ; Airen:2018nvp ; Chakraborty:2019wxe ; Shalgar:2021wlj . To see whether these differences translate to the many-body picture, the neutrino many-body problem has recently been analyzed in a three-flavor setting Siwach:prd2023 . The Hamiltonian given in Eq. (18) can be generalized to the three-flavor case:

where the generators $Q_{i^{\prime}}$ are given by

with $\lambda$’s as the Gell-Mann matrices. The auxiliary vector $\vec{B}$ is given by

Here, the oscillation frequencies are

where $\delta m^{2}=m_{2}^{2}-m_{1}^{2}$, and $\Delta m^{2}\approx|m_{3}^{2}-m_{2}^{2}|\approx|m_{3}^{2}-m_{1}^{2}|$, and $E$ is the energy of neutrino.

Due to computational limitations, only a system of up to five neutrinos could be considered in this treatment. To quantify the entanglement, the von-Neumann entanglement entropies $S(\omega_{q})$ and the components of the neutrino polarization vectors $\vec{P}_{q}$ are calculated (straightforward generalizations of the definitions in Sec. 4.2) . In the three-flavor case, two components of the total polarization vector (of the entire ensemble), namely $P_{3}$ and $P_{8}$, are conserved through the course of flavor evolution, whereas, in the two-flavor case only the total $P_{z}$ is conserved. It was found that the entanglement in the three-flavor case can be significantly larger in comparison to the two-flavor case for an initial condition with all neutrinos in the electron flavor.

Here we further investigate the entanglement in three-flavor neutrino many-body system with an initial state different from those considered in Ref. Siwach:prd2023 . Shown in Fig. 3 are the results for the time-evolution of a five-neutrino system with an initial state $\ket{\psi}=\ket{\nu_{e}\nu_{\mu}\nu_{\mu}\nu_{\tau}\nu_{\tau}}$ in both normal (NO) and inverted (IO) mass orderings. The vacuum oscillation frequencies of the neutrinos are chosen to be $\omega_{q}=q\,\omega_{p}$ and $\Omega_{q}=q\,\Omega_{p}$, with $\omega_{p}$ and $\Omega_{p}$ defined by Eq. (31), with a neutrino energy $E=10\,$MeV. $\kappa=10^{-17}\,$MeV is a suitably chosen scaling factor in terms of which all the other dimensionful quantities are defined in the numerical calculations. The ordering of the asymptotic $P_{3}$ values with respect to the frequency modes $q$ is observed to be the same in both NO and IO. However, the magnitudes are significantly different. The $P_{8}$ values on the other hand show differences in both magnitude and ordering. One can see from the entanglement entropies $S(\omega_{q})$ that the neutrino in frequency mode $q=4$ (3) is maximally entangled in NO (IO), whereas this neutrino is the least entangled in IO (NO), respectively. Further information on the entanglement can be obtained from the graphical representation of projection of polarization vectors $P_{3}$ and $P_{8}$ values in the $\hat{e}_{3}$-$\hat{e}_{8}$ plane, as given in Ref. Siwach:prd2023 .

From the probabilities of finding a neutrino in a particular mass eigenstate, $P_{\nu_{i}}$, as shown in the right panels of Fig. 3, the mixing of different mass eigenstates can be investigated. In NO, the neutrino which is maximally entangled, i.e. $q=4$, has $\sim 40\%$ contribution from first and second mass eigenstates each and $\sim 20\%$ contribution from third mass eigenstate. The least entangled $q=4$ neutrino has $\sim 80\%$ contribution from third mass eigenstate and $\sim 10\%$ contribution from first and second mass eigenstates each. In IO, the maximally entangled neutrino in $q=3$ frequency mode has $\sim 50\%$ contribution from the first mass eigenstate and $\sim 25\%$ contribution from second and third mass eigenstates each. Therefore, the results in different mass orderings differ significantly.

We notice that the entanglement in many-body neutrino systems depends substantially on the initial state. These differences enhance in three-flavor case as compared to the two-flavor case. For example, the entanglement entropy for $N=5$ neutrino system with an initial state $\ket{\nu_{e}\nu_{e}\nu_{\mu}\nu_{\mu}\nu_{\tau}}$ (see Ref. Siwach:prd2023 ) in the case of NO is maximum for the neutrino in the highest frequency mode. Whereas in the case of an initial state $\ket{\nu_{e}\nu_{\mu}\nu_{\mu}\nu_{\tau}\nu_{\tau}}$, the entropy is maximum for the second highest frequency mode $q=4$ neutrino. Several other differences can be observed in other properties and in inverted mass ordering. Therefore, to better understand the neutrino many-body system, similar calculations should be performed for a significantly larger system considering various initial states.

Many-body quantum dynamics of dense neutrino systems is an active area of research that in recent years has shown a rapid development in terms of understanding and analyzing many body correlations, instabilities, and dynamical phase transitions, with various groups attempting to investigate the problem using different types of classical and quantum computational tools. Methods adapted from quantum information have shown great potential for further studying many body effects. The ongoing quest to augment our understanding of these quantum dynamics, and collective flavor oscillations in general, stems from the importance of these phenomena in influencing the neutrino transport and nucleosynthesis in astrophysical environments like core-collapse supernovae and binary neutron star mergers, which are important open problems in the domain of high-energy/nuclear astrophysics.