Collective neutrino oscillations with tensor networks using a time-dependent variational principle

If a system of flavor-oscillating neutrinos is at high enough densities that neutrino-neutrino coherent forward scatterings are non-negligible, the system becomes a time-dependent many-body problem. An important and open question is whether the flavor evolution is sufficiently described by a mean-field approach or can be strongly affected by correlations arising from two-body interactions in the neutrino Hamiltonian, as measured by nontrivial quantum entanglement. Numerical computations of the time evolution of many-body quantum systems are challenging because the size of the Hilbert space scales exponentially with the number of particles N in the system. Thus, it is important to investigate approximate but beyond-mean-field numerical treatments at larger values of N. Here we investigate the efficacy of tensor network methods to calculate the time evolution of interacting neutrinos at larger values of N than are possible with conventional methods. In particular, we introduce the use of time-dependent variational principle methods to address the long-range (in momentum space) interactions of the neutrino Hamiltonian when including many distinct vacuum oscillation frequencies. We also define new error measures based upon the instantaneously conserved charge operators known for this Hamiltonian to determine validity of large-N tensor network calculations.


I. INTRODUCTION
Collective effects in the flavor oscillations of neutrinos in environments where large fluxes of neutrinos are present, such as core collapse supernovae, neutron star mergers, or the early universe, have been a subject of great interest over the past few decades (e.g., [1][2][3][4] and references therein). Such flavor oscillations of neutrinos could play an important role in the synthesis of elements in these environments, as well as in the supernova explosion mechanism itself [5][6][7][8][9][10][11][12][13][14][15]. Understanding collective flavor oscillation effects is needed to robustly interpret numerical simulations of these environments.
Despite the weakness of weak interactions, at sufficient density neutrino-neutrino interactions contribute substantially to the neutrino forward scattering potential, transforming collective neutrino oscillations into a quantum many-body problem. As in any interacting manybody quantum system, the dimension of the Hilbert space describing the wave function of the system grows exponentially with particle number. Consequently, the computational complexity grows exponentially as well, and it is untenable to fully solve the interacting many-neutrino system for the large densities of neutrinos present in environments where collective effects matter.
To get around this roadblock, one frequently turns to "mean-field" treatments which neglect multineutrino quantum correlations [16][17][18][19][20][21]. Assessing the reliability of the mean-field approximation in this context has been a topic of long-standing interest, explored through the use of simplified models of interacting neutrinos . Here in this paper, we conduct a comparison of advanced numerical techniques for time evolution of many-neutrino systems, and we further explore whether such time evolution brings about strong many-neutrino correlations, i.e., entanglement, signaling a deviation from mean-field approaches. We find that tensor network methods, described below, can provide a significant speedup, allowing us to reach much larger values of N for certain initial conditions. At these larger values of N , our simulations continue to find significant entanglement in the manyneutrino system.

A. Overview of past and present numerical approaches
As in many-body problems more generally, the baseline approach for collective neutrino flavor oscillations is a mean-field model, replacing two-body interactions with an average one-body potential [16][17][18][19][20][21]. To study the beyond-mean-field time evolution of interacting neu-arXiv:2202.01865v4 [hep-ph] 23 Jun 2022 trinos, the system (in the two-flavor, single angle approximation) was mapped [31,33] to the Richardson-Gaudin Hamiltonian, which was originally developed as a model of pairing solvable by the Bethe ansatz and which has since been applied to a variety of many-body systems, including atomic nuclei [42].
For neutrino numbers N ≥ 10, however, numerical solutions of the Bethe ansatz for time evolution of the many-neutrino system were empirically unstable. As a result, the authors of Ref. [39] instead utilized a fourthorder Runge-Kutta (RK4) method to integrate the timedependent Schrödinger equation for the N -body neutrino wave function. By using sparse-matrix representations of the operators constituting the Hamiltonian, neutrino systems with N up to 16 could be studied. The exponential scaling of the computational complexity eventually renders this method intractable for larger numbers of neutrinos. On the other hand, this work made a potentially useful observation, namely that the degree of quantum entanglement seemed to be larger among the neutrinos nearer to the "spectral splits" in their energy distributions, and smaller among the neutrinos away from these splits. This suggests that a numerical scheme that can zero in on specific subregions of the full Hilbert space where the entanglement primarily resides could potentially yield accurate results while scaling more favorably with N , compared to traditional integration methods.
Given the limitations of both the Bethe ansatz and direct RK4 approaches, in this paper we turn to the use of tensor network methods to model correlated neutrino wave functions and to investigate the dynamics of collective neutrino oscillations. Tensor networks provide a means for calculating dynamics using a truncated basis set with dimensions that can grow subexponentially with system size but that can nonetheless be highly entangled. In this approach, the many-body wave function is written in terms of inner products of tensors that encode pairwise entanglement.
The problem of collective neutrino oscillations involves nonlocal (in momentum space) interactions, and so it is nontrivial to determine whether the method is wellbehaved in such a way that one may practically apply these recent methods without requiring exponential growth in the 'bond dimension' [as defined in Eq. (14)] used to forward-integrate the wave function. In order to treat the nonlocal Hamiltonian, one can implement the tensor network using SWAP operations (analogous to SWAP gates in quantum computing, both defined in Ref. [35]) to "localize" the interactions [35,36,41]: a nonlocal interaction is replaced by interactions between network neighbors interlaced with the SWAP operations that make distant network neighbors appear nearby. Tests of this method have been limited, however, to just a few different neutrino momenta, typically while employing the two-beam model [43].
While intriguing properties such as phase transitions can still be learned in a model with a reduced set of neutrino momenta, including many neutrino modes leads to additional effects, such as spectral splits, which we expect to be physically relevant to astrophysical phenomena in these environments. To address these issues, in this paper we calculate the time evolution of a tensor network model of the many-neutrino wave function using recent time-dependent variational principle (TDVP) methods [44,45]. We compare our reduced-basis, tensornetwork model against two methods calculating the entire wave function: fourth-order Runge-Kutta and Lanczos propagation. (Lanczos is a different kind of reducedbasis method: while the underlying basis dimension is not changed, Lanczos efficiently computes time evolution by projecting into a small effective basis dictated by the initial state and by a few iterated applications of the Hamiltonian.) Furthermore, since these numerically exact methods use the entire basis and thus scale exponentially with the system size N , there is a limit in size after which the accuracy of tensor network methods cannot be assessed by comparison with other methods; as such, we introduce a strategy for consistency checks in our evolved wave function to guide our calculations for N beyond the abilities of RK and Lanczos methods on modern hardware. We find that tensor network calculations are potentially very useful for the study of the collective neutrino oscillation problem. In particular, for initial conditions of a spectrum of neutrinos that result in fewer different spectral split frequencies, the entanglement of the system can be efficiently described by a matrix-product state (MPS) representation of our many-body neutrino state, permitting memory-and time-efficient computations of time evolution. As depicted in Fig. 1, we find that for an initial condition with just one spectral split (i.e., a system starting with all electron-flavor neutrinos), the basis set for representing our MPS wave function can be reduced greatly, allowing for improvements in complexity. In contrast, in the case of a system with two spectral splits, though for sufficiently large N we find some speedup, the required basis set to obtain accuracy comparable to exact methods is much greater in size. We conclude the usefulness of this treatment of neutrinos depends considerably on the number of spectral splits that result from an initial condition. Nonetheless, computing the time evolution of the many-neutrino systems using a timedependent variational principle on a tensor networkthe central approach of this paper-is a promising tool for modeling and understanding the beyond-mean-field behavior of collective flavor oscillations.

B. Outline of the paper
The paper is organized as follows. In Sec. II, we introduce our toy model of collective neutrino oscillations in a dense neutrino gas. In Sec. III, we briefly summarize the methods for calculating the entire wave function after time evolution from the time-dependent Hamiltonian describing our problem, and we introduce the recent tensor network techniques for treating the same problem. In Sec. IV, we define a measure for error in our evolved wave function calculated with any method, based upon the instantaneously conserved charges of the Hamiltonian from Ref. [26], to assess how precisely the state solves the Schrödinger equation. In Sec. V, we use this new measure to assess the quality of a given calculation and decide upon appropriate time-step sizes as well as bond dimension values set, depending upon the initial condition chosen for our problem. In Sec. VI, we summarize our findings and suggest paths forward to investigate larger systems within our model. Finally, in Appendix A we provide additional details about the calculations of entire wave functions, while in Appendix B, we explain in greater detail the tensor network methods to perform time evolution with long-range interactions.

II. DEFINITIONS
To study collective neutrino oscillations from a manybody perspective in a two-flavor, single-angle model, we consider the Hamiltonian [21, 26, 29-31, 33, 39, 46] where ω denotes the vacuum oscillation frequencies of the neutrinos and µ(t) is the time-dependent, angle-averaged ν-ν interaction strength. Here, we have used the SU (2) neutrino isospin operators (in the mass basis) J ω for a neutrino of a given mode ω: with c † iω and c iω as the (fermionic) creation and annihilation operators of a neutrino of mass eigenstate |ν i in the mode ω = ∆m 2 /(2|p|), where ∆m 2 is the masssquared difference and p is the momentum of this mode. Each neutrino in this model therefore has a description as a plane wave with well-defined momentum. Such a treatment is considered adequate for capturing the coherent many-body effects in collective neutrino oscillations [24,25].
In the single-angle approximation (e.g., Ref. [47]), all the time dependence of the Hamiltonian is described by a single, angle-averaged parameter µ(t). Consequently, the flavor evolution of a neutrino in this approximation depends only on its energy and not on the direction of its momentum, reducing the computational complexity of the problem. The single-angle approximation is known to exhibit many of the same collective phenomena, such as synchronized precession and spectral swaps/splits, that are also present in the more sophisticated, multiangle treatments [48].
The many-body state |Ψ(t) of a N -neutrino system then evolves according to the Schrödinger equation 1 with the time-dependent Hamiltonian in Eq. (1): The polarization vector for a neutrino with a given ω is defined as 1 Note that there is just one affine parameter in the evolution equation, if the neutrinos are assumed to be relativistic and their emission from the source is assumed to be time-independent. The latter assumption can be justified based on an observed hierarchy in the dynamical timescales. In a core-collapse supernova environment, a neutrino would experience significant interactions with other neutrinos over an interval of 1000 km [2], or equivalently, a timescale of O(ms). This is much smaller than the O(s) timescales over which the emission characteristics (luminosities and energy spectra of different flavors) change significantly in the late-time neutrino-driven wind phase; see, e.g., the estimates in Refs. [49][50][51] for the cooling and deleptonization timescales. Throughout the text, the affine parameter is henceforth interchangeably referred to as either time (t) or radius (r).
where J ω is the corresponding isospin operator for that neutrino and Ψ is the many-body wave function of the N -neutrino system. The polarization vectors can also be obtained from the Pauli spin decomposition of the individual neutrino density matrices in the mean-field limit, or the "reduced" density matrices in the case of a many-body calculation. In each case, the decomposition is given by where 1 is the 2 × 2 identity matrix and σ j are the Pauli spin matrices. From the above expressions, one can also conclude that the probability of finding an individual neutrino in the mass eigenstate |ν 1 is i.e., the 11 matrix element of the (reduced) density matrix, ρ(ω) = Tr ω ( =ω) [|Ψ Ψ|]. Then, entanglement entropy between a neutrino with frequency ω and the rest of the ensemble can be obtained from where λ ± (ω) = 1 2 (1 ± | P (ω)|) (9) are the eigenvalues of the reduced density matrix ρ(ω).

III. METHODS
As per the bulb model [2], we take a system composed of neutrinos in definite flavor states emitted isotropically from a source; in this case, the initial many-body state has the form |Ψ = N j=1 |ν αj , where α j = e or x for each i, and evolves according to Eq. (4) with the timedependent Hamiltonian in Eq. (1). Additionally, we consider a "box spectrum" with discrete, equally spaced vacuum oscillation frequencies ω j = jω 0 , for j = 1, . . . , N (where ω 0 is an arbitrary reference frequency), such that each oscillation frequency is occupied by a single neutrino. Similarly in keeping with Ref. [33], we use the mixing angle θ = 0.161, a single-angle coupling with R ν = 32.2 ω −1 0 and µ(R ν ) = 3.62 × 10 4 ω 0 , and an initial time/radius given by µ(r) = 5 ω 0 .
Moreover, in the mean-field treatment, the evolution of the N -body neutrino system can then be described using a set of N differential equations, each describing the evolution of one neutrino. In terms of the polarization vectors P (ω), 2 the evolution equations can be written as (11) where in the mass basis B = (0, 0, −1).

B. Numerical calculations of time evolution of many-body wave functions
One can of course directly solve the time-dependent Schrödinger Eq. (4) in the Hilbert space for N neutrinos spanned by 2 N basis states of the form ω |ν αω , where each |ν αω is a flavor-spinor for neutrino frequency ω (e.g., in the flavor basis, each ν αω = ν e or ν x , resulting in 2 N combinations over N frequencies). While limited by the exponential growth of the basis dimension, this is nonetheless a useful benchmark for other methods. While large in dimension, because the Hamiltonian from Eq. (1) has interactions at most between two flavor-spinors, the matrix representation ofĤ in the e, x basis (equivalent to ↑, ↓ for ordinary spinors) is sparse.
In Ref. [39], we evolve the many-body state |Ψ(t) via the classical RK4, a textbook [64] approach to solving ordinary differential equations. The goal of this effort was to extend to N ≥ 10 earlier calculations [33] that were performed by diagonalizing H efficiently via Bethe ansatz while studying the behavior of instantaneously conserved quantities of the system.
Another computational approach to numerically timeevolving the many-body state |Ψ(t) , working still in a sparse-matrix representation, is a Lanczos propagation [65][66][67] for a time-dependent Hamiltonian [68]. The Lanczos algorithm [69], a dimensional reduction method which generates the basis vectors of an effective (Krylov) subspace by repeated application of the Hamiltonian, is widely used, particularly in nuclear structure physics [70]. In our application, the many-body state is forwardintegrated according to Eq. (4) by applying a timeevolution operator; While formally one should compute the time-evolution operator U in a Magnus expansion (see Appendix A for details), in practice we found the naive time-evolution operator exp(−iH(t)δt), was sufficient. The Lanczos algorithm aids the efficient calculation of U by exponentiating the Hamiltonian projected into a very small but effective subspace, generated by applying powers of H on the state |Ψ(t) -again, details can be found in Appendix A.
Using established methods presented in Ref. [39] for evolving a many-body state with the Hamiltonian H(t), one can verify for N = 2 to 16 that using Lanczos propagation to approximate the evolved state to order (δt) 5 with the appropriately chosen δt produces accurate results for the wave function even after evolving over many time steps. When compared with results from RK4, the value of each coefficient in the wave function, j|Ψ(t) (j = 0, . . . , 2 N − 1), was calculated with a discrepancy 10 −3 . Because these methods are in numerical agreement with one another, we do not separately show explicit results for the time evolution of the wave function in the case of the Lanczos method. Note that the truncations involved in these numerical methods can cause the normalization of the resulting wave function to change, so one may need to normalize the resulting wave function between time steps. With each of these methods, we prescribe the time step to scale in N as inversely with the scaling of the difference between the extremal eigenvalues in our Hamiltonian: δt ∼ 0.1[µ N 2 ( N 2 + 1) + ω |ω|] −1 , where µ is evaluated at the radius prior to taking this time step. Just as with the use of RK4, we implement Lanczos propagation in a sparse-matrix representation, permitting calculations of the evolved many-body wave function according to a time-dependent Hamiltonian for up to N = 16 on a personal computer. In implementing the sparse representation in our own programs, we make use of submodules from SPARSKIT [71], a Fortran90 library for performing operations with sparse matrices.
However, because the number of nonzero elements in the many-body Hamiltonian matrix grows as O(N 2 2 N ), memory limitations severely restrict the values of N that can be studied. The time required to calculate the time evolution using these methods also grows exponentially in N . Tensor networks provide a method that in principle could scale more favorably with N . We describe these methods in the next subsection and investigate how the resources needed to obtain accurate results using tensor networks scale with N in the next section.

C. Calculating matrix product state wave functions
This exponential growth in problem size motivates the use of tensor network methods; appropriately chosen tensor network representations allow for the complexity of the problem to scale instead much more slowly with system size. However, it is not clear a priori how the size of the tensor network representation needed to obtain accurate results scales with N . Determining whether the necessary size grows slower than exponentially with system size is a key goal of this paper. The bond dimension required to obtain accurate results may scale either exponentially or polynomially for our system; using the methods that follow, we will in particular investigate how methods based upon the time-dependent variational principle may scale with N in treating our problem.
Here, we briefly outline how to express the many-body wave function and operators acting upon this state in terms of MPSs, followed by a sketch of the TDVP algorithm we use, reserving greater detail for Appendix B. The mathematical language associated with this decomposition will be referred to as "matrix product" or "tensor network" formalism interchangeably for our purposes. After establishing the mathematical definitions for describing a MPS, we will briefly outline the computational process of time evolving a wave function efficiently in a closed quantum system in the matrix product formalism, using the TDVP in the tensor network formalism. To this end, we largely adopt the language used in Ref. [44], which first outlined the version of a TDVP algorithm that we use in this work, to describe MPS formalism and its use with the TDVP. For a more general review of time evolution methods utilizing MPS representations of wave functions, see, e.g., Ref. [72].

General MPS review
Let us begin by establishing the language needed to describe a MPS. For a system of N neutrinos where we bin the spectrum of ω such that all neutrinos have distinct frequencies, we may label their frequencies with index values 1, . . . , N ; in the context of tensor network formalism, we can refer to these frequencies interchangeably as "sites." (Put another way, while the MPS community frequently considers sites in reference physical locations along a lattice, we are instead considering sites in reference to definite momentum states for different neutrinos in a spectrum.) Then, given a wave function |Ψ decomposed in the flavor basis we may view Ψ α1,...,α N = ν α1 , . . . , ν α N |Ψ as a complexvalued tensor with N indices each spanning a twodimensional vector space. By N −1 iterations of Schmidt decomposition (see, e.g., Refs. [35,73,74]) starting from the leftmost indices, we may write this component as a product of site-dependent matrices ψ αj L (j): where for fixed ( vector. This general decomposition into a matrix product is also referred to as an example of a "tensor train," and specifically with our choice of direction in Schmidt decomposition is a "left-canonical" form for the tensor. Here, we call β j and D j the "bond indices" and "bond dimensions" of our tensor train. An exact representation of Ψ is obtained if we take D j = min{2 j , 2 N −j }. With these exact choices for maximal bond dimensions, our procedure requires computational resources that scale exponentially with N . To reduce the scaling of computational resources with N , we seek to allow the maximum values of D j used in our computations to grow minimally with N while maintaining a similar level of error as obtained in methods such as RK4. Note that the bond dimension can help us to assess the entanglement in our ensemble; in the case that there is exactly zero entanglement entropy at each site, we find that D j = 1 for each j permits an exact representation of the state, resulting in an independent dim-2 vector subspace for each bodyjust as we would write in the mean-field theory calculations of the ensemble state [i.e., ψ αj (j) = α j |ψ(ω j ) in this case]. Regardless of our particular choice of D j , we order the singular values of each tensor by size and keep the ≤ D largest values.

TDVP for a neutrino MPS
Numerous recent developments have been made in the community studying time evolution of MPS representations of spin systems, whose accuracy is controlled in part by the choice of maximum bond dimension D use for a given calculation. In particular, in connection with MPS density matrix renormalization group (DMRG) techniques, Refs. [44,75,76] developed a method of real-time evolution in analogy with the TDVP. This TDVP algorithm in particular readily permits calculations with a spin Hamiltonian involving nonlocal interactions, with an acceptable level of accuracy reproduced for the case of a power-law potential [44,76]. Viewing (14) as N coordinates parametrizing a MPS manifold M of the state |Ψ , the TDVP can be interpreted geometrically as a projection of the right-hand side of the Eq. (4) onto the tangent space of said manifold at a location given by Ψ, T Ψ M , resulting in the nonlinear differential equations where P TΨM is the projection operator onto the tangent space. Specifically, we may choose different projection operators such that we evolve only one or multiple tensors ψ(j) in Eq. (14) at once; their forms and the consequences of each choice are also presented in Appendix B.
The TDVP algorithm implemented with the choice of having n active sites being evolved at once is referred to as nTDVP, with n = 1, 2 generally found to be practical computationally. Furthermore, an augmentation called global subspace expansion (GSE) has more recently been made to the nTDVP algorithm [45], the combination of which we call GSE-TDVPn, where n is the number of active sites being evolved in each TDVP step. In this procedure, one includes additional singular values from global Krylov vectors, calculated as [1 − iδtH(t)] |Ψ(t) ( ∈ N), into the bonds of the MPS wave function obtained between time steps of TDVP (in analogy with DMRG to optimize for a mixture of lowest-lying energy eigenstates). It was found that this addition provided greater flexibility to the choice of appropriate time-step sizes used for time evolution in cases such as the one-axis twisting model. However, the problem of collective neutrino oscillations in principle requires not only nonlocal interactions, but also the inclusion of one-body kinetic terms in the planewave treatment of neutrinos in our toy model. As such, it is not immediately apparent that entanglement describing the evolved many-body state of our system is accurately captured by these recent methods without requiring exponential growth in the bond dimension used to forward-integrate the wave function.
Notably, each virtual bond within the tensor train need not have identical dimension D j (j = 1, . . . , N − 1). In the repeated singular value decomposition to obtain a MPS (e.g., outlined in Refs. [35,74]), the bonds closest to the ends of the ω spectrum would have dim ≤ 2 while bonds closest to the center have dimension 2 N/2 . In 1/2TDVP algorithms made available through the TeNPy library [77] as well as GSE-TDVP1/2 algorithms in the ITensor library [45,78], we can control the maximum cut-off dimension for all of the bonds, which we denote by D. However, with an initial wave function in which neutrinos are unentangled, carefully note that 1TDVP calculations prevent the bond dimensions in the initial MPS from rising at all as time evolves, and correspondingly entanglement entropy is negligible throughout the calculation. Consequently, a one-site effective Hamiltonian calculation with an entirely fixed bond dimension can only replicate the results of an exact many-body calculation by beginning a time evolution using 2TDVP or GSE-TDVP1/2 for long enough to let all bonds reach their maximum permitted dimensions before then switching to 1TDVP. Moreover, we find that the greatest flexibility in the choice of time-step size and bond dimension, which are determined by a procedure described in Sec. V, is afforded by GSE-TDVP2. Therefore, this algorithm will be the MPS time evolution method used throughout our results in that section.
While Ref. [35] finds that a bond dimension D that scales linearly with the number of neutrinos is adequate in studying a two-beam model of collective oscillations, we find that the scaling of D with N is more complex for our case. By calculating the magnitude of the deviations of the results of our calculations from satisfying Ehrenfest's theorem, as outlined in Sec. IV, we can evaluate how much this restriction of bond dimension impacts the precision of the TDVP evolution of the wave function.
Beside choosing a bond dimension cutoff D in our tensor network calculations, we must also take care in choosing a time-step size δt throughout the evolution of our many-body wave function. As we will demonstrate in our Results, the problem of determining an appropriate δt to accurately evolve our MPS wave function is not entirely straightforward; while smaller time steps may help to more accurately forward-integrate our evolution equations, there are errors from ignoring singular values in both deriving our equations and following each time step, implying that shrinking δt to be too small can introduce even greater errors, as has been described, e.g., in Ref. [72]. 3 (For a more detailed explanation, see Appendix B.) In general, one must select a way to assess the error of a method without already knowing the exact solution to the problem; we shall propose a method for our problem in the following section. However, for initial conditions where entanglement is limited (i.e., |ν e ⊗N ), one can determine an appropriate δt in tensor network calculations by evolving |Ψ(t 0 ) via 2TDVP to t = t 0 + ω −1 0 using decreasing step sizes. Comparing the evolved wave function |Ψ n ≡ |Ψ(t 0 + 2 n δt n ) obtained with each step size δt n ≡ 2 −n ω −1 0 , we choose δt = δt n for the smallest natural number n such that 1 − | Ψ n+1 |Ψ n | 2 is less than some chosen tolerance value. Comparison with results obtained using sparse matrix computations in the complete Hilbert space for N ≤ 16 suggests that a tolerance of 10 −4 is appropriate for finding a practical δt in TDVP calculations of our system. For more general initial conditions and when using GSE in addition to TDVP, we will see below that determining an appropriate choice of δt will require greater attention. In fact, by carefully checking how error accumulates with differing δt, we find that certain larger values are often preferred in accurately evolving the MPS wave function than the upper bound prescribed for RK4 or Lanczos propagation by the argument outlined in Sec. III B.

IV. CONSISTENCY CHECKS FOR NUMERICAL TREATMENTS
To compare tensor network methods to the other methods described above, we need to consider both how the resources needed to calculate the evolution scale with system size and whether the accuracy of the solution that is obtained is adequate. Because the memory used and the computation time required grow polynomially with D [44] yet the maximum physical choice of D can grow exponentially in N , it is important to characterize the accuracy of the calculation when D is modest. For small N the results of the tensor network method can be compared to the numerically exact results obtained by the other methods, but for larger N it is useful to have another method to assess the accuracy of the results. In this section we discuss consistency checks that follow from conservation laws that can be used to assess whether the results yielded by the tensor network method with a given bond dimension are accurate.
It is known that the many-body neutrino Hamiltonian in the single-angle approximation has a number of commuting invariant operators. One such operator is J z = ω J z ω , i.e., the z-component of the total neutrino isospin in the mass basis. Another set of invariants, given by [26] and used further by [29][30][31] is These invariants can be used as consistency checks in numerical calculations. We do so by using Ehrenfest's theorem, which states that the time evolution of the expectation value of an operator A is given by 4 4 Using the chain rule of differentiation, one can write for a wave function |Ψ that satisfies the Schrödinger equation. As an example, taking A = J z in the above equation gives the simple result d J z /dt = 0, since J z has no explicit time dependence. Alternatively, taking A = h ω , one could, for instance, construct the norm to quantify how well |Ψ approximately solves the Schrödinger equation-if |Ψ solves the equation exactly, then the norm must vanish (since [h ω , H] = 0). Note that this norm is evaluating an overall uncertainty of sorts, if we assume the uncertainty for each h ω 's constraint is uncorrelated to that of the rest.
One may attempt to further simplify matters by inserting the form of h ω from Eq. (16) into Eq. (18). Doing so, one obtains Using the chain rule on the left-hand side leads to a cancellation, leaving us with Since P z (ω) = 2 J z ω , one may define the Ehrenfest error measure as where Ehr ω = 0 for each ω if |Ψ(t) is the exact evolved state.
Note that applying the Ehrenfest theorem again to the latter term of the above equation, and using Eq.(16) and the fact that [h ω , H] = 0, the above relation simplifies to which is simply the Ehrenfest theorem applied to J z ω . This condition could be used for a consistency check, to ensure that the many-body wave function obtained using any numerical approximation does indeed satisfy the Schrödinger equation to an acceptable level of accuracy. In particular, we propose the use of the maximum value of as a measure to assess how accurately a particular timeevolution algorithm is calculating our many-body state.
In the case of MPS calculations of real-time evolution, we point out that there may be a handful of initial time steps needed to transition an initial product state (where D = 1 for all bonds) to an entangled state where virtual bonds have as many nontrivial singular values as the max bond dimension that we set for a given calculation. As such, we find that the algorithm requires a few time steps to "warm up" and accumulate enough singular values to well-approximate our desired wave function, and so Ehr may vary in a less well-behaved fashion until the state is evolved to t ∼ t 0 + ω −1 0 ; as a consequence, we choose the domain for evaluating max Ehr ω [Ψ(t)] to be t > t 0 +ω −1 0 .

V. RESULTS
We wish to assess the resources needed to implement TDVP methods for MPS calculations of the dynamics of our model. Our first step is to address how to determine the minimum bond dimension needed to accurately time evolve a many-body state describing a dense neutrino gas. We use the max Ehr quantity discussed in the previous section to do so. We find that max Ehr may tend to be larger for initial conditions that result in multiple spectral splits, even for methods calculating the entire wave function such as RK4. However, this error quantity turns out to be independent of N , according to calculations for N < 20; as such, we propose the use of this quantity to determine when a MPS calculation is being carried out with insufficient precision. Using max Ehr provides a method to assess the accuracy at large values of N for which comparison to the results of other numerical methods is unavailable.
Throughout this section we will consider two kinds of initial conditions for our many-body state: one where all neutrinos begin in the electron flavor state, |ν e ⊗N , and one where the half of neutrinos with the lowest ω values start in the electron flavor state and the rest start in the x flavor state, |ν e ⊗N/2 |ν x ⊗N/2 . As pointed out in, e.g., Refs. [29,39], these conditions will result in one and two spectral splits, respectively. By taking these two cases into consideration, we may observe how our MPS methods handle evolved states with differing numbers of spectral splits.
We begin with a demonstration for the use of max Ehr in the relatively simple case of N = 4 calculations using GSE-TDVP2 [45]. Knowing the results of RK4 calculations to be very precise, we may compare the values of probability P ν1 (ω N ) defined in Eq. (7) after evolving to µ(t) ω 0 and the max Ehr for the wave function Ψ(t) [over all times t > t 0 + ω −1 0 ] with D ≥ 2, shown in Fig. 2. In the cases of each initial condition, we find that there is eventually a growth in error as one tries to decrease the time-step size δt to be too small. Even for the maximum physical bond dimension D = 2 N/2 = 4, results converge to values of lim t→∞ P ν1 (ω N ) that differ from those obtained using RK4. Recall that there are not only finite time-step errors related to forwardintegrating with the effective one/two-body Schrödinger equations for each tensor in train, but also truncation and projection errors as outlined in Appendix B. Since GSE introduces new singular values, which may not be physical if the number of values exceeds 2 N/2 , there will be a truncation error after each time step, independent of the size of the time step. Therefore, there is a truncation error that could grow at least linearly with the total number of time steps. (As a consequence there can be a nontrivial truncation error when using GSE even for D = 2 N/2 .) Additionally, comparing the two initial conditions, we observe that these max Ehr values are of a certain order of magnitude for certain choices of D and dt that match with max Ehr in RK4 results, suggesting that this quantity may help diagnose if our choices of parameters for evolution of our MPS can reasonably approximate the exact many-body wave function. For the mixed-flavor initial condition, we find that there is a region of δt values for each initial condition where max Ehr is small for at least certain values of D, while for other D values always result in max Ehr that is orders of magnitude larger. For the single-flavor initial condition, there is less sensitivity overall to choice of D, but there is a shared trend of converging to the wrong results after δt 0.01, where max Ehr quickly trends upward.
Taking this understanding of max Ehr, we can more easily approach calculations with larger N . For example, we can consider the case of N = 12, testing our method with the mixed-flavor initial condition described earlier in this section. We find that a time-step size δt = 0.01ω −1 0 is appropriate for varied choices; in particular, we depict in Fig. 3 the time evolution of P ν1 , S, and Ehr for each neutrino mode with the choices of bond dimension cutoff D = 45, 50, and 64. We find that the differences in P ν1 and S for D = 64 and 50 are relatively modest, corresponding with a difference in max Ehr values that is slightly less than an order of magnitude. In comparison, there is a larger discrepancy between D = 50 and 45, both with max Ehr as well as S and P ν1 . Specifically, one can see here that the lowest values of S areperhaps counterintuitively-overestimated by the use of too small of a cutoff D. (This observation is reflected also in 2TDVP calculations without the addition of GSE.) Consequently, the values of P ν1 permitted by Eq. (8) are thus more tightly bound, resulting in another observable difference, whereby probabilities P ν1 cannot approach 0 and 1. This relationship is isolated in Fig. 4, where we show, as an example, the evolution in the discrepancies for P ν1 (ω 1 ) as a function of time for D = 50, 45, and 30 from the maximum D = 64. Correspondingly, there is a yet larger growth in max Ehr between D = 50 and 45, suggesting a tipping point in one's choice of decreasing D where our predictions become progressively less accurate.
For these relatively small values of N , our system does not yet reflect a benefit in terms of complexity in our MPS calculations with growing N . However, we may look to cases of yet larger N in order to check that the D needed to obtain results at a desired level of accuracy (according to e.g., max Ehr) does not grow exponentially in N . In Fig. 5, we again consider the predicted time evolution of entanglement entropy for each ω, starting from the mixed-flavor initial condition but with N = 18; here, the largest physical bond dimension would be D = 512, which we compare against the choices of cutoff D = 300, 200, and 150. As in the case of N = 12, we find that lowering D too far can result in overestimates of the lowest values of S(ω) at a given time. We depict an example of this effect in Fig. 6, where values of P ν1 (ω 1 )| D deviate further from P ν1 (ω 1 )| D=512 throughout the evolution as D is decreased. Furthermore, in Fig. 7, we explore how these overestimates impact the prediction of the spectral split for this system; as in smaller N calculations, we find that the location and width of the split are unaffected, while the permitted range of values for probability P ν1 are more limited. However, in contrast with smaller N calculations, we find that one can reasonably approximate our system using values of D 300. In the same vein, max Ehr values are ∼ 10 −3 for D = 150 and 200 and ∼ 10 −4 for D = 300 and 512. This result suggests that a shrinking fraction of the complete set of singular values of the MPS are required by the GSE-TDVP2 algorithm to obtain accurate results as we increase N .
Carrying out an analogous comparison of time evolution results for varied N with these initial conditions, we summarize the temporal complexity of using GSE-TDVP2 with the smallest necessary bond dimensions and largest time-step sizes to reproduce the desired level of accuracy in max Ehr when compared to the methods of , we find that the final data for the wave function deviate from the RK4 results and converge slowly to an incorrect value for Pν 1 (ωN ). Correspondingly, we find that max Ehr is relatively well-behaved, with values O(10 −4 ) until δt 10 −2 , at which point we find max Ehr to increase by orders of magnitude, for all choices of D. This increase in error for these calculations with a decrease in the choice of δt is indicative of a growth in total truncation error from the end of each step in the GSE-TDVP2 algorithm. In the right case, we find that only the use of D = 4 produces results from our MPS method that agree closely with those of RK4. In kind, we find that max Ehr is consistently orders of magnitude smaller for D = 4 than for D = 3, 4. Though obscured by the greater discrepancies for D < 4 for the latter initial condition, there is a larger discrepancy for step sizes 10 −3 and 3 × 10 −3 correlating with the sudden increases in max Ehr.
RK4 and Lanczos propagation. We show these results as well as extrapolations from these data in Fig. 1. Extrapolated fit functions for RK4 and Lanczos propagation are based on O(N 4 2 N ) complexity discussed in Sec. III B, while there is an empirical trend line e −aN 2 +bN +c with a, b, c > 0 to depict one possible extrapolation of the data obtained from MPS methods for the initial condition resulting in one spectral split, as its well-behavedness permitted a greater number of calculations to display. 5 In case of GSE-TDVP2, for an initial state |ν e ⊗N , the computation time scales more favorably. Whereas, for a mixed initial state |ν e ⊗N/2 ⊗|ν x ⊗N/2 , scaling of computation time is less easily controlled, though calculations carried out on multiple cores (not shown in figure) suggest that choosing the smallest D to reproduce max Ehr values comparable to those of RK4 and Lanczos results t t t t t t t t t t t t FIG. 3. Time evolution from µ(t) = 5ω0 to µ ω0 of entanglement entropy (top row), probability Pν 1 (middle row), and Ehrenfest error "Ehr" defined in Eq. (23) (bottom row), for each neutrino mode in an ensemble of N = 12 with an initial spectrum of |νe ⊗6 |νx ⊗6 , evolved using a time step of δt = 0.01 and a maximum bond dimension of D = 45 (left column), 50 (middle column), or 64 (right column). Note that D = 64 = 2 12/2 is the largest physical choice of maximum bond dimension for a system of 12 neutrinos in two flavors, implying there is no error due to projection of the evolution equations; additionally, for D > 32, the only bond for which singular values are being ignored are on the central virtual bond (i.e., between sites 6 and 7). Consequently, we find that Ehr is most well-behaved for this case, while Ehr can reach values that are an order of magnitude larger, particularly at early times in the evolution. Correspondingly, we can carefully observe that as we decrease D, there is an overestimation of S(ω), especially for neutrino modes with the lowest values of S. As per Eqs. (8) and (9), Pν 1 is therefore bounded only to slightly smaller values. Interestingly, there is a greater discrepancy in results between D = 45 and 50 than that between D = 50 and 64, implying that our projection error while evolving this system grows dramatically as D decreases beyond the point of D ∼ 50.
in a shrinking fraction D/2 N/2 as N increases implying better scaling than RK4 and Lanczos scaling. Note that the data are obtained considering the smallest maximum bond dimension and lowest number of time steps required for a desired accuracy up to N 20. These parameters are expected to increase with N . Consequently, the scaling might differ from the ones presented here. Though bond dimension appears to grow more than linearly in N for our system, the growth could still be polynomial. Consequently, we see that the growth in complexity with N for our MPS methods appears optimistic, requiring a shrinking fraction D/2 N/2 as N grows, in contrast with the manifestly exponential scaling of RK4 and Lanczos propagation.

VI. CONCLUSIONS
This paper addresses the computational challenge of investigating the nature of collective neutrino oscillations in a dense neutrino gas. Employing the growing body of work in the MPS community to reformulate the problem in terms of a chain of effective two-body Schrödinger equations, we have investigated how TDVP methods (in- with respect to the maximum bond dimension D = 64. As in Fig. 3, we consider a system of N = 12 evolving from the state |νe ⊗6 |νx ⊗6 at µ(t) = 5ω0. As one excludes more virtual bond singular values by decreasing the value of D, again progressively greater discrepancies in the predictions of probabilities Pν 1 are found, due to progressively greater overestimates of the neutrino mode's entanglement entropy.
cluding GSE) may help in this line of inquiry. Furthermore, we have defined measures for error for many-body calculations from the instantaneously conserved charges of our Hamiltonian, to assess how well a calculated manybody wave function reflects the solution to our evolution equations. Where there is limited entanglement in our system (i.e., especially in conditions resulting in at most one spectral split), there are some use cases for GSE-TDVP2 in which MPS methods appear to scale much more favorably with N than other numerical methods such as Runge-Kutta (RK4). However, existing TDVP methods scale much less favorably when we consider initial conditions that require D ∼ 2 N/2 , in which case the temporal and spatial complexity will scale exponentially with N anyhow-although this growth will be slower than in the cases of the other algorithms.
The methods of GSE-TDVPn are new, and it remains to be seen whether further augmentations or improvements to existing methods can be made to permit yet larger N calculations (∼ 100) in a reasonable time frame. Conversely, one may point to several initial conditions in our model as systems in which quantum entanglement grows in such a way that even MPS representations may not efficiently treat the problem. In such cases there may still be an exciting opportunity for digital or analog quantum simulation hardware to assist in studying our Hamiltonian [38]. As we ramp up to larger N calculations with techniques such as those we have discussed, we hope to learn more about the scaling behavior of correlations t ω 7 ω 8 ω 9 ω 10 t ω 7 ω 8 ω 9 ω 10 ω 11 t ω 7 ω 8 ω 9 ω 10 ω 11 ω 12 t ω 7 ω 8 ω 9 ω 10 ω 11 ω 12 t ω 7 ω 8 ω 9 ω 10 ω 11 ω 12 t ω 7 ω 8 ω 9 ω 10 ω 11 ω 12 t ω 7 ω 8 ω 9 ω 10 ω 11 ω 12 t ω 7 ω 8 ω 9 ω 10 ω 11 ω 12 t ω 7 ω 8 ω 9 ω 10 ω 11 ω 12 FIG. 5. Entanglement entropy of each neutrino mode ω as a function of time for a system of N = 18, evolving an initial state |νe ⊗9 |νx ⊗9 from µ(t) = 5ω0, using δt = 0.01. We show calculations using maximum bond dimension values of D = 150, 200, 300, and 512, respectively from top to bottom. In particular, one can notice that the neutrino modes with the lowest entanglement entropy values throughout most of the time evolution are also those whose values are most greatly overestimated as we decrease D.
in collective neutrino oscillations with larger systems.  Again, neutrinos far from the spectral split frequencies experience greater overestimation of their entanglement entropy as the bond dimension cutoff D is decreased; consequently, the discrepancy in probability Pν 1 , as measured by χ 2 , grows.

Appendix A: Details of Many-body Wave Function Calculations
In this appendix, we elaborate on methods used to obtain the time-evolved wave function from Eq. (4) using the full 2 N -dim Hilbert space. In particular, we provide greater depth on how to perform further explanation regarding the use of Lanczos propagation.
In our application, the many-body state is forwardintegrated according to Eq. (4) by applying a timeevolution operator: |Ψ(t + δt) = U (t + δt; t) |Ψ(t) . Formally, one treats the time-evolution operator U by a Magnus expansion: Final (t = 1000ω −1 0 ) mass-1 probability spectra calculated via GSE-TDVP2, evolving the initial state |νe ⊗9 |νx ⊗9 from µ(t0) = 5ω0, with δt = 0.01. As in Fig. 5, we compare the results using bond dimensions D = 2 18/2 = 512, 300, 200, and 150. Interestingly, the locations as well as the widths of the spectral splits (ω/ω0 ∼ 3, 15) are unchanged by the reduction in bond dimension. However, we find that the range of permitted values in probability Pν 1 is decreased due to the overestimation in entanglement entropy observed in Fig. 5; therefore, the probabilities for modes furthest for the spectral splits are restricted to values closer to 1/2.
At the core of the Lanczos propagation method is the use of a low-dimensional effective basis in which one can apply an approximate time-evolution operator given by Eqs. (A1)-(A4). Because we found no noticeable advantage in going beyond the first term in the Magnus expansion, this task means computing the naive evolution operator u = exp(−iH(t)δt). Rather than exponentiating H in the full basis, which would be prohibitive for large dimensions, one constructs iteratively a k-dimensional subspace, the Krylov subspace, by orthogonalizing the set of vectors {[H(t)] |Ψ(t) } for = 0, . . . , k − 1; the representation of H(t) in this subspace is provided automatically as part of the algorithm. While in many applications the Lanczos algorithm is used to approximate extremal eigenpairs of a Hamiltonian [69,70], here we use it to construct the time evolution operator by approximate spectral decomposition.
For instance, let K be the 2 N × k matrix, constructed by the orthonormalized basis vectors of the Krylov subspace, that maps from the full Hilbert space to the Krylov subspace, in which the Hamiltonian is tridiagonal, and let V be a k × k real orthogonal matrix that diagonalizes H(t) K = K H(t)K T , the approximate Hamiltonian at time t projected into the Krylov subspace. As = V KH(t)K T V T is diagonal, it is trivial to exponentiate: u = exp(−i δt). Then one transforms back from the (approximate) eigenbasis to the Krylov basis and then finally to the original space, which is the time evolution or Lanczos propagation by one time step [68]. Notably, in the Lanczos propagation method, when taking the same time steps as those from our RK4 procedure, we find that the Lanczos algorithm needs very few iterations (i.e., typically 3-4) in order to arrive at a convergent result for the time-evolved wave function |Ψ(t + δt) as the number of iterations of the algorithm is increased. Additionally, we find negligible difference in the results of the evolution whether or not the second term of the Magnus expansion Ω 2 is included. Moreover, we find an agreement between this Lanczos propagation and RK4, whereby expansion beyond order (δt) 4 is unnecessary to replicate results obtained exactly with the Bethe ansatz method. 6 In the case of a constant Hamiltonian [i.e., constant µ(t)], the Lanczos algorithm would greatly simplify the complexity of the time evolution that results from applying the many-body Hamiltonian to an initial state |Ψ(t) , as the bulk of the time steps could be taken via calculations within the scalably small Krylov subspace. However, the case of a time-dependent Hamiltonian does not share this benefit, as the eigenbasis is evolving as well as the wave function, implying that we must change our basis out of the Krylov subspace following each time step. As a consequence, the time-dependent Lanczos propagation still suffers from the same difficulty as did RK4, in which a Hamiltonian that grows exponentially in N must be applied to a wave function repeatedly to time evolve our ensemble.

Appendix B: Tensor Network Formalism
In this Appendix, we provide greater detail into the MPS treatment of time evolution for long-range interacting systems such as that described by our Hamiltonian in this paper. First, we will elaborate on how we obtain different forms of MPS descriptions for our many-body wave function. Once this procedure has been outlined in greater detail, we will expound upon how the TDVP algorithm can be performed.
When we perform the truncation in the singular values with a bond dimension cutoff D to a state such as that in Eq. (14), where each bond dimension D j is replaced with min{D j , D}, we must normalize our state by imposing on the remaining entries of the tensor train the following constraint: 6 In close analogy with the Lanczos method, we could also propose the Fer expansion (see, e.g., [79]) as another numerical technique that approximates Eq. (12). Due to the extreme similarity at order (δt) 4 to RK4, we find that computation times are very similar between the two methods. Consequently, results produced from this method are not shown.
which we refer to as "left-normalization." This constraint also fixes a "gauge," the transformation of which leaves the MPS form unchanged (i.e., under insertions of G j G −1 j between each bond 1 ≤ j < N , where G is a D j × D j invertible matrix). For Eq. (B1), we define D N ≡ 1. We can repeat this same process of Schmidt decomposition while instead starting from the rightmost indices, and we replace the tensor symbols ψ L → ψ R and the bond indices β j →β j+1 and dimensions D j →D j+1 to denote this change in method, yielding a so-called "right-canonical" form. After performing the same truncation in the bond dimensions, we impose the "rightnormalization": where we defineD 1 ≡ 1. Going forward, it will also be useful to define D 0 ,D N +1 ≡ 1 to include the cases of j = 1, N automatically. We can then use these two decompositions to write left and right blocks of the MPS wave function: with which we construct the "mixed-canonical" form 7 where for the center site j we have the D j−1 ×D j+1 matrix ψ αj C (j) = ψ αj L (j)C(j) = C(j − 1)ψ αj R (j) with a D j ×D j+1 matrix C(j) containing the singular values for the "virtual bond" between sites j and j + 1. Moreover, using these definitions of left and right blocks, we can define orthonormal projection operators The particular MPS form in Eq. (B5) is immediately useful in the 1TDVP. As described earlier in Sec. III C 2, the 1TDVP algorithm involves approximating Eq. (4) by where the projection operator onto the tangent space, P TΨM , is given by 8 in the 1TDVP method; here, the jth site is the one "active" site in a given step, evolving exactly according to where H(j) is an effective one-site Hamiltonian at site j obtained using the projection operators described above.
Notably, this equation does not permit changes in bond dimension between sites and therefore limits the growth of entanglement in the system as well if the initial state is, for example, a simple product state (i.e., D = 1). In order to observe growth of entanglement as the many-body state evolves, we require a generalization of our earlier tensor train decompositions into left, right, and center blocks where we permit the center block to include multiple sites j, . . . , k: ψ C (j : k) such that we can write the entire state as |Ψ = αj ,...,α k =e,x × |Φ L,βj−1 (1 : j − 1) |ν αj · · · ν α k |Φ R,β k+1 (k + 1 : N ) . (B11) We provide diagrammatic forms for presenting a two-site center tensor as well as left-and right-normalized onesite tensors that can be chained together by contraction over virtual bonds to obtain a complete MPS of a wave function, using a style in keeping with the diagrammatic conventions presented in Ref. [44]. For example, one can depict a MPS with a two-site center by Fig. 8(a). We are then prepared to define a tangent space projector for the case of two active sites: We can then define an effective multisite Hamiltonian by applying projection operators such as the first series of terms above to H(t) from the left [e.g., depicted in Fig. 9(b) for a two-site center]. Replacing the choice of P TΨM with the two-site projection operator defined here, we obtain the two-site TDVP (2TDVP) equations where [H(j : j + 1)ψ C (j : j + 1)] αj αj+1 βj−1βj+2 × |Φ L,βj−1 (1 : j − 1) |ν αj ν αj+1 |Φ R,βj+2 (j + 2 : N ) We visually summarize the steps of 2TDVP in Fig. 9, again using a style in keeping with that of Ref. [44] for 1TDVP. Forward-integrating the many-body state |Ψ(t) according to the Schrödinger equation with the above form for the right-hand side is the 2TDVP algorithm, whose full list of instructions is given in Ref. [44]. More compactly, in order to carry out a step of 2TDVP, one forward-integrates the effective two-site evolution equations i d dt ψ C (j : j + 1) = H(j : j + 1)ψ C (j : j + 1), and subtracts from the resulting wave function the MPS obtained by backward-integrating the effective one-site FIG. 8. The tensor elements serving as the basic building blocks of tensor trains used to represent states and operators as MPSs and MPOs, respectively. In (a) we depict a tensor for a pair of center sites j and j + 1 where 1 ≤ j < N . Symbols β and β denote the indices for wide legs by which the factors of the train are connected via tensor multiplication, while symbols α denote the indices for thin legs by which the factors of the train connect to particular basis kets or bras for the corresponding sites. In (b) and (c) we depict a left-and right-normalized tensor for a site j, respectively. In all subfigures, we use β and β to denote internal indices for bonds to left-and right-normalized tensors, respectively, while we use the convention that bond indices β0, βN ,β1,βN+1 ≡ 1 are entirely, as there is no connection to a further site for the case of an open boundary condition. In keeping with the style of Ref. [44], we also use equilateral triangles pointing rightward (leftward) to symbolize a left-(right-)normalized tensor and rectangles to symbolize center site tensors, written in symbolic form by Eqs. (B1) and (B2). A contraction over a given virtual bond index βj orβj involves a sum over index values 1, . . . , Dj or 1, . . . ,Dj, respectively, while a contraction over an external flavor/mass index αj involves a sum over e, x or masses 1, 2. Additionally, the Tensor Network Python (TeNPy) library [77] provides several functions to help set up a program that evolves a many-body state in a MPS representation via a time-dependent Schrödinger equation via TDVP. As each time step is performed in the MPS formalism via an application of the Lanczos algorithm, this procedure can be thought of as a tensor network analog of the Lanczos evolution performed with a complete many-body state in sparse matrix representation. The finite time-step error of this algorithm is of order (δt) 3 , though the use of two active sites in our tangent space projections necessitates a truncation during singular value decomposition (SVD) to reduce the timeevolved ψ C (j : j + 1) → ψ(j)ψ(j + 1) that introduces error of size constant with respect to the choice of timestep size. Also, in contrast with the unitary evolution of 1TDVP, normalization of the wave function is no longer automatically preserved with 2TDVP if a truncation is performed at the end of a time step; in this case, one may need to divide the MPS by its norm between steps. As a consequence, it is (perhaps counterintuitively) desirable for the sake of precision to keep the time step from being taken as very small if one is to use 2TDVP over many time steps, as suggested by Ref. [72].
Before we conclude this section, let us discuss a more recent augmentation to the TDVP algorithms involving another avenue for growth in bond dimension between time steps. In particular, we follow the GSE method proposed by Ref. [45] to be used prior to each time step of TDVP; notably, this method does not depend in principle upon the number of active sites used during nTDVP (where n is the number of active sites), so this method introduces an algorithm for each choice of n: GSE-TDVPn.
The first of two steps of the GSE is to gather the Krylov subspace by which we will extend the MPS |Ψ(t) . We can obtain k − 1 states to extend the bond basis of Ψ in a numerically stable fashion by replacing [H(t)] with [1 − iδtH(t)] , as a first-order expansion of U (t + δt; t) ≈ e −iδtH(t) ≈ 1 − iδtH(t) for sufficiently small δt produces smaller changes to the norm of our vectors, yielding: Empirically, one finds in using GSE-TDVP that only a small value of k ∼ 5 and relatively little accuracy in obtaining the extra Krylov states > 0 are typically needed, implying a much larger truncation parameter for SVD can be utilized in this step-let us call it ε K -than that for representing our time-evolved wave function. Now, let us introduce the second step of GSE, in which we use the basis of the Krylov subspace above to extend our MPS for |Ψ . 9 Let us start with a given left-canonical form for |Ψ as in Eq. (14) and k − 1 additional basis MPSs |Ψ ( ) ( > 0) obtained in Eq. (B15). The general goal of this step is to incorporate singular values from the Krylov basis as we rewrite |Ψ as a MPS in rightcanonical form via N − 1 steps of SVD as mentioned earlier in this section. Starting from j = N and working our way to j = 1, we perform SVD at each site j as the orthogonality center of |Ψ : ψ C (j) → C(j − 1)ψ R (j) and define a projection operator P R (j) ≡ 1 − ψ R (j) † ψ R (j). Then, to ensure the additional bond bases of our final MPS are orthogonal to those of the original MPS, we project the tensors at the orthogonality center for each Ψ ( ) : ψ that any truncation parameter that we may use here ε M in mixing the Krylov states while neglecting small singular values needs neither to be the same as ε K from the earlier GSE step of obtaining the Krylov subspace nor to correspond to the truncation error of whatever D max we use during the TDVP steps.) Finally, we use the resulting right-orthonormal tensorψ R (j) to extend ψ R (j) via ψ R (j) ⊕ψ R (j).
Extending our time-evolved state in this fashion between TDVP time steps has proved useful in the case of a GSE-TDVP1 calculation of the real-time evolution of the one-axis twisting model [45]. More generally, it was proposed that GSE-TDVP1 allows the user to enlarge the tangent space of the MPS manifold before each time step, thus permitting growth in bond dimension even in models involving various kinds of non-neighbor interactions; even without using two active sites during TDVP, the extra Krylov states allow calculations to grow the bond dimension over time. This augmentation of the TDVP method may come with additional benefits, such as permitting smaller bounds on bond dimension under certain circumstances and use of relatively large time steps compared to other methods.