Exact solutions for the time-evolution of quantum spin systems under arbitrary waveforms using algebraic graph theory

A general approach is presented that offers exact analytical solutions for the time-evolution of quantum spin systems during parametric waveforms of arbitrary functions of time. The proposed method utilises the \emph{path-sum} method that relies on the algebraic and combinatorial properties of walks on graphs. A full mathematical treatment of the proposed formalism is presented, accompanied by an implementation in \textsc{Matlab}. Using computation of the spin dynamics of monopartite, bipartite, and tripartite quantum spin systems under chirped pulses as exemplar parametric waveforms, it is demonstrated that the proposed method consistently outperforms conventional numerical methods, including ODE integrators and piecewise-constant propagator approximations.


Introduction
Understanding the dynamics and control of quantum spin systems, indispensable in spectroscopy, sensing, quantum computing and information processing, is among the most challenging areas of research in current science. The diverse fields of applications include high-resolution magnetic resonance spectroscopy and imaging [1], terahertz technologies [2,3], and control of trapped ions [4], cold atoms [5] and NV-centers in diamond [6,7]. In many of these applications, sophisticated manipulations of quantum systems are achieved using pulses in the form of radio-frequency, microwave, or laser. A good understanding of the evolution of quantum spin systems during these events is crucial for enabling new methodologies.
A large class of these pulses are parametric, i.e. the waveform can be represented as a function dependent on certain time-dependent and time independent parameters. Numerical solutions of the spin dynamics during such pulses can be presented via solutions of ODEs e.g. using adaptive Runge-Kutta method [8], or methods relying on the approximation of matrix exponentials and Fokker-Planck formalism [9,10], but exact solutions of quantum systems under arbitrary parametric pulses are less explored. Although for frequency-swept pulses, an important class of parametric pulses, some approximated [11] and exact [12,13] analytical solutions for the time evolution of single spin-1 {2 particles have been presented in the literature. Other approaches that can offer similar solutions include integrable multi-state Landau-Zener Models [14,15] In this article we present very general exact analytical solutions for spin dynamics during parametric waveforms of arbitrary functions of time using the path-sum approach [16,17]. The approach relies on the algebraic and combinatorial properties of walks on graphs and here yields the complete description of the time-dependent propagator matrix in terms of the waveform and spin system parameters. The proposed approach is general and applicable to any user-defined pulses, arbitrarily constructed of parametric components e.g. chirps, hyperbolic secants, polynomials, Fourier series, etc. The method can be utilised to describe the time-evolution of quantum systems under the effect of parametric pulses. This paper is structured as follows: in section 2 the underlying theory for the equation of motion of quantum spin systems along with a general background on the path-sum approach is presented. In section 3 we present in detail the analytical solution for the propagation of monopartite systems driven by arbitrary time-dependent pulses in the path-sum formalism. In section 4 the method is extended to larger systems of spin-1 {2 particles and explicitly worked out in the case of bipartite and tripartite systems, demonstrating the applicability of the path-sum formalism to diverse Hamiltonian structures. Finally, in section 5 we present numerical results and assess the computational performances of Matlab codes implementing the path-sum formalism, using 6-parameter chirped pulses as exemplar waveforms, in the context of nuclear magnetic resonance (NMR) and the presence of scalar interactions.
and the Hamiltonian of the system under a pulse can be written as where Ω is the spin resonance offset, as commonly used in nuclear magnetic resonance (NMR), and β x and β y are real and imaginary components of the pulse The state of the system, ρ can be expressed as where gptq " rg 1 ptq, g 2 ptq, g 3 ptqs P R 3 is a time-dependent unit vector representing the position of the spin on a Bloch sphere. Considering 9 gptq "´iH ptqgptq, we can write¨9 Here time-dependencies of g 1 , g 2 , g 3 , β x , and β y are removed for simplicity. For a single spin-1 2 the basis set can also be written using shift operators where L`" L x`i L y and L´" L x´i L y . Using this basis set, the Hamiltonian under an arbitrary waveform can be written as where the bar¯indicates complex conjugation. Equation (5) using this basis set can be written as Switching between two representations of the Hamiltonian in Equations (6) and (9) can be achieved with a transformation matrix U as where and H SZ and H C are the Hamiltonians in the shift and Cartesian basis respectively.

Path-sum approach
The system eq. (9) is solved exactly with the method of path-sum, which we here present only briefly. We refer the reader to [16,17] for further details. Conceptually, the method relies on three mathematical pillars: i) The evolution operator solution of the quantum equation of motion is the resolvent of the Hamiltonian H if the usual multiplication between entries of H is replaced with another product, called Volterra composition [18].
ii) This resolvent is formally given by a sum over all walks on a graph G H whose adjacency matrix is the Hamiltonian operator and which thus encodes the discrete structure of the quantum state space.
iii) Sums of walks are given exactly by a continued fraction of finite depth and breadth involving a finite number of progressively simpler terms and stemming from an algebraic structure of sets of walks. Remarkably, each term of this continued fraction corresponds to a simple cycle or a simple path on the graph, i.e. to walks that do not visit any vertex more than once [19].
Combining these three observations, one can evaluate any time-ordered evolution operator exactly from a finite number of Volterra compositions and inverses between entries of the Hamiltonian as they appear successively along simple cycles and simple paths on the quantum state space. This conceptual framework has already been established in its full generality, so that when implementing path-sum concretely only writing the continued fraction needs to be done. This is because it changes according to the situation since the structure of the Hamiltonian dictates what the quantum state space graph G H looks like and thus which simple cycles/paths occur on it.
We begin by recasting eq. (9) in matrix form as d dt U "´iH ptqU ptq with H the Hamiltonian matrix of eq. (9) and U the corresponding evolution operator, which is the time-ordered exponential of´iH , As stated earlier, we can put this in resolvent form upon using the Volterra composition, a special case of the ‹-product [20]. remark 1. Fundamentally, the ‹-product is defined differently and on a much wider class of distributions of two variables, which include piecewise-smooth functions but also the Dirac delta and all its distributional derivatives. From there, it is also defined on matrices of distributions [20]. These technicalities are necessary to rigorously prove, among other things, that the Dirac delta distribution is the identity element with respect to the ‹-product and to perform ‹-inversions. Here we admit such results and proceed with the simplest definition of the ‹-product as Volterra did before the inception of distributions [18].
To define the resolvent form of U ptq, let f pt 1 , tq and gpt 1 , tq be two smooth functions of two variables. We define f ‹ g as For functions of less than two time variables, the variable must be treated as the left one, e.g. if hpt 1 q depends only on one variable, then ph ‹ gqpt 1 , tq " hpt 1 q ş t 1 t gpτ, tqdτ , pf ‹ hqpt 1 , tq " ş t 1 t f pt 1 , τ qhpτ qdτ . One can then show that [17] U ptq " where I ‹ :" 1 ‹ˆI " δpt 1´t qI is the identity matrix times the Dirac delta distribution 1 ‹ " δpt 1´t q and p.q ‹´1 designates a matrix inverse where all ordinary multiplications are taken instead to be Volterra compositions. To alleviate the notation and because it is related to Green's functions, the above resolvent p1 ‹ I´H q ‹´1 pt 1 , tq is denoted G pt 1 , tq. Then U ptq " ş t 0 G pτ, 0qdτ. Now path-sum guarantees not only that G pt 1 , tq is well defined, but also that any entry of G is given from ‹-products and ‹-inverses of entries´iH . This is illustrated in details in Section 3.

Solutions for general parametric waveforms
We consider the Hamiltonian of eq. (8), with the corresponding graph G H illustrated in Figure 1. The path-sum formulation of U 22 indicates that U 22 ptq "  .
Here all´i factors originates from the´i in´iH . The ‹-inverses in the above expression, e.g. in p1 ‹´p´i qΩq ‹´1 , stem from sums over all possible repetitions of a simple cycle c. Indeed, the observation 1`c`c.c`c.c.c`¨¨¨" ř n c n " p1´cq´1 is true at the formal level for cycles and walks [19], which implies its validity for actual weighted walks [17].
We can now evaluate the above path-sum explicitly. First, since Ω depends on strictly less than two time variables we have [16] 1 ‹´p´i qΩ˘‹´1pt 1 , tq " δpt 1´t q´iΩe´i Ωpt 1´t q .
Second, since β depends on a single time, for any function f pt 1 , tq we have With these observations, we obtain In these expressions, B t 1 ,t pΩ{2πq :" ş t 1 t e´i Ωτ βpτ qdτ is the Fourier transform with respect to τ of βpτ q`Θpτ´t 1 q´Θpτ´tq˘evaluated in Ω{2π, Θp.q being the Heaviside function with the convention that Θp0q " 1. is not a feature of the path-sum method itself, but rather of the peculiar form of the Hamiltonian H , in particular that Ω is time-independent.
The end result is remarkably simple and holds for all pulse shapes β, provided they are smooth in the mathematical sense. Numerical evaluation of the exact path-sum solution eq. (14) is straightforward using a discretized version of the ‹-product. Analytically speaking, the path-sum solution is best evaluated from its Neumann expansion, a non-perturbative series representation that is super-exponentially (and thus unconditionally) convergent and lends itself to analytical calculations. In order to facilitate the notation, we introduce Then, since the Neumann expansion of the path-sum solution is here displaying only the first two orders. Higher order terms are analytically accessible and can reach any desired accuracy. To illustrate this, consider the m th Neumann order approximation In this notation U p8q 22 ptq is the exact solution, which we emphasize that it can be evaluated numerically directly without resorting to Neumann series. Indeed, Volterra compositions map to ordinary matrix products when time is discretized [21] so that the ‹-inverse of eq. (14) becomes the ordinary inverse of a well-conditioned triangular matrix, directly yielding the numerical evaluation of U p8q 22 ptq. Nonetheless, the Neumann approximations U pmq 22 ptq are useful in that they provide analytical closed-form expressions for any desired accuracy.
Using path-sum we can evaluate the other entries of G similarly.
and G 33 ptq " Ğ G 11 ptq. Off-diagonal terms are given by sums over simple paths on the graph of Figure 1. Here we have for example , which leads to and G 32 ptq " Ğ G 12 ptq. Here U 22 and G 22 are given by Equations (14) and (16) respectively. Similarly we get and U 21 " ş t 0 G 21 pτ qdτ . By symmetry, U 23 ptq " Ğ U 21 ptq. We conclude with and, by symmetry, U 13 ptq " Ğ U 31 ptq. Just as for the diagonal terms, all of the path-sum results for the off-diagonal terms yield unconditionally convergent Neumann expansions, e.g.
. Solutions in SUp2 N q

Liouville-von Neumann equation
The solution of the Liouville-von Neumann equation in the Hilbert space can be written using a time-dependent unitary propagator as where U ptq is the solution of the differential equation In general, regardless of whether the Hamiltonian commutes with itself or not, the solution is eq. (12).

Solution for monopartite systems
The Hamiltonian of a monopartite quantum spin-1 {2 system with a resonance offset Ω can be written as Since the Hamiltonian is 2ˆ2 with no zero entry, the corresponding graph G H is the complete graph on 2 vertices and the path-sum formulation of the corresponding evolution operator is [16] ş t 0 G 11 pτ, 0qdτ and

System of M interacting spin-1 2 s
In an M -spin system, the general form of spin operator for the i th spin can be expressed as with a Pauli matrix in position i. The total Hamiltonian can be written as a sum of the Hamiltonian for spin resonance offset (H Ω ), pulse (H P ), and some scalar interaction, known as J-coupling in NMR (H J ) and

Solution for bipartite systems
The Hamiltonian for a bipartite quantum spin-1 {2 system with offsets Ω 1 and Ω 2 and a coupling constant J can be expressed as where For the sake of simplicity regarding the equations that give the evolution operator, and without loss of generality, let us change the energy gauge to with I the 4ˆ4 identity matrix. Of course, this Hamiltonian leads to the very same dynamics as of the original H of eq. (23). Now we rely on path-sum's scale invariance to express the solution of this spins system. Mathematically speaking, this 'scale invariance' refers to pathsum's continued validity for all partitions of the Hamiltonian matrix into (possibly non-contiguous) blocks. In other terms, there is a path-sum formulation of the evolution operator U ptq in terms of any chosen partition of the Hamiltonian into blocks of any size, including non-square ones. Physically, this means that path-sum is able to formulate the dynamics of a quantum system in terms of the isolated dynamics of any chosen collection of its subsystems. We may also use partitions exhibiting mappings from one Hamiltonian to another, so as to demonstrate the similarity in the time evolutions they effect.
To illustrate the direct relation between monopartite and bipartite dynamics, we partition the Hamiltonian H 1 ptq as follows where we defined u :" p1, 1q and is the 2ˆ2 submatrix of H ptq formed by entries H 1 ptq i,j with i, j " 2, 3. In these expressions, h 1 ii " h ii´p 1{2qπJ. We chose this partition so has to exhibit an exact mapping from the bipartite system to an effective monopartite case with non-Abelian energy H II for the degrees of freedom spanned by states p0, 1, 0, 0q and p0, 0, 1, 0q. This mapping is best seen on the graph G H 1 , shown in Figure 2, which encodes this partition of H 1 and is structurally similar to that of one the spin case shown in Figure 1. The partition in itself plays no special role in the analytical solution however, and the problem could be solved using the ordinary entries of H 1 directly. Now let U ptq be the evolution operator of the bipartite system-i.e. the solution of (20) with Hamiltonian H 1 -and let G " 9 U . Using the same partition for U as for H 1 , let U II be the 2ˆ2 submatrix of U formed by the entries U ptq i,j with i, j " 2, 3. Then U II ptq " ş t 0 G II pτ qdτ with G II the Green's function given by the path-sum expression where I ‹ " 1 ‹ˆI " δpt 1´t qI is the two by two identity matrix times 1 ‹ " δpt 1´t q and P " u T .u is the 2ˆ2 matrix full of 1. As predicted above, G II is similar to (14) for G 22 in the monopartite case. Note that the additional iH 1 II ptq term would also have been present in the one spin case had there been a non-zero H 22 diagonal term in the Hamiltonian. The path-sum for G II involves a matrix ‹-inverse, meaning that G II enters a non-Abelian linear Volterra integral equation of the second type, with unconditionally convergent analytical representation given by the matrix-valued Neumann series of ‹-powers of its kernel.
Furthermore, one can use path-sum's scale invariance to get scalar-valued expressions for any entry of G II rather than manipulating matrices. For example, entry pG II q 1,1 " G 2,2 is given by the path-sum where bpt 1 , tq " B Bt 1ˇBt 1 ,t`p Ω 1`Ω2 q{4π˘ˇˇ2. Vector G 1,II " pG 1,2 , G 1,3 q then fol-lows as Similarly we have The remaining entries are succinctly obtained from path-sum in matrix form aŝ e´i H II pt´τ q .ˆβ pτ qβpτ q βpτ qβpτ q˙.ˆU 1 pτ q U 1,4 pτ q U 4,1 pτ q U 4 pτ q˙d τ.

Solution for tripartite systems
The Hamiltonian for a tripartite quantum spin-1 {2 system with offsets Ω 1 , Ω 2 , and Ω 3 and coupling constants J 12 , J 13 , and J 23 can be expressed as This system presents no further difficulty than the monopartite and bipartite cases. We may exploit different partitions of the quantum state space to exhibit mappings either from the tripartite to the bipartite or monopartite cases, or to an altogether different situation such as a mathematically advantageous structure, all thanks to path-sum scale invariance. Here we choose to map the Hamiltonian to a path-graph, whose path-sums yield continued fractions with a single branch. Let |iy designate the canonical basis states, e.g. |1y " p1, 0, 0, 0, 0, 0, 0, 0q T . Let V I " spanp|1yq, V II " spanp|2y, |3y, |5yq, V III " spanp|4y, |6y, |7yq and V IV " spanp|8yq. At the scale formed by these vector spaces, the Hamiltonian takes the form where u " p1, 1, 1q and M "¨1 The corresponding graph is a path-graph on 4 vertices, illustrated on Figure 3.
To succinctly present the path-sum expression for the evolution operator let G " 9 U , define I ‹ :" Iˆ1 ‹ the 3-by-3 identity matrix times the Dirac Delta distribution 1 ‹ " δpt 1´t q and let

II III
here all inverses are to be understood as ‹-inverses, i.e. inverses with respect to the ‹-product, the above presentation being chosen so as to reveal clearly the continued fraction nature the path-sum expressions. Now we have access to all entries and block of entries of the Green's function G , with the path-sum expressions while the off-diagonal blocs are This provides the exact, fully analytical formulation of the evolution of the tripartite systems under a time-dependent driving. These expressions can all be expanded analytically into closed form approximations via truncated unconditionally convergent Neumann series. They can also be directly evaluated numerically when time is discretized using strategies presented in the following Section.

Numerical computations
A Matlab implementation of the path-sum approach, presented in this work, has been made freely available [22]. Here we present an assessment of its computational performances compared with the standard piecewise-constant propagator approximation, based on the identity and adaptive Runge-Kutta method, using ode45 in Matlab.

Exemplar waveform: chirped pulse
In this article we employ a general expression of βptq for a chirped pulse [8], which can be expressed using 6 parameters: amplitude (ω 1 ), bandwidth (∆F ), duration (τ p ), overall phase (φ 0 ), time offset (δ t ), and frequency offset (δ f ). The time envelope (amplitude profile) of a chirped pulse can be expressed using a super-Gaussian distribution where n is a smoothing factor. The frequency sweep function can be written as which correspond to the phase The maximum amplitude of a chirped pulse (ω 1,max ) can be calculated using three parameters as where Q 0 is the adiabaticity factor at time τp {2, where t P r0, τ p s and δ t " τp {2.
As the effective flip angle of a chirped pulse approaches 180˝asymptotically as the ω 1 increases, for most practical purposes a value of Q 0 is chosen to satisfy adiabatic condition, while avoiding excessive pulse amplitudes. Figure 4 shows comparisons between numerical and analytical results for single spin-1 {2 under a chirped pulse. Figure 5 shows the performance of several analytical approximations as in eq. (18), for a chirped pulse with realistic parameters, versus the true solution. Figure 6 shows the time-evolution of a spin-1 2 experiencing an adiabatic inversion on a Bloch sphere, obtained from the path-sum solution described in Section 3.1.

Numerical implementation of path-sums
Let I " r0, T s be the time interval over which the numerical solution is sought. Let tt i P Iu 0ďiďN´1 be the discrete times at which the solution is to be numerically evaluated. For simplicity, we take these time points to be equally  (14) and (16) (dashed red line), for a monopartite system with Ω " 2π 7kHz under a chirped pulse, described in Section 5.1, with ω 1 " 2πˆ1545 rad/s, ∆F " 30kHz, τp " 10ms, φ 0 " 0, δt " 5ms, δf " 0, and n " 30. Two curves show perfect overlap and hence are indistinguishable. Inset shows the normalised temporal amplitude profile 2|βptq|{ω 1 .
spaced by ∆t. For a smooth functionf pt 1 , tq over I 2 , we define the triangular matrix F with entries With discretized time, the Volterra composition ‹ of eq. (13) turns into an ordinary matrix product where G is the triangular matrix built from g. Most importantly, this observation extends to ‹-resolvents Therefore, in practical numerical computations with ∆t ! 1 the ordinary inverse of the triangular matrix I´∆t F approximates the ‹-resolvent of f . Furthermore, matrix I´∆t F is necessarily well conditioned since its diagonal entries 1´∆t F i,i can be made arbitrarily close to 1 by choosing ∆t small enough. Consequently, numerical evaluation of any path-sum requires only multiplying and inverting well-conditioned triangular matrices [21,24] We improve upon this strategy by noting that eq. (34) corresponds to using the rectangular rule of integration. Using the trapezoidal or averaged Simpson rule [25] instead leads to much more accurate results. Standard numerical analysis indicates that a code using trapezoidal quadrature on N time points and time step ∆t " T {N should have an accuracy scaling as Op∆t´2q and a computational cost of OpN 2 q. Similarly, a code relying on the average Simpson rule of integration should have an accuracy scaling as Op∆t 3 q for a computational cost of OpN 2 q. In contrast, the piecewise-constant propagator approximation of eq. (28) is expected to have accuracy Op∆tq and linear computational cost OpN q.
Seeking more flexibility in the calculations, the path-sum codes allow for subdivisions of the time interval I into N I smaller intervals. Over each of the subintervals, the evolution operator is evaluated from its path-sum formulation on N p discrete time points. The evolution operator at the end of the interval is passed as a seed to the next interval, over which it is once again evaluated from its path-sum. This offers an hybrid approach between piecewise-constant propagator approximation (PCPA) (N p " 1, N I " 1) and pure path-sum pN p " 1, N I " 1). Overall, this hybrid approach evaluates the solution at a total of N " N p N I time points. Tuning N p and N I independently allows the user to trade accuracy for speed and vice-versa. For example, keeping N p moderate while choosing N I " 1 allows for faster evaluations than N p " 1, N I " 1 and is thus better suited should the user require the solution at a very large number of time points. At the opposite, more accurate results will be obtained by making N p large while N I can be as low as 1. In practise, we found that in most-though not all-situations a pure path-sum approach (N p " 1, N I " 1) offers the best trade-off of speed versus accuracy. The codes designed to produce the entire discretized time-evolution of the density matrix taking the spin system and the waveform parameters as inputs. Figure 7 shows the output density matrix for a bipartite (4ˆ4) system undergoing an adiabatic inversion using a chirped pulse.

Accuracy evaluation
Let ρ M ptq be the density matrix as evaluated by a method M (e.g. pathsum or PCPA) and let ρ r be the reference density matrix as evaluated within machine precision (relative and absolute tolerances set to 10´1 3 for each entry) by the standard solver ode45. We evaluate the accuracy of ρ M ptq by evaluating  Figure 7: Example of an output produced by the path-sum code using Simpson quadrature in the case of a bipartite system with relative error E M " 3.28ˆ10´9 using eq. (35). Each element of the figure represents the real part of the density matrix elements ρ i,j ptq as a function of time t (in µs). Spin system parameters: Ω 1 " 2πˆ700 rad/s, Ω 2 " 2πˆ600 rad/s, J " 150 Hz; waveform parameters: τp " 1ms, ∆F " 50kHz, φ 0 " 0, ω 1 " 2πˆ6.31 krad/s, δt " 0.5ms, δ f " 0 and n " 20.
the deviation from 1 of its normalised Frobenius scalar product with ρ r which is the relative error on ρ M with respect to the reference solution. Here }A} F :" TrpA : Aq designates the Frobenius norm of matrix A. As constructed above, the relative error E M evaluates to 0 if ρ M ptq " ρptq at all times.
These tests show that for relative errors E M on the order of 10´3 to 10´7largely sufficient to simulate actual experiments-path-sum codes outperform the piecewise-constant propagator approximation (see in particular Figure 8) and even ode45. The latter observation is surprising given that: i) hear we are dealing with small systems for which ode45 is highly optimized; ii) ode45 is allowed to choose the number and position of its computation points while pathsum is constrained by our implementation to work on equally spaced points-a known disadvantage; iii) ode45's computation time degrades rapidly when asked to produce the solution at a large number of evaluation points.
Further improvements to the current path-sum implementation can be made by using optimized, instead of linearly spaced, time points or working with an orthogonal polynomial basis. Additionally, C or Python implementations can offer significant speed-up compared to the present Matlab one, as evidenced by path-sum based codes for computing Heun functions [21,26].

Numerical outlook for large systems
As the size of the Hamiltonian increases, we expect the difference of computation times to grow in favor of codes evaluating the path-sum solutions, in particular, for codes exploiting path-sum's scale invariance property [16] or using Lanczos path-sum [24,27,28]. Lanczos path-sum is a type of pre-conditioning procedure for path-sum that finds a tridiagonal time-dependent matrix T ptq whose line 1 column 1 entry of the time-ordered exponential is the same as   (3) representation driven by a chirped pulse with parameters given in the caption of Table 1.
Here path-sum codes solve the whole 3ˆ3 system of Bloch equations irrespective of the initial state, ψp0q. In contrast, the PCPA-based algorithm evolves ψptq from one given ψp0q, an easier task.  Table 3: Comparisons of computational performances in the bipartite cases with parameters given in the caption Figure 7. " v T .U ptq.w.

Method
While the above equality is exactly satisfied for a tridiagonal matrix T of the same size as U , an excellent approximation of v T .U ptq.w may be achieved us- ing a truncation of T that is much smaller than U . This procedure, standard in time-independent Lanczos methods, partially alleviates the problem posed by the exponential expansion of the quantum space with the number of spins. The ‹-Lanczos approach heavily exploits the Hamiltonian's sparsity and can effectively reduce the size of the matrices involved while maintaining numerical accuracy by truncating T ptq [24]. Once a suitable truncation of T is determined, path-sum is exploited to compute its time ordered exponential. Because T and its truncations are all tridiagonal the path-sums are finite scalar-valued ‹-continued fractions with a single branch similar in structure to those presented in this work.

Conclusion
In this work we used the path-sum formalism and showed the wide applicability of this approach to solve, both analytically and numerically, the spin dynamics of quantum spin systems driven by time dependent forces. Analytically speaking, the method offers closed form expressions in terms of ‹-resolvents, which may nonetheless themselves be transcendent mathematical functions. However, these ‹-resolvents are always available from their analytical unconditionally converging Neumann series expansions or numerically, to any desired accuracy. Using the time-discretization of Volterra compositions, numerical implementations of path-sums require only multiplying and inverting well conditioned triangular matrices. Via a diverse set of examples, we demonstrated that the resulting numerical path-sum integrator consistently outperforms piecewiseconstant propagator approximation and ODE integration methods, which fur-thermore do not provide any analytical insights in contrast with the path-sum approach.
The proposed method has the potential to find applications in a variety of areas where having access to exact solutions of the time evolution of quantum spin systems is beneficial, including geometric [29][30][31][32][33] and adiabatic optimal control [34][35][36][37] methods, in addition to optimal control of quantum spin systems [38].