# Trapped-Ion Quantum Simulation of Collective Neutrino Oscillations

Valentina Amitrano  Dipartimento di Fisica, University of Trento, via
Sommarive 14, I–38123, Povo, Trento, Italy  INFN-TIFPA Trento Institute of
Fundamental Physics and Applications, Trento, Italy  Alessandro Roggero
Dipartimento di Fisica, University of Trento, via Sommarive 14, I–38123, Povo,
Trento, Italy  INFN-TIFPA Trento Institute of Fundamental Physics and
Applications, Trento, Italy  Piero Luchi  Dipartimento di Fisica, University
of Trento, via Sommarive 14, I–38123, Povo, Trento, Italy  INFN-TIFPA Trento
Institute of Fundamental Physics and Applications, Trento, Italy  Francesco
Turro  Dipartimento di Fisica, University of Trento, via Sommarive 14,
I–38123, Povo, Trento, Italy  INFN-TIFPA Trento Institute of Fundamental
Physics and Applications, Trento, Italy  Luca Vespucci  Dipartimento di
Fisica, University of Trento, via Sommarive 14, I–38123, Povo, Trento, Italy
ECT*, European Center for Theoretical Studies in Nuclear Physics and Related
Areas, Strada delle Tabarelle 286, Trento, Italy  Francesco Pederiva
Dipartimento di Fisica, University of Trento, via Sommarive 14, I–38123, Povo,
Trento, Italy  INFN-TIFPA Trento Institute of Fundamental Physics and
Applications, Trento, Italy

(August 28, 2024)

###### Abstract

It is well known that the neutrino flavor in extreme astrophysical
environments changes under the effect of three contributions: the vacuum
oscillation, the interaction with the surrounding matter, and the collective
oscillations due to interactions between different neutrinos. The latter adds
a non-linear contribution to the equations of motion, making the description
of their dynamics complex. In this work we study various strategies to
simulate the coherent collective oscillations of a system of $N$ neutrinos in
the two-flavor approximation using quantum computation. This was achieved by
using a pair-neutrino decomposition designed to account for the fact that the
flavor Hamiltonian, in the presence of the neutrino-neutrino term, presents an
all-to-all interaction that makes the implementation of the evolution
dependent on the qubit topology.

We analyze the Trotter error caused by the decomposition demonstrating that
the complexity of the implementation of time evolution scales polynomially
with the number of neutrinos and that the noisy from near-term quantum device
simulation can be reduced by optimizing the quantum circuit decomposition and
exploiting a full-qubit connectivity. We find that the gate complexity using
second order Trotter-Suzuki formulae scales better with system size than with
other decomposition methods such as Quantum Signal Processing. We finally
present the application and the results of our algorithm on a real quantum
device based on trapped-ions qubits.

##  I Introduction

Quantum computation can provide an enormous advantage for the physical
description of many-body quantum systems due to the fact that it does not
necessarily require an exponential scaling of computational resources as the
size of the system increases [1, 2]. The prospect of employing quantum devices
to study standard model physics has led to a world wide effort to design
algorithms and apply them to currently available quantum platforms (see e.g.
[3] for a recent review).

Given a quantum system in some fixed initial state $\ket{\Psi_{0}}$, its time
evolution under the Hamiltonian $H$ is given by the action of the real-time
evolution operator:

$U(t)=e^{-iHt}\,,$ (1)

transforming the state according to the time-dependent Schrödinger equation
$\ket{\Psi(t)}=U(t)\ket{\Psi_{0}}$. In general, a direct approach based on
this description faces an exponentially growing cost on classical computers as
the size of the system is increased due to both the enormous memory
requirements to encode the states of the system and to operational cost needed
to perform matrix multiplications. Important exceptions to this behavior are
found, for instance, in stabilizer states [4, 5] or in systems with low levels
of bipartite entanglement [6, 7].

According to the Deutsch model [8] of a quantum computer, given a system in a
pure state $\ket{\varphi_{0}}$, a quantum algorithm consists of a unitary
transformation $U$ which produces a certain final state
$\ket{\varphi_{f}}=U\ket{\varphi_{0}}$ according to quantum mechanical rules.
It is then possible to perform quantum measurements yielding the probability
of finding the systems in a given state of a given basis, probability that
constitutes the result of the calculation. The Solovay–Kitaev theorem [9]
demonstrates that there is a finite set of quantum gates, which can
approximate with arbitrary accuracy any unitary transformation $U$. In this
sense, the Deutsch model is universal.

An interesting many-body system amenable for exploring simulations of the time
evolution on a quantum computer is that of collective flavor oscillations of
neutrinos caused by forward neutrino-neutrino scattering. These are predicted
to occur in extreme astrophysical environments like core-collapse supernovae,
neutron star mergers and the early universe [10, 11, 12, 13, 14, 15, 16]. The
description of flavor oscillations is a crucial aspect of such studies, since
the physics of matter under extreme conditions is strongly flavor-dependent
[17, 18] and moreover the energy spectrum is different for different neutrino
flavors [19]. A star with a sufficiently large mass ($\gtrapprox
8\,\text{M}_{\odot}$) undergoes a gravitational collapse that can result in a
neutron star or a black hole. During the collapse it emits a very large amount
of energy ($\sim 10^{53}$ erg) in the form of a large number of neutrinos
($\sim 10^{58}$). In a few seconds the $10\%$ of the gravitational mass of the
star is converted into neutrinos flowing with an energy $E_{\nu}\simeq(10\div
30)$ MeV. The evolution of this neutrino sea plays a fundamental role in
supernova collapse phenomena. They are in fact the main carriers of the lepton
number within the reactions taking place in the interior. Moreover, they are
responsible for the loss of entropy, and they can increase the instability of
the star by eventually generating the explosion [20]. A simple diagram showing
the regions where different neutrino processes are active in a core collapse
supernovae, assumed to have spherical symmetry, is shown in Fig. 1. Collective
neutrino oscillations are generally expected to be dominant in a range of
intermediate distance from the core ($\sim 100$ km) where the density of
neutrinos is large while the external lepton electron density is not
sufficiently large to suppress flavor oscillations [20]. The external shell,
in which the neutrino density is lower, is instead dominated by vacuum
oscillations and interactions with the surrounding matter leading to the MSW
effect [21].

A full description of the dynamical evolution of flavor in these processes is
hindered by the large computational cost required to carry out simulations
with large numbers of interacting neutrinos. A common approach adopted to
circumvent the problem is to use a mean-field approximation for the equation
of motions allowing to study large-scale systems with complex geometries [22,
23, 24, 25]. A full-treatment of correlation effects in the complete many-body
evolution can however be attained in relatively small systems with $O(10)$
neutrinos [26, 27, 28], in situations of large symmetry [29, 30, 31, 32]
and/or small levels of bipartite entanglement using tensor network methods
[33, 34, 35]. Another class of approaches that retains some of the
correlations while maintaining numerical efficiency are semi-classical
methods, recently applied to the neutrino problem in Ref. [36].

Quantum simulations offer an alternative to explore out-of-equilibrium flavor
dynamics in regimes that are not accessible by these classical approaches.
Early calculations on small systems with up to $4$ neutrinos have been carried
out on both digital quantum computers [37, 38] and quantum annealers [39]
showing that both careful algorithm design and error mitigation techniques
have to be considered when tackling this challenging problem on current
generation devices. In this work we propose an efficient quantum algorithm to
describe the evolution of the flavor state of a many-neutrino system using a
digital quantum computer and paying attention to the optimal quantum gate
decomposition and to the complexity of the quantum circuit needed for the
simulation. We analyze in detail the scaling of the Trotter error, the number
of operation needed to perform the evolution at a fixed error and the
complexity of the quantum gate decomposition of the evolution operator paying
attention to the machine-aware compilation that has to consider the topology
of the qubit system.

Figure 1: Sketch of the environment in a core-collapse supernova assuming
spherical symmetry. Near the proto-neutron star the large local neutrino
density causes an important effect from pair-wise neutrinos scattering while
at large radii the interaction with electrons prevails instead. Moreover we
assume a narrow cone of forward peaked emitted neutrinos.

In Sec. II we present the description of the physical system of several
neutrinos used in this work. In Sec. III we analyze the decomposition of the
unitary propagator in Eq. (1) containing an all-to-all interaction. In order
to do that, we exploit the pair property of the potential showing some of the
advantages of a full-qubit connectivity which guarantees greater freedom in
the decomposition and therefore less complexity of the final implementation.
In Sec. IV we present the optimal decomposition technique we used to found the
best quantum circuit for the trapped-ion quantum device used to perform the
real quantum simulation (a Quantinuum System Model (QSM) H1-2 trapped ion
device). In Sec. V we report the results obtained from the real quantum
simulations for a single Trotter step propagation using different time steps
and for multi-step long time evolution in the cases of $N=4$ and $N=8$
neutrinos. Finally, in Sec. VI we present a detailed derivation of the Trotter
error scaling introduced by the decomposition of $U(t)$ and the time
discretization. We also compare the gate complexity to that of other
decomposition methods such as qubitization. Finally, we demonstrate the
usefulness of the decomposition proposed in Sec. III which ensures that the
number of necessary operations follows a low degree polynomial with the number
of particles.

##  II SU(2) model of a many-neutrino system

The Hamiltonian of a many-neutrino system including the effect of forward
scattering is equivalent to an all-to-all coupled spin system and thus
represents an interesting many-body quantum problem governed by the weak
interaction. The first approximation we make is to consider only two neutrino
flavors in the description: the electron flavor $\nu_{e}$ and a single heavy
flavor $\nu_{x}$ which is a combination of $\mu$ and $\tau$ neutrinos. This
can be justified if the mixing angle $\theta_{13}=0$ as shown in [40]. In this
way the general flavor state of each neutrino, given by the two-flavor
superposition $\ket{\Psi}=\alpha\ket{\nu_{e}}+\beta\ket{\nu_{x}}$, can be
fully encoded in the state of a qubit by means of the mapping:

$\ket{\nu_{e}}\longmapsto\ket{0},\quad\ket{\nu_{x}}\longmapsto\ket{1}\,.$ (2)

The flavor Hamiltonian in this basis can be decomposed into three main terms
[20, 41]: (1) a 1-body contribution describing vacuum mixing $H_{vac}\coloneqq
H^{(1)}$ which takes into account the flavor oscillations of each neutrino due
to the misalignment between flavor states and mass eigenstates; (2) a second
1-body part describing the coupling to external matter leading to the MSW
effect and (3) the neutrino-neutrino interaction term $H_{\nu\nu}:=H^{(2)}$
generated by forward scattering. In this work we will neglect the second
contribution, since we are interested in describing the dynamics in the
coherent oscillation-dominated region, and only focus on a simplified
Hamiltonian of the form

$\begin{split}H=H^{(1)}+H^{(2)}&=\sum_{i=0}^{N-1}h_{i}+\sum_{i<j}^{N-1}h_{ij}\\\
&=\sum_{i=0}^{N-1}\boldsymbol{b}\cdot\boldsymbol{\sigma}_{i}+\sum_{i<j}^{N-1}J_{ij}\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}\,,\end{split}$
(3)

where we used bold symbols to denote 3-dimensional vectors. The vectors
$\boldsymbol{\sigma}_{i}=(\sigma_{i}^{(x)},\sigma_{i}^{(y)},\sigma_{i}^{(z)})$
are formed by the Pauli matrices acting on the $i$-th neutrino. Here and in
the following we suppress the identity operators acting on the other spins,
for example
$\sigma^{(x)}_{1}=\mathbb{1}\otimes\sigma^{(x)}\otimes\mathbb{1}\otimes\dots\otimes\mathbb{1}$.
In the flavor basis, the vector $\boldsymbol{b}$ in the first term describes
vacuum mixing of neutrinos with the same energy and is given explicitly by

$\boldsymbol{b}=\frac{\delta
m^{2}}{4E}(\sin(2\theta_{\nu}),0,-\cos(2\theta_{\nu}))\;.$ (4)

In this expression $\delta m^{2}=m_{2}^{2}-m_{1}^{2}$ is the square mass
difference between mass eigenstates, which is of order $10^{-4}\text{
eV}^{2}$, $\theta_{\nu}$ is the mixing angle that we took equal to
$\theta_{\nu}=0.195$, and $E$ is the neutrino energy. The two-body interaction
term is described by the $SU(2)$ invariant product of Pauli matrices in which,
also in this case, the tensor products with the other particles are implicit.
For example

$\sigma^{(x)}_{0}\sigma^{(x)}_{2}=\sigma^{(x)}\otimes\mathbb{1}\otimes\sigma^{(x)}\otimes\mathbb{1}\otimes\dots\otimes\sigma^{(x)}\otimes\mathbb{1}\;.$
(5)

The coupling constant of the pair interaction can be written explicitly as
follows

$\begin{split}J_{ij}&=\frac{\sqrt{2}G_{F}}{V}(1-\cos(\theta_{ij}))\coloneqq\frac{\mu}{N}(1-\cos(\theta_{ij}))\,,\end{split}$
(6)

and depends on the relative angle of propagation

$\cos(\theta_{ij})=\frac{\boldsymbol{p}_{i}\cdot\boldsymbol{p}_{j}}{\|\boldsymbol{p}_{i}\|\|\boldsymbol{p}_{j}\|}\;,$
(7)

where $\boldsymbol{p}_{i}$ is the momentum of $i$-th neutrino. This means that
the neutrinos that interact the most are those that propagate in directions
with a larger relative angle. In Eq. (6) we introduced the energy scale
$\mu=\sqrt{2}G_{F}n_{\nu}$ where $G_{F}$ if the Fermi constant, $n_{\nu}=N/V$
the neutrino number density and $V$ the volume of the system. We choose the
neutrino energy in order to obtain the same coupling constant for the one body
and two body energies

$\frac{\mu}{N}=\frac{\delta m^{2}}{4E}\;,$ (8)

and measure time in units of $\mu^{-1}$.

Figure 2: Exact evolution of a system of $N=4$ neutrinos with initial state
$\ket{\Psi_{0}}=\ket{0011}$. The panel (a) shows the evolution of the
expectation value $\expectationvalue{Z_{i}}$ while the panel (b) shows the
evolution of the flavor inversion probability.

As done in previous work [37], we take a simple grid of angles

$\theta_{ij}=\arccos(0.9)\frac{\absolutevalue{i-j}}{(N-1)}\,,$ (9)

meant to reproduce a narrow cone of forward peaked neutrinos in accordance to
the geometry displayed in Fig. 1. With this choice of angular distribution,
and for even $N$, the neutrino Hamiltonian in Eq. (3) turns out to be
symmetric under the particle exchange

$\begin{split}\nu_{k}&\longleftrightarrow\nu_{N-1-k}\;,\end{split}$ (10)

for $k=0,\dots,N/2$. The initial state of the system, used for all the
simulations presented in this work, consists of setting the first $N/2$
neutrinos in the $\ket{\nu_{e}}$ flavor state and the other $N/2$ in the
$\ket{\nu_{x}}$ state. In this way the initial state is symmetric under the
composition of particle exchange and flavor inversion. For $N=4$ neutrinos,
for instance, the initial state is given by:

$\ket{\Psi_{0}}=\ket{\nu_{0}\nu_{1}\nu_{2}\nu_{3}}=\ket{0011}=\ket{\nu_{e}\nu_{e}\nu_{x}\nu_{x}}\;.$
(11)

One can obtain the exact time evolution by directly performing a matrix
multiplication $\ket{\Psi(t)}=U(t)\ket{\Psi_{0}}$. In the spin basis, the
flavor content of an individual neutrino is obtained from the expectation
value of the Pauli matrix $\sigma^{(z)}_{i}=Z_{i}$

$\expectationvalue{Z_{i}(t)}=\expectationvalue{Z_{i}}{\Psi(t)}\;,$ (12)

where the tensor products are implicit for the other particles. In a similar
way, the flavor inversion probability $P_{i}(t)$ can be expressed as

$P_{i}(t)=\frac{|\expectationvalue{Z_{i}(0)}-\expectationvalue{Z_{i}(t)}|}{2}\,.$
(13)

We display the exact evolution of both quantities in Fig. 2: the top panel
shows results for $\expectationvalue{Z_{i}(t)}$ while the bottom panel
displays the inversion probabilities. As expected from the exchange symmetry
in Eq. (10) and the asymmetric choice of initial state $\ket{\Psi_{0}}$, the
flavor evolution of neutrinos $0$ and $1$ are the mirror image of neutrinos
$3$ and $2$ (respectively). This is reflected in the equivalence of inversion
probabilities for these neutrinos (bottom panel of Fig. 2). Due to the
presence of this symmetry, in the rest of this work we will show results for
inversion probabilities only.

##  III Implementation of the time evolution operator

In order to carry out a quantum simulation one always needs two ingredients:
(1) a state encoding map and (2) a way to map operators into quantum gates.
Due to the two-flavor approximation presented in Sec. II, the flavor state of
neutrinos can be directly encoded into a qubit according the map in Eq. (2).
In the case of a digital quantum simulations the operator $U(t)$ must then be
decomposed into a sequence of quantum gates from a fixed set. In this way the
initial state enconding the flavor state $\ket{\Psi_{0}}$ is evolved under a
sequence of unitary transformations overall implementing the real-time
propagator. A projective measurement of the final state of the qubits
eventually allows to extract flavor observables from the simulation. In the
computational basis, the operator $U(t)$ is represented by a $2^{N}\times
2^{N}$ unitary matrix which must be decomposed into the gate set provided by
the machine, usually composed by single or 2-qubit elementary gates. One way
to decompose this operator is to explicitly exploit the actual interaction of
the physical system which occurs in pairs and which can therefore be
implemented by considering only a pair of qubits at a time. The approach
followed in Ref. [37] uses the exact pair propagator. This can be obtained in
two steps, first one symmetrizes the 1-body term and expressing the total
Hamiltonian as a sum of 2-body terms, namely:

$H=\sum_{i<j}^{N-1}\bigg{(}\boldsymbol{b}\cdot\frac{(\boldsymbol{\sigma}_{i}+\boldsymbol{\sigma}_{j})}{N-1}+J_{ij}\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}\bigg{)}:=\sum_{i<j}^{N-1}H_{ij}\;.$
(14)

Then the total propagator can be approximated by a product of pair
propagators:

$U(t)\approx\prod_{i<j}^{N-1}e^{-iH_{ij}t}\,.$ (15)

The approximation introduces an error of order $\mathcal{O}(t^{2})$ due to the
non-commutativity of the symmetrized 2-body terms, namely
$\commutator{H_{ij}}{H_{ik}}\not=0$. A first improvement of this
implementation can be made considering that the entire 1-body and 2-body terms
commute, namely $\commutator{H^{(1)}}{H^{(2)}}=0$ and therefore their
separation does not introduce any errors:

$\begin{split}U(t)=e^{-iH^{(2)}t}e^{-iH^{(1)}t}:=U_{2}(t)U_{1}(t)\,.\end{split}$
(16)

Subsequently, the 2-body term alone can be efficiently implemented as a pair
decomposition:

$U_{2}(t)\approx\widetilde{U}_{2}(t)=\prod_{i<j}^{N-1}e^{-ih_{ij}t}:=\prod_{i<j}^{N-1}u_{ij}(t)\,.$
(17)

The implementation in Eq. (15) might lead to an increase in the error due to
the lack of commutativity between individual one and two-body contributions.
However, the approach can still be useful in the case of a particle-dependent
external field, in which $\boldsymbol{b}\to\boldsymbol{b}_{i}$ as it occurs
when different neutrinos have different energies $E\to E_{i}$ (a necessary
ingredient to observe spectral splits [20, 28]).

###  III.1 Qubit connectivity and pair ordering

The Hamiltonian in Eq. (3) contains an all-to-all interaction term. Using the
implementation of $U_{2}(t)$ from Eq. (17) we have to make all qubits interact
with all the others during the simulation at least once. This fact implies
that in a quantum computation the circuit implementing the sequence of
$U_{2}(t)$ operators for each particle pair must be adapted to the particular
topology of the specific quantum device employed for the simulation. As shown
in Ref. [37] it is possible to construct $U_{2}(t)$ with only linear
connectivity and with a gate depth of $N$ using a swap network (the same
scheme was later adopted for tensor-network simulations in Ref. [33]). The
algorithm consists in applying the $u_{ij}$ propagator to a qubit pair
followed by a $SWAP$ gate which exchanges the qubit state. The new unitary
$w_{ij}$ is thus

$w_{ij}=SWAP\times u_{ij}\;,$ (18)

where in the computational basis the $SWAP$ unitary is

$SWAP=\begin{pmatrix}1&0&0&0\\\ 0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&1\\\
\end{pmatrix}\;.$ (19)

For example, the swap network required for the case of $N=4$ is depicted in
Fig. 3.

Figure 3: Swap network for $N=4$ qubits that implements the 2-body propagator
$U_{2}(t)$ using a chain of linearly connected qubits. Vertical lines
represent qubit interaction and crosses the SWAP operation. The pair ordering
influences the total error, here we use the map
$(q_{0},q_{1},q_{2},q_{3})=(\nu_{0},\nu_{2},\nu_{1},\nu_{3})$.

It should be pointed out that the error introduced by the approximation in Eq.
(17) also depends on the order in which the pairs interact. This is due to the
dependence of the error in Eq. (17) on sums of the commutators
$\commutator{h_{ij}}{h_{kl}}$ taken with a given order (see Appendix A for
additional details). It is then worth looking for an optimal ordering allowing
to maximize the cancellations between the commutators minimizing the
decomposition error. Note that, once the network structure shown in Fig. 3 is
fixed, the ordering can be chosen by varying the initial encoding of each
neutrino into the qubits. The structure of the swap network, dictated by the
available qubit topology, imposes a constraint on the possible ordering of the
pair propagators. In principle, any ordering could be achieved by adding
additional $SWAP$ gates but at the cost of increasing the complexity in terms
of depth and number of 2-qubit gates of the scheme. For our Hamiltonian, and
in the case of $N=4$, we found that the best interaction order would be the
one described by the network in Fig. 4. Such network cannot be expressed using
a swap network with $4$ layers if restricted to the use of linear qubit
connectivity, since at least $5$ layers would be needed. With all-to-all
connectivity, however, this algorithm can be implemented using one less layer.
This is due to the fact that each layer is full, in the sense that the maximum
number of possible operations at the same time is performed.

Figure 4: Implementation of the 2-body propagator $U_{2}(t)$ on a qubit system
with all-to-all connectivity. This scheme can implement the optimal ordering
of pair to minimize the decomposition error. For our system this corresponds
to the qubit to neutrino map
$(q_{0},q_{1},q_{2},q_{3})=(\nu_{0},\nu_{1},\nu_{2},\nu_{3})$.

Determining the optimal ordering for large systems is not feasible in general
as this would require a super-exponential cost in the system size $N$. For
large systems a randomization procedure for the order could prove valuable to
control the error [42].

Figure 5: The top panels show the inversion probability after single Trotter
step for different values of the time steps $dt$. Panel (a) shows the
evolution for neutrinos $\nu_{0}$ and $\nu_{3}$ and panel (b) for $\nu_{1}$
and $\nu_{2}$. The dotted curve is the exact evolution, the solid blue line
the one obtained by applying the implementation proposed in this work (Eq.
(16) and Eq. (17)) together with the scheme in Fig. 4. The dashed orange and
red lines are the evolution using the implementation proposed in Ref. [37] and
described by Eq. (15) and the scheme from Fig. 3. Panel (c) shows the error in
spectral norm between the exact propagator and the two Trotter decompositions.

We analyze the effect of different orderings on a single time-step in Fig. 5.
The two top panels show the evolution of the inversion probability as a
function of the time-step $dt\in[0,40]\,\mu^{-1}$ for the same initial state
as in Fig. 2 and for different implementation of the propagator: the dotted
curve is the exact evolution, the solid blue line the one obtained by applying
the implementation proposed in this work (Eq. (16) and Eq. (17)) together with
the optimal ordering (OO) from Fig. 4. The dashed orange and red lines are the
evolution using the implementation proposed in Ref. [37] and described by Eq.
(15) and the swap network (SN) scheme from Fig. 3. Panel (a) shows results for
neutrinos $\nu_{0}$ and $\nu_{3}$ while panel (b) for $\nu_{1}$ and $\nu_{2}$.
The results highlight the preservation of the exchange symmetry from Eq. (10)
of the scheme introduced in this work and afforded by the separation between
the one and two-body contributions in Eq. (16). Note that this property,
observed for $N=4$, does not hold for general system sizes. For instance, we
were not able to find an ordering choice preserving this property while
keeping the minimum number of layers for the case of $N=8$.

Panel $(c)$ of Fig. 5 shows the error, for a single Trotter step, of the two
approximations for the time evolution operator as a function of the time-step
$dt$. We calculate the error using the infinity (spectral) norm

$\varepsilon(dt)=\norm{\widetilde{U}_{2}(dt)-U_{2}(dt)}_{\infty}\,.$ (20)

The results displayed in the top panels show that the error in the inversion
probability is lower for the (SN) approximation when $dt<30\,\mu^{-1}$. This
is mostly an effect of choosing a particular initial state $\ket{\Psi_{0}}$
and this specific observable. For general initial states and observables, the
error displayed in panel $(c)$ shows that the (OO) approximation has indeed
the smallest worst case error for all time-steps.

Figure 6: The top panel shows the time evolution of the inversion probability
for neutrino $\nu_{1}$ using a time step $dt=16\mu^{-1}$. The dotted black
curve is the exact evolution, the orange dashed curve Trotter decomposition
with the best possible swap network (SN) while the solid blue curve uses the
optimal ordering (OO) achievable using the full-connectivity.

The main advantage of using the optimal order is evident when the goal of the
simulation is to describe the evolution of the system for a long total time
$T$ in which the accumulation of the error is dominant. In Fig. 6 we plotted
the time evolution of the inversion probability for the neutrino $\nu_{1}$
(set in the electron flavor $\ket{0}$ at the beginning) for a long total time
$T=1200\,\mu^{-1}$ using a time step $dt=16\,\mu^{-1}$. For each time, we
apply sequentially the approximate propagator

$\ket{\Psi(kdt)}=\widetilde{U}_{2}(dt)^{k}U_{1}(dt)^{k}\ket{\Psi_{0}}\,,$ (21)

where $\widetilde{U}_{2}(dt)$ is defined in Eq. (17). We employ two different
orderings: the best possible one that is achievable with a swap network (SN)
with linear connectivity as presented in Fig. 3 above, and denoted by the
dashed orange curve in panel (a) of Fig. 6, and the optimal ordering (OO) one,
achievable by allowing all-to-all connectivity shown in Fig. 4, presented as
the solid blue curve in Fig. 6. As we can see, the results obtained using (OO)
are much more stable than with (SN) and allow to reach long evolution times
even with large time steps.

##  IV Optimized quantum circuit

In order to implement the unitary propagator we need to decompose it as a
sequence of elementary gates from the universal gate set used by the quantum
machine. The 1-body part is trivial because it is the tensor product of the
same single-qubit gate applied to each qubit:

$U_{1}(dt)=\bigotimes_{i=0}^{N-1}\exp\left(-i\boldsymbol{b}\cdot\boldsymbol{\sigma}_{i}dt\right)\;.$
(22)

The two-body part of the propagator, as approximated in Eq. (17), is the
product of pair terms of the form:

$u_{ij}(dt)=e^{-idtJ_{ij}(X\otimes X+Y\otimes Y+Z\otimes Z)}\,.$ (23)

In the case of full qubit-connectivity we do not need to add the $SWAP$ gate
after the pair interaction (see Eq. (18)) and we can directly implement the
$u_{ij}(dt)$ operator. This results in a decomposition with a smaller number
of single-qubit gates. In fact using the result in Ref. [43], the optimal
CNOT-based decomposition for the $SU(2)$ invariant unitary operator $u_{ij}$
in Eq. (23) can be written as

$\Qcircuit@C=0.5em@R=1em{&\qw\targ{}\gate{R_{z}(\phi)}\ctrl{1}\qw\targ{}\gate{R_{z}(\frac{\pi}{2})}\qw\\\
\gate{R_{z}(-\frac{\pi}{2})}\ctrl{-1}\gate{R_{y}(-\phi)}\targ{}\gate{R_{y}(\phi)}\ctrl{-1}\qw\qw}$
(24)

where we introduced the angle $\phi=-(2dtJ_{ij}+\pi/2)$.

$\Qcircuit@C=0.5em@R=1em{&\gate{U_{q}(-\frac{\pi}{2},\frac{\pi}{2})}\multigate{1}{ZZ}\gate{U_{q}(\frac{\pi}{2},\pi)}\gate{R_{z}(2\alpha-\frac{3}{2}\pi)}\multigate{1}{ZZ}\gate{U_{q}(-\frac{\pi}{2},\frac{\pi}{2})}\qw\multigate{1}{ZZ}\gate{U_{q}(\frac{\pi}{2},\pi)}\qw\\\
\qw\ghost{ZZ}\gate{U_{q}(2\alpha,\frac{\pi}{2})}\gate{R_{z}(-\frac{3}{2}\pi)}\ghost{ZZ}\gate{U_{q}(2\alpha,\frac{\pi}{2})}\gate{R_{z}(-\frac{\pi}{2})}\ghost{ZZ}\qw\qw}$
Figure 7: Decomposition of the two-body propagator $u_{ij}(dt)$ in Eq. (23)
using native gates and setting $\alpha=-dtJ_{ij}$.

When performing a quantum simulation, it is always advisable to write the
circuit using the physical gates actually implemented by the machine and
optimize it. This procedure allows to reduce the number of gates necessary for
the simulation. In this work we simulated the real-time evolution of a $N=4$
and $N=8$ neutrino systems using the Quantinuum System Model (QSM) H1-2
trapped ion device, similar to the machine described in Ref. [44]. As other
trapped ion-based quantum devices it offers the advantage of a full
connectivity between the qubits and high fidelity: a single-qubit error of
$\sim 10^{-4}$ and a two-qubit gate error of $\sim 10^{-3}$. The native
single-qubit gates used are:

$R_{z}(\lambda)=\begin{pmatrix}e^{-i\lambda/2}&0\\\
0&e^{i\lambda/2}\end{pmatrix}$ (25)

and:

$U_{q}(\theta,\varphi)=\begin{pmatrix}\cos{\theta/2}&-ie^{-i\varphi}\sin{\theta/2}\\\
-ie^{i\varphi}\sin{\theta/2}&\cos{\theta/2}\end{pmatrix}\,.$ (26)

Currently only two possible values $\theta=\pi,\pi/2$ are natively available
on the device. This means that up to $4$ physical gates are required to
compile a general $SU(2)$ unitary. The 2-qubit gate is:

$ZZ=\exp\left(-i\frac{\pi}{4}Z\otimes Z\right)=\begin{pmatrix}1&0&0&0\\\
0&i&0&0\\\ 0&0&i&0\\\ 0&0&0&1\end{pmatrix}\,.$ (27)

We can replace, in the circuit from Eq. (24), the following $ZZ$-based
implementation of the $CNOT$ gate:

$\Qcircuit@C=0.5em@R=0.7em{&\ctrl{1}\qw\raisebox{-25.00003pt}{=}\qw\multigate{1}{ZZ}\qw\gate{R_{z}(-\frac{\pi}{2})}\qw\\\
\targ\qw\gate{U_{q}(-\frac{\pi}{2},\frac{\pi}{2})}\ghost{ZZ}\gate{U_{q}(\frac{\pi}{2},\pi)}\gate{R_{z}(-\frac{\pi}{2})}\qw.}$
(28)

Then, using commutation properties, all possible simplifications including
cancellation and gate mergers, and finally translating all single-qubit gates
into native gates we obtain the circuit displayed in Fig. 7 that implements
the two-body part evolution for one pair of neutrinos using only native gates.
The angle parameter $\alpha$ is defined, for every qubit pair, as
$\alpha=-dtJ_{ij}$.

A single Trotter step is implemented by first applying the one-body part Eq.
(22), made up of single-qubit gates, followed by the optimal ordered sequence
of pairs interactions. A scheme for the circuit of a complete single Trotter
step is presented in Fig. 8 where each 2-qubit gate is implemented by the
circuit in Fig. 7. The evolution for a total time $T$ using a time step $dt$
is implemented by sequentially applying the scheme in Fig. 8 for $r=T/dt$
times and moving to the beginning the full 1-body part. In these cases the
overall gate count for the circuit can be further optimized using an inverted
interaction order in alternate steps (that means performing the last
interactions $(1,2)+(0,3)$ at the beginning, which is also an equally optimal
ordering) resulting in a reduction of $3N/2$ $ZZ$ gates for each step. Using
this alternating scheme can also reduce the overall approximation error since
it becomes then equivalent to a second order Trotter step with time-step $2dt$
(see Sec. VI for additional details).

Figure 8: Scheme for implementing a single Trotter step. We first apply the
single-qubit gates corresponding to the one-body propagator $U_{1}(dt)$ and
then the pair propagator using the optimal ordering implementing $U_{2}(dt)$.

##  V Results of Quantum Simulations

In this section we report the results obtained from the simulations carried
out on the Quantinuum trapped-ion device. We perform two distinct simulation:
Frst, we approximate the full time evolution with a single Trotter step with
different values of $dt$, as done in Ref. [37]; Second, we perform a multi-
step simulation fixing $dt=4\,\mu^{-1}$. The single step results are obtained
using both a system with $N=4$ neutrinos, equivalent to the one studied in
Ref. [37] on a superconducting device and in Ref. [39] on a quantum annealer,
as well as for a larger system with $N=8$. Similarly to this work, the
simulation on superconducting circuits required a number of qubits equal to
$N$ while, using the strategy described in Ref. [39], in order to map the
problem into a quantum annealer we would need at least $2^{N+2}$ physical
qubits in the ideal case of all-to-all connectivity (corresponding to $64$
qubits for $N=4$ and $1024$ qubits for $N=8$). However, due to the limited
connectivity in these devices, a large number of auxiliary qubit needs to be
employed. With the Pegasus topology implemented in the D-Wave Advantage
system, the $N=4$ simulation presented in Ref. [39] required $\approx 2000$
physical qubits and, by extrapolating the resource requirements estimated in
Fig. 9 there, a simulation with $N=8$ neutrinos will require $\approx 10^{5}$
physical qubits and is likely out of reach in the near future.

For each simulation, the circuit is repeated $M=200$ times in order to collect
statistics for the measurement outcomes. The results of the simulation are
then analyzed by calculating statistical confidence intervals using the
Bayesian approach already employed in Ref. [37]. In this section the size of
the squares representing the results in the figures cover the $68\%$
confidence interval, while errorbars give the $90\%$ interval.

###  V.1 Single Trotter step propagation

The initial flavor states for $N=4$ and $N=8$ are chosen to contain a mixture
of both $e$ and $x$ flavors as

$\rvert\Psi_{0}^{(4)}\rangle=\ket{0011}\;,\quad\quad\rvert\Psi_{0}^{(8)}\rangle=\ket{00001111}\;,$
(29)

respectively, and are prepared using single qubit gates $U_{q}$. We then apply
the Trotter step according to the diagram displayed in Fig. 8 (and it’s
generalization to the $N=8$ case) using the gate decomposition in Fig. 7. This
requires a number of two-qubit $ZZ$ gates equal to $18$ (for $N=4$) and $84$
(for $N=8$). Since the trapped-ion device has more than $4$ qubits at our
disposal, in order to minimize the execution cost of the $N=4$ experiments we
always use two sets of $4$ qubits at the same time, each one performing the
circuit sequence corresponding to two different values of $dt$. Thus, the
results presented below are all obtained from experiments using $8$ qubits
with a two-qubit gate depth given by either $9$ or $21$ (for the $N=4$ and
$N=8$ simulations respectively).

Figure 9: Single Trotter step evolution for the inversion probability starting
with the initial state $\rvert\Psi_{0}^{(4)}\rangle=\ket{0011}$. Panel (a) for
neutrinos $\nu_{0}$ and $\nu_{3}$ and (b) for $\nu_{1}$ and $\nu_{2}$. The
dotted line represent the ideal results using the exact propagator while the
solid black line indicates the ideal result obtained using one Trotter step.
The results obtained from experiments on QSM H1-2 are represented by square
data points (corresponding to a $68\%$ confidence interval) while the error
bar indicates a $90\%$ confident interval. Figure 10: Single Trotter step
evolution for the inversion probability starting with the initial state
$\rvert\Psi_{0}^{(8)}\rangle=\ket{00001111}$. Panel (a) shows result for the
pair $(\nu_{0},\nu_{7})$, panel (b) for $(\nu_{1},\nu_{6})$, panel (c) for
$(\nu_{2},\nu_{5})$ and panel (d) for $(\nu_{3},\nu_{4})$. Curves and data
points follow the same convention used in Fig. 9.

Fig. 9 shows the results for the inversion probability obtained for the
simulation of a system of $N=4$ neutrinos while the Fig. 10 for a system of
$N=8$ neutrinos. In both figures the pairs of neutrinos related by the
exchange symmetry from Eq. (10) are represented in the same panel: neutrinos
$\nu_{0}$ and $\nu_{7}$ in panel (a), $\nu_{1}$ and $\nu_{6}$ in panel (b),
$\nu_{2}$ and $\nu_{5}$ in panel (c) and $\nu_{3}$ and $\nu_{4}$ in panel (d).
In the limit of negligible machine errors, the implementation of the
propagator proposed in this work guarantees perfect symmetry under particle
exchange in the case of $N=4$ and the results obtained on the real device
respect almost always this symmetry within a confidence interval of $68\%$ and
always in the interval of $90\%$. Due to our choice of using $8$ qubits to
perform two simulations for $N=4$ in parallel, these errors might in principle
be affected by the cross-talk between the two simulations. However, by
analyzing the results of $3$ parallel executions over $12$ qubits using time
steps $dt=16,\,24,\,32\,\mu^{-1}$ we found them to be compatible with the
results shown in Fig. 9 indicating that cross-talk effects are minimal.

For $N=8$ neutrinos our implementation of the propagator respects the symmetry
of particle exchange up to $dt\approx 24\,\mu^{-1}$ as can be seen by the
theoretical results shown as solid lines in Fig. 10. The real data respects
the same symmetry within the $90\%$ confidence interval. The results obtained
for the evolution of a single step are very promising and much more compatible
with the theoretical ones than those obtained in the work [37].

In order to compare the results obtained in this work using a trapped-ion
device and with optimal ordering (OO) to the previous results obtained using a
swap network (SN) on the IBMQ Vigo supercondicting device[37], we present in
Tab. 1 the values for $\chi^{2}$ of the inversion probabilities of each
neutrino for the $N=4$ simulation. The quality of the results is assessed by
measuring the distance from the theoretical prediction of the inversion
probability of neutrino $i$ at time $t$, $P^{(th)}_{i}(t)$, of the computed
results $P_{i}(t)$ by means of the function:

$\chi_{i}^{2}=\frac{1}{10}\sum_{k=1}^{10}\frac{\left(P_{i}(kdt)-P^{(th)}_{i}(kdt)\right)^{2}}{\delta
P_{i}(kdt)^{2}}\;,$ (30)

where $\delta P_{i}(t)$ is the estimated variance (taken to coincide with the
$68\%$ confidence interval). The last row shows results from the
superconducting IBMQ device after error mitigation. In our current simulation
we did not attempt to mitigate errors and the report values are the bare
results. It can be noticed that there is an important increase of the fidelity
in the results presented in this work. This is a combined effect of the higher
gate fidelity provided by the trapped-ion device and the reduction in complex
one qubit rotations afforded by the different Trotter decomposition from Eq.
(16) adopted here. The number of general $SU(2)$ single qubit operation in a
single step for $N=4$ is in fact reduced from $40$ with the decomposition
adopted in Ref. [37] to $36$ in the present work. More importantly, the $40$
rotations adopted in the previous work have arbitrary angles due to the
combination of the pair propagator and the $SWAP$ gate and always require 3
elementary rotations to be implemented. Exploiting instead the full
connectivity and the optimal decomposition in Fig. 7, about 2/3 of the
unitaries consist of rotations of angles which are multiples of $\pi/2$. This
helps reduce the effect of coherent errors in the final results. Finally, the
results of the simulations carried out on IBMQ were obtained using a much
larger statistical sample (8192 circuit repetitions instead of 200), in order
to compare more directly the new results with the ones obtained there, we also
report in the last line of Tab. 1 the estimated $\chi^{2}$ we would have
expected to see if we reduced the statistics of the bare IBMQ results. The
same procedure cannot be consistently performed for the mitigated results,
since such estimate is also affected by systematic errors (for more details
see Refs.[37, 45, 46]). However, the strong effect of an increased gate
fidelity is still evident.

| $\nu_{0}$ | $\nu_{1}$ | $\nu_{2}$ | $\nu_{3}$
---|---|---|---|---
OO + QSM H1-2 (bare) | 1.64 | 1.59 | 1.22 | 0.64
SN + IBMQ Vigo (bare) | 424.31 | 527.64 | 545.28 | 502.60
SN + IBMQ Vigo (mit) | 71.35 | 73.64 | 126.38 | 142.72
SN + IBMQ Vigo (bare*) | 10.36 | 12.88 | 13.31 | 12.27
Table 1: Values of $\chi^{2}$ for each neutrino calculated on the results
obtained from the propagation of a single Trotter step. The results denoted
IBMQ Vigo are taken from Ref. [37] while the QSM H1-2 ones are from the
present work.

###  V.2 Multiple Trotter steps for $N=4$

In order to reach long simulation times while keeping the error under control,
the standard approach is to divide the full interval into time steps which are
then approximated using a short-time approximation. For a system of $N=4$
neutrinos, initially in the state $\ket{\Psi_{0}}$, this can be done through
the sequential application of the scheme depicted in Fig. 8 for a number of
time steps $k$ obtaining the final state

$\ket{\Psi(kdt)}=\widetilde{U}_{2}(dt)^{k}U_{1}(dt)^{k}\ket{\Psi_{0}}\,.$ (31)

For the results shown in this section we used $k\in[1,10]$ with a time step
$dt=4\,\mu^{-1}$. Also in this case, in order to exploit the specific features
of the trapped-ions machine, we have carried out the simulations in pairs in
which the first 4 neutrinos are evolved up to a time $T=kdt$ and the last 4 up
to $T=(k+1)dt$. This means that the circuit to simulate the first two times
$T=1$ and $T=2$ contains 18 $ZZ$ gates applied to the first 4 qubits and 36 to
the last 4. We summarize the gate counts needed to simulate the system up to a
certain evolution time $T=kdt$, divided into number of two-qubit $ZZ$ gates
and one-qubit $SU(2)$ unitaries, in Tab. 2. For each step we need $6\times 3$
two qubit gates while, exploiting cancellations between neighboring pair
propagators, the number of general single qubit rotations needed scales as
$32\times k+4$.

$\\#$ of steps $k$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$\\#$ of $ZZ$ gates | 18 | 36 | 54 | 72 | 90 | 108 | 126 | 144 | 162 | 180
$\\#$ of $SU(2)$ gates | 36 | 68 | 100 | 132 | 164 | 196 | 228 | 260 | 292 | 324
Table 2: Number of $SU(2)$ gates and $ZZ$ gates to evolve the system to a
fixed time $T=kdt$ to produce the results in 11.  Figure 11: Real time
evolution of the system with $N=4$ neutrinos for time step $dt=4\,\mu^{-1}$
and for a total time of $T=kdt$ with $k\in[1,10]$ using the QSM H1-2 trapped-
ion device. The top panel shows neutrinos $\nu_{0}$ and $\nu_{3}$ and the
bottom one $\nu_{1}$ and $\nu_{2}$. Dashed blue lines are the ideal evolution
that, using this small time steps, is almost the same of the Trotter
approximated propagation.

In this quantum machine measurements are carried out by default at the end of
the circuit. Our simulations include two independent circuits applied on two
separate blocks of qubits where the first ends before the second. This means
that the first 4 qubits are subject to idle errors which, in the case of
trapped-ions, are mainly related to dephasing. One way to improve the
simulation could be to use a Dynamical Decoupling approach [47] consisting in
applying to idle qubits a sequence of single-qubit gates that have the overall
effect of an identity.

##  VI Error scaling and gate cost

In this section we analyze the scaling of the Trotter error and the gate cost
of the circuit necessary to evolve a system of neutrinos as a function of the
number $N$ of neutrinos. To compare the different decomposition techniques, we
look at the scaling corresponding to a description of a system in which we fix
the neutrino density $n_{\nu}=N/V$ and we increase the dimension of the
system. Additional details on the results presented here can be found in the
Appendices.

The implementation in Eq. (17) approximates the two-body part of the
propagator $U_{2}(dt)=e^{-iH^{(2)}dt}$ using a standard first-order Trotter
formula:

$\mathcal{L}_{1}(dt)=\prod_{i<j}^{N-1}e^{-ih_{ij}dt}\;.$ (32)

By leveraging the results in Ref. [48] (see Appendix A for additional
details), we can bound the spectral norm error in this approximation as

$\varepsilon_{1}(dt):=\norm{\mathcal{L}_{1}(dt)-e^{-iHdt}}_{\infty}=\mathcal{O}(dt^{2}\mu^{2}N)\;.$
(33)

Using the union bound to obtain $\varepsilon_{1}(T)\leq r\varepsilon_{1}(dt)$
with $dt=T/r$, one finds that also the total number of steps $r_{1}$ required
to evolve for a final time $T$ while keeping the total error below $\epsilon$
scales linearly with system size. We can then bound the gate cost
$\mathcal{C}_{1}$ of a quantum circuit implementing the full evolution in
terms of the number of general $SU(4)$ 2-qubit gates needed to implement all
the steps (which, on the trapped ion system used for this work, could be
decomposed as shown in Fig. 7). For $N$ neutrinos we have $N(N-1)/2$
interaction terms in each Trotter steps so that the gate cost using a first
order product formula scales with

$\begin{split}\mathcal{C}_{1}&=\frac{N(N-1)}{2}r_{1}=\mathcal{O}\bigg{(}\frac{T^{2}\mu^{2}N^{3}}{\epsilon}\bigg{)}\,.\end{split}$
(34)

Since $SU(4)$ transformations can be decomposed with at most $3$ entangling
gates [43], as we did for the native gate set available on the QSM H1-2 system
in Fig. 7, the total count for two-qubit gates is given by $3\mathcal{C}_{1}$.

Decomposition type | Single step error | Number of steps | Circuit complexity
---|---|---|---
First order Trotter | $\mathcal{O}(dt^{2}\mu^{2}N)$ | $\mathcal{O}\left(\frac{T^{2}\mu^{2}}{\epsilon}N\right)$ | $\mathcal{O}\left(\frac{T^{2}\mu^{2}}{\epsilon}N^{3}\right)$
Second order Trotter | $\mathcal{O}(dt^{3}\mu^{3}N)$ | $\mathcal{O}\left(\frac{T^{3/2}\mu^{3/2}}{\sqrt{\epsilon}}\sqrt{N}\right)$ | $\mathcal{O}\left(\frac{T^{3/2}\mu^{3/2}}{\sqrt{\epsilon}}N^{5/2}\right)$
Qubitization | - |  $\mathcal{O}$ $\left(T\mu N+\log({1/\epsilon})\right)$ |  $\mathcal{O}$ $\left(T\mu N^{3}+\log({1/\epsilon})N^{2}\right)$
Table 3: Asymptotic scaling of the error, needed steps and number of 2-qubit
operations to evolve a system until $T$ keeping the error below $\epsilon$ as
a function of the number of particles $N$ and for different methods of
propagator decomposition.

It is possible to obtain a more accurate approximation by using a second order
Trotter-Suzuki formula which can be expressed compactly as

$\mathcal{L}_{2}(dt)=\mathcal{L}_{1}\left(\frac{dt}{2}\right)\mathcal{L}^{\dagger}_{1}\left(-\frac{dt}{2}\right)\;.$
(35)

As shown in more detail in Appendix B, one can show that also in this case the
single-step error scales at most linearly with $N$. In particular

$\varepsilon_{2}(dt)\leq\mathcal{O}(dt^{3}\mu^{3}N)\;.$ (36)

In this case, the number of steps needed to guaranteeing a total error
$\epsilon$ for simulation up to $T=rdt$ scales as

$r_{2}=\mathcal{O}\left(\left(T\mu\right)^{3/2}\sqrt{\frac{N}{\epsilon}}\right)\;.$
(37)

Since the number of $SU(4)$ gates is now $N(2N-3)/2$ (assuming, as before, $N$
even) the total gate cost is now

$\mathcal{C}_{2}\leq\frac{N(2N-3)}{2}r_{2}=\mathcal{O}\left(\frac{(T\mu)^{3/2}}{\sqrt{\epsilon}}N^{5/2}\right)\;.$
(38)

The scaling with the number of neutrinos $N$ has substantially improved within
this scheme. The use of even higher order formulas could allow to reach an
almost optimal scaling $\mathcal{C}(N^{2+\delta})$ for $\delta\ll 1$, but with
possibly much larger constant prefactors.

Finally, we want to comment on the prospect of using more modern approaches to
simulate the time evolution operator, such as qubitization [49, 50]. This
scheme also approximates the propagator over small time intervals but,
contrary to Trotter-Suzuki based approaches, is able to reach an optimal
scaling in both $T$ and the error $\epsilon$ for the number of steps

$r_{Q}=\mathcal{O}\left(T\alpha_{H}+\log\left(\frac{1}{\epsilon}\right)\right)\;.$
(39)

Here, $\alpha_{H}$ is a suitable norm of the Hamiltonian operator which for
the two-body neutrino Hamiltonian is given by

$\alpha_{H}=3\frac{\mu}{N}\sum_{i<j}^{N}\left(1-\cos(\theta_{ij})\right)=\mathcal{O}(\mu
N)\;.$ (40)

For a general angular distribution, the gate cost for each segment scales as
the number of different coefficients and thus the gate cost for the algorithm
scales as

$\mathcal{C}_{Q}=\mathcal{O}\bigg{(}T\mu
N^{3}+N^{2}\log(\frac{1}{\epsilon})\bigg{)}\bigg{)}\;.$ (41)

For a fixed evolution time $T$ and target error $\epsilon$ the second-order
Trotter-Suzuki scheme scales then better than a qubitization based approach.
This is not a special property of the neutrino system and has been noted
already in other applications (e.g. simulations of the Schwinger model [51,
52]) and is related to the fact that qubitization does not exploit the
commutation properties of the terms that form the Hamiltonian.

Figure 12: Number of 2-qubit operations $\mathcal{C}$, in function of $N$,
needed to evolve the system of a total time $T=40\,\mu^{-1}$ with an error
$\varepsilon_{1}(T)\leq\epsilon=0.15$. Solid blue line corresponds to first-
order bound and green band to the real scaling achievable using different pair
ordering. The olid orange line refers to the second-order bound instead and
the yellow band is the corresponding real achievable complexity. Using the
optimal decomposition, for the system of $N=4$ neutrinos, the number of
2-qubit universal operation needed is equal to $10N(N-1)/2=60$ and corresponds
to $180$ $CNOT$ or $ZZ$ gates that we actually used to obtains the last point
in Fig. 11.

We summarize the results on the bounds for both the number of time steps and
the circuit (gate) complexity of a simulation with $N$ neutrinos, maximum time
$T$ and error tolerance $\epsilon$ in Tab. 3.

The actual implementation cost is several orders of magnitude lower than what
predicted by the theoretical bounds. Moreover, as demonstrated in Sec. III,
the cost can be further reduced by using a good decomposition of the
propagator which guarantees a smaller error and therefore allows to use a
greater value of time step $dt$. For a total evolution time $T=40\mu^{-1}$ and
target error $\epsilon=0.15$, we show in Fig. 12 the theoretical bounds gate
count $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ for the first and second order
decomposition as solid blue and orange lines respectively. We have also
numerically determined the actual range of gate counts required for this
simulation as we vary the order in the Trotter decomposition and we show it in
Fig. 12 as the shaded green and yellow bands for first and second formulas
respectively. The real complexity is calculated using a linear accumulation
error, that is we search $dt$ such that $r=T/dt$ guarantees:

$r\norm{\widetilde{U}_{2}(T/r)-U_{2}(T/r)}\leq\epsilon=10^{-2}\,.$ (42)

where $\widetilde{U}_{2}(T/r)=\mathcal{L}_{1}(T/r)$ or
$\widetilde{U}_{2}(T/r)=\mathcal{L}_{2}(T/r)$ as defined in Eq.ns (32) and
(35) respectively. Using the exact accumulation of error:

$\norm{\widetilde{U}_{2}(T/r)^{r}-U_{2}(T)}\leq 10^{-2}$ (43)

the number of needed steps could be potentially smaller.

##  VII Conclusion

In this work we considerably extended the work of Ref. [37] in terms of
designing efficient ways to simulate through a quantum computer the time
evolution of a system of several neutrinos interacting with an all-to-all
Hamiltonian, demonstrating that the complexity of the algorithm scales in a
polynomial way with the number of particles. In particular, we have shown that
the choice of an optimal order in the decomposition can decrease the Trotter
error and therefore the number of steps necessary to simulate the evolution of
the system for a certain time. This is possible thanks to the use of a machine
with full connectivity between the qubits. Devices of this type allow to
implement the evolution in circuits with a smaller number of single and
2-qubit gates, potentially decreasing the computational error. The proposed
algorithm was tested on the Quantinuum H1-2 quantum computer to carry out the
evolution of a system of $N=4$ and $N=8$ neutrinos. The device employed sports
very good fidelity gates, and the results obtained, even without error
mitigation, have remarkably small errors even in the largest circuits we
implemented, containing up to and $616$ $SU(2)$ rotations and $342$ two qubit
gates on 8 qubits. In the previous study from Ref. [37], error mitigation
using zero noise extrapolations was instead crucial to obtain acceptable
results, even for circuits with a significantly lower number of operations
(cf. Tab 1 here).

Future improvements could be obtained using various error mitigation
techniques, such as symmetry protection [53], virtual distillation [54] and
symmetry verification [55]. Considering the quality of the results obtained in
this work we conclude that it would be possible both to increase the number of
neutrinos in the system, and to evolve it for a greater total time $T$. In
order to describe more realistic and phenomenologically rich neutrino systems,
it will be important to extend the algorithms presented here in order to
simulate neutrinos with different energies, which can be modeled with a
particle-dependent external field. Furthermore, another issue to investigate
is the description of collective oscillations in the presence of electrons,
thus including the matter part in the Hamiltonian.

###### Acknowledgements.

This research used resources of the Oak Ridge Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under Contract DE-
AC05-00OR22725. F.T. is supported by the Q@TN grant ANuPC-QS. P.L. is
supported by the Q@TN grant ML-QForge. This work was prepared in part by LLNL
under Contract No. DE-AC52-07NA27344 with support from the Laboratory Directed
Research and Development Grant No. 19-DR-005.

## References

  * Feynman [1982] R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
  * Klco  _et al._ [2022] N. Klco, A. Roggero, and M. J. Savage, Reports on Progress in Physics 85, 064301 (2022).
  * Gottesman [1997] D. Gottesman, CalTech Ph.D. Thesis (1997), arXiv:quant-ph/9705052 .
  * Gottesman [1998] D. Gottesman, in  _22nd International Colloquium on Group Theoretical Methods in Physics_ (1998) arXiv:quant-ph/9807006 .
  * Vidal [2003] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
  * Schuch  _et al._ [2008] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 100, 030504 (2008).
  * Deutsch [1985] D. Deutsch, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 400, 97 (1985).
  * Nielsen and Chuang [2001] M. A. Nielsen and I. L. Chuang, Phys. Today 54, 60 (2001).
  * Pantaleone [1992a] J. Pantaleone, Physics Letters B 287, 128 (1992a).
  * Pantaleone [1992b] J. Pantaleone, Phys. Rev. D 46, 510 (1992b).
  * Qian and Fuller [1995a] Y.-Z. Qian and G. M. Fuller, Physical Review D 51, 1479 (1995a).
  * Samuel [1996] S. Samuel, Phys. Rev. D 53, 5382 (1996).
  * Pastor and Raffelt [2002] S. Pastor and G. Raffelt, Phys. Rev. Lett. 89, 191101 (2002).
  * Balantekin and Yüksel [2005] A. B. Balantekin and H. Yüksel, New Journal of Physics 7, 51 (2005).
  * Wu and Tamborra [2017] M.-R. Wu and I. Tamborra, Phys. Rev. D 95, 103007 (2017).
  * Qian and Fuller [1995b] Y.-Z. Qian and G. M. Fuller, Phys. Rev. D 52, 656 (1995b).
  * Qian  _et al._ [1993] Y.-Z. Qian, G. M. Fuller, G. J. Mathews, R. W. Mayle, J. R. Wilson, and S. E. Woosley, Phys. Rev. Lett. 71, 1965 (1993).
  * Janka [2012] H.-T. Janka, Annual Review of Nuclear and Particle Science 62, 407 (2012).
  * Duan  _et al._ [2010] H. Duan, G. M. Fuller, and Y.-Z. Qian, Annual Review of Nuclear and Particle Science 60, 569 (2010).
  * Smirnov [2005] A. Y. Smirnov, Physica Scripta 2005, 57 (2005).
  * Pantaleone [1992c] J. Pantaleone, Physics Letters B 287, 128 (1992c).
  * Duan  _et al._ [2006a] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys. Rev. Lett. 97, 241101 (2006a).
  * Duan  _et al._ [2006b] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys. Rev. D 74, 105014 (2006b).
  * Chakraborty  _et al._ [2016] S. Chakraborty, R. Hansen, I. Izaguirre, and G. Raffelt, Nuclear Physics B 908, 366 (2016).
  * Cervia  _et al._ [2019] M. J. Cervia, A. V. Patwardhan, A. B. Balantekin, S. N. Coppersmith, and C. W. Johnson, Phys. Rev. D 100, 083001 (2019).
  * Rrapaj [2020] E. Rrapaj, Phys. Rev. C 101, 065805 (2020).
  * Patwardhan  _et al._ [2021] A. V. Patwardhan, M. J. Cervia, and A. B. Balantekin, Phys. Rev. D 104, 123035 (2021).
  * Birol  _et al._ [2018] S. Birol, Y. Pehlivan, A. B. Balantekin, and T. Kajino, Phys. Rev. D 98, 083002 (2018).
  * Martin  _et al._ [2022] J. D. Martin, A. Roggero, H. Duan, J. Carlson, and V. Cirigliano, Phys. Rev. D 105, 083020 (2022).
  * Xiong [2022] Z. Xiong, Phys. Rev. D 105, 103002 (2022).
  * Roggero  _et al._ [2022] A. Roggero, E. Rrapaj, and Z. Xiong, arXiv preprint arXiv:2203.02783 (2022).
  * Roggero [2021a] A. Roggero, Phys. Rev. D 104, 103016 (2021a).
  * Roggero [2021b] A. Roggero, Phys. Rev. D 104, 123023 (2021b).
  * Cervia  _et al._ [2022] M. J. Cervia, P. Siwach, A. V. Patwardhan, A. Balantekin, S. Coppersmith, and C. W. Johnson, arXiv preprint arXiv:2202.01865 (2022).
  * Lacroix  _et al._ [2022] D. Lacroix, A. B. Balantekin, M. J. Cervia, A. V. Patwardhan, and P. Siwach, arXiv preprint arXiv:2205.09384 (2022), https://doi.org/10.48550/arxiv.2205.09384.
  * Hall  _et al._ [2021] B. Hall, A. Roggero, A. Baroni, and J. Carlson, Physical Review D 104, 063009 (2021).
  * Yeter-Aydeniz  _et al._ [2022] K. Yeter-Aydeniz, S. Bangar, G. Siopsis, and R. C Pooser, Quantum Information Processing 21, 1 (2022).
  * Illa and Savage [2022] M. Illa and M. J. Savage, arXiv preprint arXiv:2202.12340 (2022), https://doi.org/10.48550/arXiv.2202.12340.
  * Balantekin and Fuller [1999] A. Balantekin and G. Fuller, Physics Letters B 471, 195 (1999).
  * Pehlivan  _et al._ [2011] Y. Pehlivan, A. Balantekin, T. Kajino, and T. Yoshida, Physical Review D 84, 065008 (2011).
  * Childs  _et al._ [2019] A. M. Childs, A. Ostrander, and Y. Su, Quantum 3, 182 (2019).
  * Vatan and Williams [2004] F. Vatan and C. Williams, Physical Review A 69, 032315 (2004).
  * Pino  _et al._ [2021] J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, M. S. Allman, C. H. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, C. Ryan-Anderson, and B. Neyenhuis, Nature 592, 209 (2021).
  * Roggero  _et al._ [2020a] A. Roggero, C. Gu, A. Baroni, and T. Papenbrock, Phys. Rev. C 102, 064624 (2020a).
  * Roggero  _et al._ [2020b] A. Roggero, A. C. Y. Li, J. Carlson, R. Gupta, and G. N. Perdue, Phys. Rev. D 101, 074038 (2020b).
  * Das  _et al._ [2021] P. Das, S. Tannu, S. Dangwal, and M. Qureshi, in  _MICRO-54: 54th Annual IEEE/ACM International Symposium on Microarchitecture_ (2021) pp. 950–962.
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Physical Review X 11, 011020 (2021).
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, Phys. Rev. Lett. 118, 010501 (2017).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * Shaw  _et al._ [2020] A. F. Shaw, P. Lougovski, J. R. Stryker, and N. Wiebe, Quantum 4, 306 (2020).
  * Rajput  _et al._ [2021] A. Rajput, A. Roggero, and N. Wiebe, arXiv preprint arXiv:2109.03308 (2021), 10.48550/ARXIV.2109.03308.
  * Tran  _et al._ [2021] M. C. Tran, Y. Su, D. Carney, and J. M. Taylor, PRX Quantum 2, 010323 (2021).
  * Huggins  _et al._ [2021] W. J. Huggins, S. McArdle, T. E. O’Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean, Physical Review X 11, 041036 (2021).
  * Bonet-Monroig  _et al._ [2018] X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. O’Brien, Physical Review A 98, 062339 (2018).

##  Appendix A First order Trotter error for the interaction Hamiltonian

The 2-body part of the neutrino Hamiltonian in Eq. (3) of the main text is a
sum of $\Gamma=N(N-1)/2$ terms:

$H^{(2)}=\sum_{i<j}^{N}h_{ij}:=\sum_{K=1}^{\Gamma}h_{K}\,,$ (44)

where $h_{ij}=J_{ij}\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}$ and
$J_{ij}=\mu(1-\cos(\theta_{ij}))/N$. Note that we count here the neutrinos
from $1$ to $N$ instead of from $0$ to $N-1$. Using the first order Trotter
decomposition we can implement the propagator with a product formula:

$U_{2}(dt)\approx\mathcal{L}_{1}(dt)=\prod_{i<j}^{N}e^{-ih_{ij}dt}=\prod_{K=1}^{\Gamma}e^{-ih_{K}dt}\,.$
(45)

Using the result from Proposition 9 of Ref. [48] we bound the first order
Trotter error as:

$\begin{split}\epsilon_{1}&=\norm{\mathcal{L}_{1}(dt)-e^{-iHdt}}\\\
&\leq\frac{dt^{2}}{2}\sum_{K=1}^{\Gamma}\norm{\sum_{L=K+1}^{\Gamma}\commutator{h_{K}}{h_{L}}}\,,\end{split}$
(46)

where in our case $K$ and $L$ correspond to pair indices $K=(i,j)$ and
$L=(k,l)$. The sum inside the norm in the expression above can be expressed
explicitly as

$\begin{split}\sum_{L=K+1}^{\Gamma}\commutator{h_{K}}{h_{L}}&=\delta_{ik}\sum_{l=j+1}^{N}\commutator{h_{ij}}{h_{kl}}\\\
&\quad\quad\quad+\sum_{k=i+1}^{N}\sum_{l=k+1}^{N}\commutator{h_{ij}}{h_{kl}}\\\
&=\sum_{l=j+1}^{N}\commutator{h_{ij}}{h_{il}}\\\
&\quad\quad\quad+\sum_{k=i+1}^{N}\sum_{l=k+1}^{N}\commutator{h_{ij}}{h_{kl}}\;,\end{split}$
(47)

where we have separated the sum for $L>K$ in two contributions: those where
the first index of $K$ is the same as the first index of $L$ and those for
which the first index of $L$ is greater than the first index of $K$. The last
contribution can be simplified by splitting the sum around the index $j$ and
realizing that, for the commutator $\commutator{h_{ij}}{h_{kl}}$ to be
different from zero at least one index in $(k,l)$ needs to match an index in
$(i,j)$. The result reads

$\begin{split}\sum_{k=i+1}^{N}\sum_{l=k+1}^{N}\commutator{h_{ij}}{h_{kl}}&=\sum_{k=i+1}^{j-1}\commutator{h_{ij}}{h_{kj}}\\\
&\quad\quad\quad+\sum_{l=j+1}^{N}\commutator{h_{ij}}{h_{jl}}\;.\end{split}$
(48)

The Trotter error over a small time step $dt$ thus reads

$\begin{split}\epsilon_{1}\leq\frac{dt^{2}}{2}\sum_{i<j}^{N}&\left\|\sum_{k=i+1}^{j-1}\commutator{h_{ij}}{h_{kj}}\right.\\\
&\quad+\left.\sum_{l=j+1}^{N}\left(\commutator{h_{ij}}{h_{il}}+\commutator{h_{ij}}{h_{jl}}\right)\right\|\,.\end{split}$
(49)

The commutators between different two-body Hamiltonians can be computed
straightforwardly as

$\begin{split}\commutator{h_{ij}}{h_{ik}}&=J_{ij}J_{ik}\commutator{\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}}{\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{k}}\\\
&=2iJ_{ij}J_{ik}\boldsymbol{\sigma}_{i}\cdot(\boldsymbol{\sigma}_{j}\wedge\boldsymbol{\sigma}_{k})\;,\end{split}$
(50)

where we used $\boldsymbol{a}\wedge\boldsymbol{b}$ to denote the standard
cross product in three dimensions. Using the cyclic permutation equivalence of
the cross product, the fact that the coupling matrix $J_{ij}$ is positive and
the bound
$\norm{\boldsymbol{\sigma}_{i}\cdot(\boldsymbol{\sigma}_{j}\wedge\boldsymbol{\sigma}_{k})}_{\infty}\leq
4$ for all $i,j,k\in[1,N]$ one arrives at

$\varepsilon_{1}\leq
4dt^{2}\sum_{i<j}^{N}J_{ij}\norm{\sum_{l=j+1}^{N}(J_{il}-J_{jl})+\sum_{k=i+1}^{j-1}J_{kj}}\,.$
(51)

For a specific choice of angular distributions the sums can be computed
straightforwardly. However, in order to obtain a general bound on the error we
can introduce $\Theta\coloneqq\max_{i,j}\left[1-\cos(\theta_{ij})\right]$ and
obtain the upperbound

$\varepsilon_{1}\leq
12dt^{2}\mu^{2}\frac{\Theta^{2}}{N^{2}}\binom{N}{3}=\mathcal{O}\left(dt^{2}\mu^{2}N\right)\;,$
(52)

Using a fixed Trotter time step $dt$ and evolving the system until a total
time $T$ using $r=T/dt$ steps the total additive error can be bound by:

$\varepsilon_{1}(T)\leq
r\varepsilon_{1}(dt)=12\frac{T^{2}}{r}\mu^{2}\frac{\Theta^{2}}{N^{2}}\binom{N}{3}=\mathcal{O}\bigg{(}\frac{T^{2}\mu^{2}N}{r}\bigg{)}\,.$
(53)

In order to have a total error less that some target error
$\varepsilon_{1}(T)\leq\epsilon$ we need a number of steps that scales still
linearly with $N$:

$r\leq 12\frac{T^{2}\mu^{2}\Theta^{2}}{\epsilon
N^{2}}\binom{N}{3}=\mathcal{O}\bigg{(}\frac{T^{2}\mu^{2}N}{\epsilon}\bigg{)}\,.$
(54)

##  Appendix B Second order Trotter error for the interaction Hamiltonian

We present here the extension of the analysis presented in Appendix A above to
the second order Trotter formula. The approximation to the propagator
$U_{2}(dt)$ now reads

$U_{2}(dt)\approx\mathcal{L}_{2}(dt)=\prod_{L=\Gamma}^{1}e^{-i\frac{dt}{2}h_{L}}\prod_{K=1}^{\Gamma}e^{-i\frac{dt}{2}h_{K}}\;,$
(55)

where we used the multi-index notation $K=(i,j)$, $L=(k,l)$ and
$\Gamma=N(N-1)/2$ as before. Using the result of Proposition 10 from Ref. [48]
one can bound the second order Trotter error by:

$\begin{split}\varepsilon_{2}\leq&\frac{dt^{3}}{12}\sum_{K}^{\Gamma}\norm{\sum_{L>K}^{\Gamma}\sum_{M>K}^{\Gamma}\commutator{h_{L}}{\commutator{h_{M}}{h_{K}}}}\\\
&+\frac{dt^{3}}{24}\sum_{K}^{\Gamma}\norm{\sum_{L>K}^{\Gamma}\commutator{h_{K}}{\commutator{h_{K}}{h_{L}}}}\;.\end{split}$
(56)

In order to bound the second term, we use a similar procedure to the one
adopted in the first order case by expanding the sums and keeping
contributions $\commutator{h_{ij}}{\commutator{h_{ij}}{h_{kl}}}$ with one of
the $(k,l)$ indices matching one of the $(i,j)$ indices. The bound can be
found by using then the expression for the nested commutator

$\begin{split}\commutator{h_{ij}}{\commutator{h_{ij}}{h_{ik}}}&=J_{ij}^{2}J_{ik}\commutator{\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}}{\commutator{\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}}{\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{k}}}\\\
&=-4J_{ij}^{2}J_{ik}\boldsymbol{\sigma}_{i}\cdot(\boldsymbol{\sigma}_{j}\wedge(\boldsymbol{\sigma}_{j}\wedge\boldsymbol{\sigma}_{k}))\,,\end{split}$
(57)

together with the bound
$\norm{\boldsymbol{\sigma}_{i}\cdot(\boldsymbol{\sigma}_{j}\wedge(\boldsymbol{\sigma}_{j}\wedge\boldsymbol{\sigma}_{k}))}\leq
8$. The results for the second term in Eq. (56) reads

$\frac{dt^{3}}{24}\sum_{K}^{\Gamma}\norm{\sum_{L>K}^{\Gamma}\commutator{h_{K}}{\commutator{h_{K}}{h_{L}}}}\leq
4dt^{3}\frac{\mu^{3}\Theta^{3}}{N^{3}}\binom{N}{3}\;.$ (58)

To determine an upper bound for the first term in Eq. (56) we used the same
procedure for both sums inside the norm. In this way four terms are obtained,
each of which can be bounded with triangular inequalities by estimating the
number of non-zero terms. After some calculations the final result reads

$\frac{dt^{3}}{12}\sum_{K}^{\Gamma}\norm{\sum_{L>K}^{\Gamma}\sum_{M>K}^{\Gamma}\commutator{h_{L}}{\commutator{h_{M}}{h_{K}}}}\\\
\leq\frac{dt^{3}}{12}32\frac{\mu^{3}}{N^{3}}\Theta^{3}\left(21\binom{N}{4}+6\binom{N}{3}\right)\\\
$ (59)

So the total error reads

$\begin{split}\varepsilon_{2}&\leq
dt^{3}\frac{\mu^{3}\Theta^{3}}{N^{3}}\bigg{[}20\binom{N}{3}+56\binom{N}{4}\bigg{]}\\\
&=\mathcal{O}\bigg{(}dt^{3}\mu^{3}N\bigg{)}\,,\end{split}$ (60)

showing a linear scaling with $N$. The upper bound error after $r=T/dt$ steps
is

$\begin{split}\varepsilon_{2}(T)&\leq\frac{T^{3}}{r^{2}}\frac{\mu^{3}\Theta^{3}}{N^{3}}\bigg{[}20\binom{N}{3}+56\binom{N}{4}\bigg{]}\\\
&=\mathcal{O}\bigg{(}\frac{T^{3}}{r^{2}}\mu^{3}N\bigg{)}\end{split}$ (61)

The total Trotter steps $r$ needed to evolve up to $T$ keeping the total error
under $\epsilon$ scales as

$\begin{split}r&\leq\frac{(T\mu\Theta)^{3/2}}{\sqrt{\epsilon}N^{3/2}}\sqrt{20\binom{N}{3}+56\binom{N}{4}}\\\
&=\mathcal{O}\bigg{(}\frac{T^{3/2}\mu^{3/2}\sqrt{N}}{\sqrt{\epsilon}}\bigg{)}\,.\end{split}$
(62)