{'fraction_non_alphanumeric': 0.03366041224549829, 'fraction_numerical': 0.013237240770701575, 'mean_word_length': 4.398593466424682, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Trapped-Ion Quantum Simulation of Collective Neutrino Oscillations

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., Ref. [3] for a recent review). Given a quantum system in some fixed initial state |Ψ 0 , its time evolution under the Hamiltonian H is given by the action of the real-time evolution operator: transforming the state according to the time-dependent Schrödinger equation |Ψ(t) = U (t) |Ψ 0 . In general, a direct approach based on this description faces an exponentially growing cost on classical computers as the size of the system increases due to both the enormous memory requirements to encode the states of the system and to the 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 |ϕ 0 , a quantum algorithm consists of a unitary transformation U which produces a certain final state |ϕ f = U |ϕ 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, prob-ability 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 ( 8 M ) 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 (∼ 10 53 erg) in the form of a large number of neutrinos (∼ 10 58 ). In a few seconds the 10% of the gravitational mass of the star is converted into neutrinos flowing with an energy E ν (10 ÷ 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 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 distances from the core (∼ 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 one to study large-scale systems with complex geometries [22][23][24][25]. A full treatment of correlation effects in the complete manybody 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]. Semiclassical methods are another class of approaches that retain some of the correlations while maintaining numerical efficiency; these methods were 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 four 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 operations 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 machineaware compilation that has to consider the topology of the qubit system. We note that to obtain physical information about the many-neutrino system one does not need to simulate the evolution of all the 10 58 emitted neutrinos, but it is sufficient to limit the simulation only to a number of neutrinos in a space region that is causally connected. Furthermore, collective oscillations can operate on timescales much smaller than the total propagation time needed for a neutrino to leave the system starting from the neutrino sphere. Finally, even smallscale simulations can provide crucial information about the conditions required for specific collective modes to appear as well as the nature of correlations responsible for them (see, e.g., Refs. [28,30,[32][33][34]).
In Sec. II we present the description of the physical system of several neutrinos used in this work. In Sec. III we 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 pairwise neutrinos scattering, while at large radii the interaction with electrons prevails instead. Moreover, we assume a narrow cone of forward peaked emitted neutrinos.
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 to find the best quantum circuit for the trappedion 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 multistep 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 ν e and a single heavy flavor ν x which is a combination of µ and τ neutrinos. This can be justified if the mixing angle θ 13 = 0 as shown in [40]. In this way the general flavor state of each neutrino, given by the two-flavor superposition |Ψ = α |ν e + β |ν x , can be fully encoded in the state of a qubit by means of the mapping: The flavor Hamiltonian in this basis can be decomposed into three main terms [20,41]: (1) a one-body contribution describing vacuum mixing H vac := 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 one-body part describing the coupling to external matter leading to the MSW effect and (3) the neutrino-neutrino interaction term H νν := H (2) generated by forward scattering. In this work we 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 [41] where we used bold symbols to denote three-dimensional vectors. The vectors i ) are formed by the Pauli matrices acting on the ith neutrino. Here and in the following we suppress the identity operators acting on the other spins, for example, σ The simple structure of the neutrino-neutrino interaction term originates from treating neutrinos as plane waves with definite momentum and the weak interaction as a contact term in coordinate space. In order to extend this formulation to take into account, in a more realistic way, the spatial localization of neutrinos, especially important in inhomogeneous systems [42], it is possible to consider instead the evolution of neutrino wave packets. Extensions along these lines are left for future work. In the flavor basis, the vector b in the first term describes vacuum mixing of neutrinos with the same energy and is given explicitly by In this expression δm 2 = m 2 2 − m 2 1 is the square mass difference between mass eigenstates, which is of order 10 −4 eV 2 , θ ν is the mixing angle that we took equal to θ ν = 0.195, and E is the neutrino energy. The twobody 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, The coupling constant of the pair interaction can be written explicitly as follows,  and it depends on the relative angle of propagation, where p i is the momentum of the ith 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 µ = √ 2G F n ν , where G F if the Fermi constant, n ν = 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, and measure time in units of µ −1 . As done in previous work [37], we take a simple grid of angles, meant to reproduce a narrow cone of forward peaked neutrinos in accordance with the geometry displayed in (3) turns out to be symmetric under the particle exchange for k = 0, . . . , 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 |ν e flavor state and the other N/2 in the |ν 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 One can obtain the exact time evolution by directly performing a matrix multiplication |Ψ(t) = U (t) |Ψ 0 . In the spin basis, the flavor content of an individual neutrino is obtained from the expectation value of the Pauli matrix σ 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 We display the exact evolution of both quantities in Fig. 2: The top panel shows results for 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 |Ψ 0 , the flavor evolution of neutrinos ν 0 and ν 1 is 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). Because of the presence of this symmetry, in the rest of this work we 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. Because of the two-flavor approximation presented in Sec. II, the flavor state of neutrinos can be directly encoded into a qubit according to the map in Eq. (2). In the case of 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 encoding the flavor state |Ψ 0 is evolved under a sequence of unitary transformations, overall implementing the realtime propagator. A projective measurement of the final state of the qubits eventually allows us to extract flavor observables from the simulation. In the computational basis, the operator U (t) is represented by a 2 N × 2 N unitary matrix which must be decomposed into the gate set provided by the machine, usually composed by single-or two-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 one-body term and expresses the total Hamiltonian as a sum of two-body terms, namely, Then, the total propagator can be approximated by the product of pair propagators: The approximation introduces an error of order O(t 2 ) due to the noncommutativity of the symmetrized two-body terms, namely, [H ij , H ik ] = 0. A first improvement of this implementation can be made by considering that the entire one-body and two-body terms commute, namely, H (1) , H (2) = 0, and therefore their separation does not introduce any errors: Subsequently, the two-body term alone can be efficiently implemented as a pair decomposition: 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 and can cause an explicit breaking of the symmetry under particle exchange. However, the approach can still be useful in the case of a particle-dependent external field, in which b → b i , as it occurs when different neutrinos have different energies E → E i (a necessary ingredient to observe spectral splits [20,28]).

A. 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 at least once during the simulation. This fact . SWAP network for N = 4 qubits that implements the two-body propagator U2(t) using a chain of linearly connected qubits. Vertical lines represent the qubit interaction and crosses the SWAP operation. The pair ordering influences the total error; here, we use the map (q0, q1, q2, q3) = (ν0, ν2, ν1, ν3).
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 (SN)-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 where in the computational basis the SWAP unitary is For example, the SWAP network required for the case of N = 4 is depicted in Fig. 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 [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 us to maximize the cancellations between the commutators minimizing the decomposition error. Note that, once the four-layer 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 or additional layers but at the cost of increasing the complexity in terms of depth and number of two-qubit gates of the scheme. For our Hamiltonian, and in the case of N = 4, we find that the best interaction order would be the one described by the Figure 4. Implementation of the two-body propagator U2(t) on a qubit system with all-to-all connectivity. This scheme can implement the optimal ordering of pairs to minimize the decomposition error. For our system this corresponds to the qubit to neutrino map (q0, q1, q2, q3) = (ν0, ν1, ν2, ν3).
network in Fig. 4. Such ordering cannot be expressed using a SWAP network with four layers (as shown in Fig. 3) if restricted to the use of linear qubit connectivity.
In fact, one can easily show that at least five layers would be needed. With all-to-all connectivity, however, this algorithm can be implemented using only three layers, as shown in Fig. 4. 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. Determining the optimal ordering for large systems is not feasible, in general, as this would require a superexponential cost in the system size N . For large systems a randomization procedure for the order could prove valuable to control the error [43,44]. We analyze the effect of different orderings on a single time step in Fig. 5. A similar study, made for a different Hamiltonian, about the Trotter error dependence on the ordering can be found in Ref. [45]. The two top panels show the evolution of the inversion probability as a function of the time step dt ∈ [0, 40] µ −1 for the same initial state as in Fig. 2 and for different implementation of the propagator: The dotted curve is the exact evolution, and the solid blue line is the one obtained by applying the implementation proposed in this work [Eqs. (16) and (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 SN scheme from Fig. 3. Panel (a) shows results for neutrinos ν 0 and ν 3 while panel (b) shows those for ν 1 and ν 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. Furthermore, the breaking of exchange symmetry in the N = 4 system for the (SN) results shown in Fig. 5 is to be ascribed to the Trotter breakup from Eq. (15), which for the (SN) ordering gen- for ν1 and ν2. The dotted black curve is the exact evolution, and the solid blue line is the one obtained by applying the implementation proposed in this work [Eqs. (16) and (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.
erates an explicit symmetry breaking. 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 spectral norm (i.e., the maximum singular value of the matrix) The results displayed in the top panels show that the error in the inversion probability is lower for the (SN) approximation when dt < 30 µ −1 . This is mostly an effect of choosing a particular initial state |Ψ 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. 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  time, we sequentially apply the approximate propagator where U 2 (dt) is defined in Eq. (17). We employ two different orderings: (1) the best possible ordering that is achievable with a SN with linear connectivity as presented in Fig. 3 above, and denoted by the dashed orange curve in Fig. 6, and (2) the optimal ordering, 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 us 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 one-body part is trivial because it is the tensor product of the same single-qubit gate applied to each qubit: The two-body part of the propagator, as approximated in Eq. (17), is the product of pair terms of the form 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. [46], the optimal CNOT-based decomposition for the SU (2) invariant unitary operator u ij in Eq. (23) can be written as where we introduced the angle φ = −(2dtJ ij + π/2). 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 us to reduce the number of gates necessary for the simulation. In this work we simulate the real-time evolution of N = 4 and N = 8 neutrino systems using the QSM H1-2 trapped-ion device, similar to the machine described in Ref. [47]. Like other trapped-ion-based quantum devices, it offers the advantage of a full connectivity between the qubits, high fidelity: a single-qubit error of ∼ 10 −4 and a two-qubit gate error of ∼ 10 −3 , and the possibility to perform parallel operations between different qubits in different zones simultaneously. The native single-qubit gates used are as follows: and: Currently, only two possible values θ = π, π/2 are natively available on the device. This means that up to four physical gates are required to compile a general SU (2) unitary. The two-qubit gate is In the circuit from Eq. (24), we can replace the following ZZ-based implementation of the CNOT gate: Then, using (1) commutation properties, (2) all possible simplifications including cancellation and gate mergers, and (3) translating all single-qubit gates into native gates, we obtain the circuit displayed in Fig. 7 that implements the two-body part of the evolution for one pair of neutrinos using only native gates. The angle parameter α is defined, for every qubit pair, as α = −dtJ ij . A single Trotter step is implemented by first applying the one-body part in Eq. (22), made up of single-qubit gates, followed by the optimal ordered sequence of pair interactions. A scheme for the circuit of a complete single Trotter step is presented in Fig. 8 where each twoqubit 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 full one-body part to the beginning. In these cases the overall gate count for the circuit can be further optimized using an inverted interaction order in alternate steps [this 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. This alternating scheme can also reduce the overall approximation error since it then becomes equivalent to a second order Trotter step with time step 2dt (see Sec. VI for additional details).

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: First, 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 multistep simulation fixing dt = 4 µ −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-and 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 qubits need to be employed. With the Pegasus topology implemented in the D-Wave Advantage system, the N = 4 simulation presented in Ref. [39] required ≈ 2000 physical qubits; by extrapolating the resource requirements estimated in Fig. 9 there, a simulation with N = 8 neutrinos will require ≈ 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].
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 its 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 trappedion device has more than four qubits at our disposal, in order to minimize the execution cost of the N = 4 experiments, we always use two sets of four 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 eight qubits with a two-qubit gate depth given by either 9 or 21 (for the N = 4 and N = 8 simulations, respectively). Figure 9 shows the results for the inversion probability obtained for the simulation of a system of N = 4 neutrinos while Fig. 10 shows those 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 ν 0 and ν 7 in panel (a), ν 1 and ν 6 in panel (b), ν 2 and ν 5 in panel (c), and ν 3 and ν 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 almost always respect this symmetry within a confidence interval of 68% and always in the interval of 90%. Because of our choice of using eight 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 three parallel executions over 12 qubits using time steps dt = 16, 24, 32 µ −1 , we found them to be compatible with the results shown in Fig. 9, indicating that Figure 8. Scheme for implementing a single Trotter step. We first apply the single-qubit gates corresponding to the onebody propagator U1(dt) and then the pair propagator using the optimal ordering implementing U2(dt).  cross-talk effects are minimal.
For N = 8 neutrinos our implementation of the propagator respects the symmetry of particle exchange up to dt ≈ 24 µ −1 as can be seen by the theoretical results shown as solid lines in Fig. 10. The real data respect 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 Ref. [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 superconducting device [37], in Table I we present the values for χ 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 following function:   where δP i (t) is the estimated variance (taken to coincide with the 68% confidence interval) and 10 is the number of simulated points, for each neutrino, used to calculate χ 2 . The last third 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 onequbit rotations afforded by the different Trotter decomposition from Eq. (16) adopted here. The number of general SU (2) single-qubit operations 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 three 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 π/2. This helps reduce  Table I. Values of χ 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.
the effect of coherent errors in the final results. Finally, the results of the simulations carried out on IBMQ are obtained using a much larger statistical sample (8192 circuit repetitions instead of 200), in order to more directly compare the new results with the ones obtained there; in the last line of Table I, we also report the estimated χ 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 an estimate is also affected by systematic errors (for more details see Refs. [37,48,49]). However, the strong effect of an increased gate fidelity is still evident.

B. 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 |Ψ 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 For the results shown in this section, we used k ∈ [1, 10] with a time step dt = 4 µ −1 . Also in this case, in order to exploit the specific features of the trapped-ion machine, we have carried out the simulations in pairs in which the first four neutrinos are evolved up to a time T = kdt and the last four 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 four qubits and 36 to the last four. We summarize the gate counts needed to simulate the system up to a certain evolution time T = kdt, divided into the number of two-qubit ZZ gates and one-qubit SU (2) unitaries, in Table II. For each step we need 6 × 3 two-qubit gates while, exploiting cancellations between neighboring pair propagators, the number of general single-qubit rotations needed scales as 32 × k + 4. 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 four qubits are subject to idle errors which, in the case of trapped ions, are mainly related to dephasing. One way to improve the simulation is to use a Dynamical Decoupling approach [50] consisting in applying to idle qubits a sequence of single-qubit gates that have the overall effect of an identity.  Table II. Number of SU (2) gates and ZZ gates to evolve the system to a fixed time T = kdt to produce the results in Fig.  11.

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 ν = N/V and 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 twobody part of the propagator U 2 (dt) = e −iH (2) dt using a standard first order Trotter formula: By leveraging the results in Ref. [51] (see Appendix A for additional details), we can bound the spectral norm error in this approximation as Using the union bound to obtain ε 1 (T ) ≤ rε 1 (dt) with dt = T /r, one also finds that the total number of steps r 1 required to evolve for a final time T while keeping the total error below scales linearly with system size. We can then bound the gate cost C 1 of a quantum circuit by implementing the full evolution in terms of the number of general SU (4) two-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 the gate cost using a first order product formula scales with Since SU (4) transformations can be decomposed with at most three entangling gates [46], 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 3C 1 . It is possible to obtain a more accurate approximation by using a second order Trotter-Suzuki formula which can be expressed compactly as As shown in more detail in Appendix B, one can also show that in this case the single-step error scales at most linearly with N . In particular, In this case, the number of steps needed to guarantee a total error for simulation up to T = rdt scales as Decomposition type Single-step error Number of steps Circuit complexity Table III. Asymptotic scaling of the error, needed steps, and number of two-qubit operations to evolve a system until T , keeping the error below as a function of the number of particles N and for different methods of propagator decomposition.
Since the number of SU (4) gates is now N (2N − 3)/2 (assuming, as before, N even) the total gate cost is now The scaling with the number of neutrinos N has substantially improved within this scheme. The use of even higher order formulas could allow one to reach an almost optimal scaling C(N 2+δ ) for δ 1, but with possibly much larger constant prefactors.
Finally, we comment on the prospect of using more modern approaches to simulate the time evolution operator, such as qubitization [52,53]. 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 for the number of steps Here, α H is a suitable norm of the Hamiltonian operator which, for the two-body neutrino Hamiltonian, is given by 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 For a fixed evolution time T and target error , the second order Trotter-Suzuki scheme then scales better than a qubitization-based approach. This is not a special property of the neutrino system; it has been noted already in other applications (e.g., simulations of the Schwinger model [54,55]) and is related to the fact that qubitization does not exploit the commutation properties of the terms that form the Hamiltonian. 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 in Table III.
The actual implementation cost is several orders of magnitude lower than what is 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 us to use a greater value of time step dt. For a total evolution time T = 40 µ −1 and target error = 0.15, we show in Fig. 12 the theoretical bounds' gate count C 1 and C 2 for the first and second order decompositions as solid blue and orange lines, respectively. We also numerically determine the actual range of gate counts required for this simulation as we vary the order in the Trotter decomposition, which we show 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 where U 2 (T /r) = L 1 (T /r) or U 2 (T /r) = L 2 (T /r) as defined in Eqs. (32) and (35), respectively. Using the exact accumulation of error, the number of needed steps could potentially be 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 us to implement the evolution in circuits with a smaller number of single-and two-qubit gates, potentially decreasing the computational error. The proposed algorithm was tested on the Quantinuum H1-2 quantum with an error ε1(T ) ≤ = 0.15. The solid blue line corresponds to the first order bound and the green band to the real scaling achievable using different pair ordering. The solid 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 two-qubit universal operations needed is equal to 10N (N − 1)/2 = 60 and corresponds to 180 CNOT or ZZ gates, which we actually use to obtain the last point in Fig. 11. computer to carry out the evolution of a system of N = 4 and N = 8 neutrinos. The device employed had 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 616 SU (2) rotations and 342 two-qubit gates on eight 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. Table I here).
Future improvements could be obtained using various error mitigation techniques, such as symmetry protection [56], virtual distillation [57], and symmetry verification [58]. Considering the quality of the results obtained in this work, we conclude that it would be possible to both 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.
where h ij = J ij σ i · σ j and J ij = µ(1 − cos(θ 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: Using the result from Proposition 9 of Ref. [51] we bound the first order Trotter error as: 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 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 [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 commutators between different two-body Hamiltonians can be computed straightforwardly as where we used a∧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 σ i · (σ j ∧ σ k ) ≤ 4 for all i, j, k ∈ [1, N ] one arrives at (A8) 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 Θ := max i,j [1 − cos(θ ij )] and obtain the upperbound 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: (A10) In order to have a total error less than some target error ε 1 (T ) ≤ we need a number of steps that scales still linearly with N : 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 where we used the multi-index notation K = (i, j), L = (k, l) and Γ = N (N − 1)/2 as before. Using the result of Proposition 10 from Ref. [51] one can bound the second order Trotter error by: 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 [h ij , [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 together with the bound σ i · (σ j ∧ (σ j ∧ σ k )) ≤ 8. The results for the second term in Eq. (B2) reads To determine an upper bound for the first term in Eq. (B2) 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 So the total error reads ε 2 (dt) ≤ dt 3 µ 3 Θ 3 N 3 20 showing a linear scaling with N . The upper bound error after r = T /dt steps is The total Trotter steps r needed to evolve up to T keeping the total error under scales as (B8)