Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR Pierre-Louis Giscard, Christian Bonhomme


I. INTRODUCTION
The unitary evolution operator U(t , t ) describing the time dynamics of quantum systems is defined as the unique solution of Schrödinger's equation with quantum Hamiltonian H, i.e., (h = 1) and such that U(t = t, t ) = Id is the identity matrix at all times. Evidently, this operator plays a crucial role at the heart of quantum mechanics, including for spin dynamics in nuclear magnetic resonance (NMR) [1][2][3]. As is typically the case in NMR, the Hamiltonian may be time dependent and might furthermore not commute with itself at various times, H(t )H(t ) − H(t )H(t ) = 0 for t = t. In this situation, the evolution operator no longer has a simple calculable form in terms of the Hamiltonian, e.g., it cannot even be evaluated via direct diagonalization of H. Rather, U is formally described by the action of a time-ordering operator on the Dyson series representation of the quantum evolution [4], a formulation which does not permit concrete calculations to be carried out. As a consequence, only approximate expressions of U(t , t ) are obtained and these are only accurate for short times. A major breakthrough in the description and understanding of solid-state NMR was the inception of average Hamiltonian theory [5]. This relies exclusively on the Magnus expansion [6] of U(t , t ). However, higher-order terms of the series remain highly cumbersome to write explicitly so that practically, only low orders of the expansion are usable. Most importantly, Magnus expansion suffers from severe and incurable divergences as already mentioned by Magnus [6] and Fel'dman [7]. In the more specific case of periodically time-dependent Hamiltonian, such as those encountered in magic angle spinning (MAS) experiments, it is well known that Magnus expansion suffers from a further two limitations, i.e., the stroboscopic detection of the NMR events, and the impossibility to take into account more than one characteristic period.
In the case of periodic Hamiltonian, Floquet theory dictates that the evolution operator takes on the form U(t, 0) = P(t ) exp(Ft ), with P(t ) a periodic time-dependent matrix and F a constant matrix, both of which are determined perturbatively when working analytically [8], or otherwise via numerical procedures [9,10]. Floquet formalism was first used by Shirley [11], who applied it to the case of a linearly polarized excitation in magnetic resonance and to give low orders analytical expressions for the Bloch-Siegert effect [12].
We also mention numerical methods: (i) Fer and Magnus-Floquet hybrids proposed recently as potential expansions for the evolution operator [13,14], (ii) Zassenhaus and Suzuki-Trotter propagator approximations [15][16][17]. The expansions presented above all suffer from various drawbacks including the divergence of the series at long times; the perturbative nature of the numerical or theoretical approach; the nonavoidable propagation of errors at long time; the failure to find exact solutions even for small, one spin- 1 2 , 2 × 2 Hamiltonians. See also Appendix A for more literature background.
In this contribution, the path-sum method is applied to NMR Hamiltonians to determine the corresponding evolution operators U(t , t ). The rigorous underpinnings of this approach were laid out in [18,19] within the general mathematical framework of systems of coupled linear differential equations with nonconstant coefficients. So far, no physical applications of these works have been presented. Consequently, they remained unnoticed outside of a specialized mathematical community, and their applicability to longstanding questions pertaining to quantum systems driven by time-dependent Hamiltonians went completely unrecognized. It thus appears important to introduce path sum to the physics community via illustrative examples bearing directly on currently open problems. Overall, it appears that this work is a different approach to the problem of simulating quantum dynamics induced by time-varying Hamiltonians since Magnus' 1954 seminal results.
Path sum is firmly established on three fundamentally novel concepts, insofar never applied within the quantum physics framework: (i) the representation of U(t , t ) as the inverse of an operator with respect to certain * product; (ii) a mapping between this inverse, and sums of weighted walks on a graph; and (iii) fundamental results on the algebraic structure of sets of walks which exactly transform any infinite sum of weighted walks on any graph into a single branched continued fraction of finite depth and breadth with finitely many terms. Taken together, these three results imply that, for finite-dimensional Hamiltonians, any entry or block of entries of U(t , t ) has an exact, unconditionally convergent analytical expression that always involves a finite number of terms. We emphasize that throughout this work, the time t is and remains a continuous variable, in particular path sum does not rely on time discretization. As a corollary, path sum yields a nonperturbative formulation of U(t , t ), as will be illustrated below with the Bloch-Siegert effect. Further properties of path sums ensure its scalability to multispin systems, for example, allowing it to recover the exact dynamics of an entire system from the separate, isolated, evolutions of any chosen collection of its subsystems. In its general form, path sum is best understood as a method to exactly and analytically solve systems of coupled linear differential equations with nonconstant coefficients. This paper is structured as follows. We first present the mathematical background of [18] culminating in the path-sum formulation of quantum dynamics. In a second part, we detail applications in connection with general quantum theory and then more specifically with NMR. The first one provides the general solution of Schrödinger's equation to all 2 × 2 timedependent Hamiltonians, a problem of current and central importance to quantum computing. As an example of application, we solve for the celebrated Bloch-Siegert dynamics of a linearly polarized radio-frequency (rf) excitation with no approximation at all. The validity of the path-sum analytical solution is demonstrated over the entire driving range, and physical interpretations for the various terms of the solution are provided.
We then show that path sum is invariant against scale transformations in the quantum state space, making it scalable to large quantum systems. Thanks to this, we consider N like spins coupled by the homonuclear dipolar coupling and spin diffusion under MAS. The effects of MAS frequency and chemical shift offsets are illustrated analytically on an organometallic molecule exhibiting 42 protons.

II. QUANTUM EVOLUTION AND WALKS ON GRAPHS
Quantum systems with discrete degrees of freedom such as spin systems obey a discrete analog to Feynman's path integrals. To illustrate this, define one history of a quantum system as a temporal succession of orthogonal quantum states h : |s 1 → |s 2 → |s 3 · · · , each transition |s i → |s i+1 happening at a specified time t i . Overall, the history h acquires a complex weight which is the product of the weights of all the transitions in the history. The weight of an individual transition |s i → |s i+1 is dictated by the Hamiltonian as s i+1 |H(t i )|s i .
A natural representation of such discrete histories is as walks on a graph. To see this, let G t be the graph such that each vertex v i corresponds to one member |s i of an orthonormal basis for the entire state space and give the directed edge v i → v j the time-dependent weight s j |H(t )|s i . In this picture, a system history as defined earlier is a walk on G t and H(t ) is the adjacency matrix of G t . Because the Hamiltonian is time dependent, the graph itself is dynamical [see Figs. 1(a) and 1(b) for an example]. Now, just as for Feynman's path integrals, the exact evolution of the system is obtained from the superposition of all its possible histories. Equivalently, every element s j |U(t )|s i of the evolution operator U(t ) is given by the sum over all walks from v i to v j on G t , including all possible jumping times for each transition between vertices. While individual walks are the discrete counterpart of Feynman diagrams, their algebraic structure is much better understood. Indeed, walks essentially behave as the natural integers [19], in particular they can be uniquely factored into products of prime walks: the simple cycles (C) and paths (P) of the graph which do not visit any vertex more than once. Since, by nesting simple cycles and paths into one another there is a unique way of reconstructing any walk, summing over all walks is achieved upon summing over all possible nestings of simple cycles and paths. For example, in a graph with a single simple cycle c 1 , all closed walks from a vertex α to itself are of the form c n 1 , i.e., c 1 repeated n times. Therefore, the sum of all such walks is formally n c n 1 = 1/(1 − c 1 ) (Fig. 2). In case another cycle c 2 is accessible to a walker while walking along c 1 , then the Step-by-step evaluation of G K 2 ,11 (t , t ) (dashed rectangle) showing the finite character of the continued fraction. The sum is performed on the simple cycles (C) of lengths 1 and 2 (these are indicated in red and other edges are indicated in dashed gray lines). At each step of the continued fraction, a vertex is removed (gray cross) and we work on subgraphs of less and less complexity. The calculation of entry U 21 from G K 2 ,11 is also illustrated; it includes a single term with two * products as the graph has a single simple path (P) from 2 to 1 (red arrow). (e) A pictorial representation of the "descending ladder principle." The calculation starts at the top of the ladder with each * inverse leading to the step below and ending in all cases on the ground. unconditionally convergent when calculated numerically. The same principles apply regardless of whether the Hamiltonian depends on time or not, in the former case, however, the theory relies on two-times functions f (t , t ) that multiply via a noncommutative convolutionlike product This means that for general time-dependent Hamiltonians the continued fraction formulation for U(t ) involves products and resolvent with respect to the * multiplication and that the order of traversal of the edges along the cycles must be respected. The * resolvents such as (1 * − f ) * −1 with 1 * ≡ δ(t − t ) the Dirac delta distribution are solutions to linear Volterra equations of the second kind. They concentrate most of the analytical complexity of the problem, rarely having a closed-form expression in terms of algebraic mathematical functions. Yet, such * resolvent can be always represented analytically by the superexponentially convergent Neumann expansion [18] 1 * + n>0 f * n .

A. General solution
Determining the dynamics of two-level systems driven by time-dependent Hamiltonians is still an open problem, which continues to be a very active area of research [8,[20][21][22][23][24][25][26][27][28][29]. This is because of the experimental relevance of such systems, their role as theoretical models, and the need to master the internal evolution of qubits undergoing quantum gates [27,30]. The most general two-level Hamiltonian is In this expression we only require that h ↓↑ (t ), h ↑↓ (t ), h ↑ (t ), and h ↓ (t ) be bounded functions of time over the interval [t, t ] of interest. So far, no analytical expression has been found for the corresponding evolution operator U(t ) defined as the unique solution of Eq. (1) with the Hamiltonian of Eq. (3). It is known that particular choices for H(t ) lead to evolution operators that involve higher transcendental functions [21,31]. Thus, the best possible result for the general U(t ) is that each of its entries be described as solving a defining equation, and that an analytical mean of generating this solution be presented. This is exactly what path sum achieves for all timedependent two-level systems. Following the process of Fig. 1 we get while the "usual" U(t ) is actually U(t, 0). The two-times Green's functions G ↑ := (1 * − K ↑ ) * −1 and G ↓ := (1 * − K ↓ ) * −1 solve linear Volterra equations of the second kind, e.g., for G ↑ , and similarly for G ↓ . The kernel K ↑ of the above equation is while kernel K ↓ entering G ↓ is obtained upon replacing up arrows by down arrows and vice versa in Eq. (6). Should a closed-form expression for the solution of the Volterra equation be out of reach, e.g., because it is a transcendental function as is typically the case [21], the solution is at least analytically available from its Neumann expansion; in the case Eq. (5) it is G ↑ = 1 * + n>0 K * n ↑ . If every entry of the Hamiltonian is a bounded function of time, this series representation converges superexponentially and uniformly [18]. Alternatively, Volterra equations can also be solved numerically [32].

Background
The Bloch-Siegert Hamiltonian, here denoted H BS (t ), is possibly the simplest model to exhibit nontrivial physical effects due to time dependencies in the driving radio-frequency fields. The detailed study of these effects is of paramount importance in the broad field of quantum computing, as they have a deleterious impact on qubit driving and stored quantum information [33]. The Hamiltonian reads as In these expressions, the coupling parameter β is the amplitude of the radio-frequency field. Continuing research over the last decades has produced perturbative expressions for the Bloch-Siegert shifts and evolution operator, starting from Shirley's seminal work [11]. Beyond the rotating-wave approximation, which omits the field's counter-rotating terms and is limited to near-resonant ω ∼ ω 0 ultraweak couplings β/ω 1, one of the most successful approaches used a combination of Floquet formalism and almost degenerate perturbation theory [34]. Still, this could only approximate the temporal dynamics in the vicinity of resonance in the weak-coupling regime β/ω 0.6.
In the case of quantum systems driven by large-amplitude fields β/ω > 1 to β/ω 1 [23], these approaches are no longer sufficient. Yet, such systems are of current fundamental interest, as short associated electromagnetic pulses can manipulate qubits on a large bandwidth. Recently, several methods have thus been designed to overcome the limitations of the theoretical treatment [28,[35][36][37]. These are based on various unitary transformations leading to approximate analytical expressions over an extended driving range. Although these methods are clearly beyond perturbation theory and what Floquet formalism can realistically achieve, they still neglect terms corresponding to multiphoton transitions. Although not dominant, these terms lead to real features in the true Bloch-Siegert dynamics that are as yet unaccounted for [35]. These are visible in qubit driving and time-optimal quantum control experiments, for which determining the Bloch-Siegert dynamics exactly is thus still full of promises [38]. In spite of the theoretical efforts, a nonperturbative truly analytical solution at all orders over the entire coupling range, on and off resonance, is ultimately lacking.

Path-sum solution
Although this is not required by the path-sum method, we pass in the interaction picture to alleviate the notation, yielding the Bloch-Siegert Hamiltonian as Since in the rotating frame, h ↑ (t ) = h ↓ (t ) = 0, the graph K 2 of Fig. 1 has no self-loops and Eqs. (4) and (5) thus give In spite of the apparent divergences in the resonant case ω 0 → ω, the kernels K ↑ and K ↓ are actually well defined in this limit where they simplify to and × e −iωt cos(ωt ).
The peculiar mathematical nature of the resonant limit ω 0 →ω is responsible for the appearance of the term 2ω(t − t ) which is proportional to time in the kernel. The quantity G ↑ as obtained from K ↑ has no closed form, rather it is a hitherto unknown higher special function. It is nonetheless analytically available thanks to the Neumann expansion G ↑ = δ + k K * k ↑ = δ + K ↑ + K ↑ * K ↑ + · · · , which is unconditionally convergent [18]. This observation holds for all N × N time-dependent Hamiltonians treated by path sum.
The Neumann expansion is well suited to analytical computations: observe that at order n of this series, G (n) ↑ = δ + n k=1 K * k is simply given as Equivalently, it is sufficient to integrate the last term of the series at order n − 1, namely, K * n−1 ↑ , to get G (n) ↑ : These integrals are all analytically available and easily accessible: we reached order 13 in a minute on an ordinary laptop treating all parameters as formal variables [39]. Vastly faster computations are achieved upon assigning parameter values before performing the integrals. These calculations give (here displaying the first two orders on resonance ω 0 = ω) Of particular interest for qubit-driving experiments is the evolution of the transition probability P ↑ →↓ (t ) := |U ↓↑ (t )| 2 between the two levels [30,35,40]. This quantity is usually found perturbatively using Floquet theory [11] as Magnus series again suffer from divergences [41]. It is here easily accessible, U ↓↑ being given by Eq. (9). We find that P ↑→↓ (t ) takes on the form of a Fourier-type series with S 2k and C 2k functions of β and t, a representation of which is analytically available (see Appendix B). This form of P ↑→↓ (t ) is due to the path-sum integral of Eq. (9), which resembles a Fourier transform. We emphasize that this is not a general feature of path sum nor of 2 × 2 Hamiltonians, but solely of the present Hamiltonian with linearly polarized driving.
We plot on Fig. 3 the transition probabilities P (3) ↑→↓ (t ), P (7) ↑→↓ (t ), and P (13) ↑→↓ (t ) as calculated analytically from the third, seventh, and thirteenth orders of the Neumann expansion of the exact path-sum solution from the weak-to the ultrastrong-coupling regimes and always in the resonant case ω 0 = ω. Here this situation was chosen because (i) it is mathematically the most difficult to approach exactly owing to the peculiar form of K ↑ which slows down convergence; and (ii) it yields "compact" expressions more suitable for a "concise" presentation (Appendix B). Higher-order terms of the Neumann expansion are readily and analytically available, enabling precise evaluation of P ↑ →↓ (t ) up to any desired target time. Recall that, as discussed above, P (13) ↑→↓ (t ) is actually a single analytical formula involving all even frequencies' sines and cosines up to sin(54ωt ) and cos(54ωt ) with coefficients up to β 54 . We stress here that the Neumann expansion of the path-sum solution is profoundly different from a Taylor . As seen here, each of these formulas are equally valid throughout the coupling regimes, from weak to ultrastrong. Whenever longer times are desired, higher orders of the path-sum solution are readily available analytically. Also shown in (a) and (b) are the second-order Floquet theory [11] (solid green line and green points). The Floquet result is not shown in subsequent figures, where it is wildly inaccurate. Parameters: two-level system driven by the Bloch-Siegert Hamiltonian of Eq. (8) [11,12] starting in the |↑ state at t = 0. series representation, as is, e.g., manifest even at order 0 [see Eq. (12) below and Appendix B].
The fact that the same expression for P (n) ↑ →↓ (t ) is an equally good approximation to the exact transition probability in all parameter regimes, i.e., from β/ω 0 1 to β/ω 0 1 is a signature that the path-sum approach is nonperturbative. For the same reason, we observe that P (n) ↑ →↓ (t ) captures roughly the same number of spin flips in time regardless of β: empirically order 3 reproduces 1-2 flips, order 7 gets 2-3 flips, order 13 captures 4-5. At the opposite, Floquet theory, which is inherently perturbative, only works for β/ω 1 [11], while the diverging Magnus series is limited to very short times.
In Fig. 4 we show the off-resonance ω 0 = ω dynamics of the analytical transition probability P (4) ↑→↓ (t ) obtained from the fourth-order Neumann expansion. Irrespectively of the coupling strength, at any fixed finite order n, P (n) ↑→↓ (t ) is reliable for longer times as we get farther from resonance, for which convergence of the Neumann expansion is slowed by the presence of a linear term in K ↑ . Once again, this is purely a feature of the Bloch-Siegert Hamiltonian and not of the path-sum approach.

Physical insights
Now that we have analytical formulas for the transition probability without the rotating-wave approximation, we may gain insights into the Bloch-Siegert dynamics. For example, we can calculate the spin-flip duration t s f , i.e., the time at which P ↑→↓ (t ) first peaks close to 1 when on resonance ω 0 = ω. Analysis of Eq. (10) with, e.g., the analytic expressions of Appendix B show that C 0 (β, t ) is the dominant contribution to t s f in the weak-coupling regimes β/ω 1 2 , while the C 2k>0 and S 2k functions describe further oscillations smaller by a factor of at least β 2 . Extracting t s f from C 0 leads to Shown here are the numerical solution (dashed black line) and the fully analytical formula for the Neumann expansion of the exact path-sum solution at the fourth order P (4) ↑→↓ (t ) (solid blue line). Parameters: two-level system driven by the Bloch-Siegert Hamiltonian of Eq. (8) [11,12] starting in the |↑ state at t = 0. This is remarkably close to the results obtain from numerical calculations (see Fig. 5). Mathematically, Eq. (11) assumes β/ω < 2 and accurately reflect its discrete jumps are immediately available, however, they cannot be expressed in terms of radicals anymore and are not reproduced here owing to length concerns. Also of interest are the changes affecting the dynamics of the transition probability P ↑→↓ (t ) as β/ω is increased from the weak to strong regimes. For example, in the ultraweakcoupling regime β/ω 1, the path-sum solution reproduces small oscillations around the Floquet calculations which are present in the numerical solution (see Fig. 6). In fact, these small oscillations are already captured by the order 0 of the Neumann expansion of the path-sum solution (!), for which and This shows that the oscillations missed by earlier treatments have a linearly growing amplitude at short times on the order of β 2 t, originate purely from the counter-rotating terms, and never truly vanish as long as β = 0. The diverging parabola in β 2 t 2 reflects the humble beginning of the Rabi oscillation, unsurprisingly missed by order 0. As β is increased, the small oscillations compete with the background Rabi oscillations, thereby giving rise to intricate intermediary effects seen in Fig. 3. This competition also explains why P ↑→↓ (t ) does not always peak at 1, as it results from a complicated superposition of oscillatory terms, in agreement with Eq. (10). We conclude the discussion on physical insights into the Bloch-Siegert dynamics by studying coherent destruction of tunneling (CDT) [10] in the strong coupling β/ω 0 1. This situation is well suited to the use of a general property of Neumann series that allows for arbitrary accelerations of their convergence in the presence of dominant terms [42]. Note that this procedure is always available when expanding pathsum solutions, and is thus not specific to the Bloch-Siegert Hamiltonian.
Concretely, we get a closed-form expression for the evolution operator U(t ) at the zeroth order of the accelerated Neumann expansion of the path-sum solution that leads to perfect or near-perfect fits for any physical quantity of interest both on and off CDT resonances. See Appendix C for details of the calculations. For example, the return probability to the | ↑ state is found to be This formula becomes exact when either ω 0 → 0 or β → 0, as expected from the acceleration procedure. In general, it provides excellent approximations when β/ω 0 is large [see Figs. 7(a)-7(c)]. The remaining integral in P (acc,0) ↑→↑ (t ) has no closed form but can be evaluated explicitly via an infinite series of sines and cosines with Bessel coefficients (Appendix C). This expansion also indicates that the time average of the return probability is which is exactly 1 2 on CDT resonances where J 0 (4β/ω) = 0, consistent with the current understanding of CDT. To be more precise, let us study CDT directly by considering the states |ψ ± = 1 √ 2 (| ↑ ± | ↓ ). The probability of transition between these states, denoted P ψ − →ψ + (t ), is found in the situation where ω 0 (β/ω) 1/2 , as (Appendix C) This expression flawlessly reproduces the numerical solution in its finest details, details which had hitherto not been captured with such accuracy [35]. Minimizing the time average of this formula confirms that the CDT condition is exactly J 0 (4β/ω) = 0, i.e., this is not changed by the nonperturbative corrections. Mathematically, the reason for this is simple: the J 0 function is quadratically dominant over the other terms of the Bessel-series expansion of Eq. (15) because it stems from the sole term of that expansion which does not depend on τ in both integrals. While these results are as expected from the standard theory of CDT, it not so for all physical quantities. Consider, for example, the expectation value of σ x for a system initially prepared in the | ↑ state. As observed by [43], σ x presents anomalous fluctuations on CDT resonances, a fact that was interpreted as a hallmark of and resulting from a crossing Floquet state. This interpretation is in fact not correct. Indeed, at order 0 of the accelerated expansion of the path-sum solution we get (Appendix C) when ω 0 (β/ω) 1/2 , with H 1 (. . .) the first Struve function. Remarkably, the difference n between the location of the nth zero of J 0 (. . .) and of the nth zero of Eq. (17) tends asymptotically to 0 as n ∼ 1/(2π n) for n 1. This asymptotics develops quite quickly: while 1 0.4, already 2 0.03. The fact that the anomalous fluctuations in the expectation value of σ x peak at the zeros of Eq. (17) rather than on CDT resonances is confirmed by the numerical simulations. This analysis indicates that while σ x does indeed seem to fluctuate the most on CDT resonances, it is in fact not true and the phenomenon driving these fluctuations is subtly different from that behind CDT.
These results demonstrate the power of various expansions of the path-sum solution, enabling very precise and hitherto unequaled analytical analysis of subtle phenomena, e.g., P ψ − →ψ + (t ) is on the order of 10 −5 on CDT resonances and is fitted to within machine precision by the formulas provided. This is not because of special features of the Bloch-Siegert Hamiltonian. Rather, the path-sum approach is generally valid for any driving field, as showed by the general solution provided in Sec. III A. This same solution is valid for dissipative non-Hermitian operators [44], and will always be amenable to analytic Neumann and accelerated Neumann expansions, should it lack a closed form.

A. Few-body, (N>2 )-level Hamiltonians
The path-sum approach is by no means limited to two-level systems: e.g., solutions to all time-dependent 3 × 3 and 4 × 4 Hamiltonians are readily available and will be presented in a future work. The number of steps in the exact solution is always finite and the terms involved get progressively simpler because of the "descending ladder principle" [see Fig. 1(e)].
For many-body systems N 1, a further problem appears, namely, the exponential growth in the size of the Hamiltonian.
While path sum does not, in itself, solve the challenges posed by this well-known scaling, it offers tools to manage it via its scale-invariance properties, which we now briefly present as we will use it to treat a many-body molecular system from NMR.

B. Scale invariance
Path sums stem from formal resummations of families of walks. This principle does not depend on what those walks represent. In particular, it remains unchanged by the nature of the evolving system. To exploit this observation, consider a more general type of system histories made of temporal successions of orthogonal vector spacesh : V 1 → V 2 → V 3 . . . . Physically, such histories can describe an evolving subsystem, such as a group of protons in a large molecule. Mathematically, they correspond to walks on a coarse-grained representation of the quantum state space, a subgraphG t of G t . To see this, take a complete family of orthogonal spaces, i.e., i=1 V i = V , where V is the entire quantum state space. To each V i associate a vertex v i and give the edge v i → v j the time-dependent weight P V j · H(t ) · P V i . Here, P V k is the projector onto V k . Observe then that these edge weights are generally non-Abelian. Yet, because path sums fundamentally retain the order and time of the transitions in histories when performing resummations of walks, this setup poses no further difficulty. It follows that the submatrix P V j · U(t , t ) · P V i of the evolution operator is again given as a matrix-valued branched continued fraction of finite depth and breadth. While the shape of this fraction depends on the particular choice of vector spaces, its existence and convergence properties do not. If the vector spaces are chosen so that the shape of the fraction itself is unchanged, and such a choice is always possible, then the path-sum formulation is truly invariant under scale changes in the quantum state space.
An immediate consequence of scale invariance is that there is always a path-sum calculation rigorously relating the global evolution of a system to that of any ensemble of its subsystems, such as clusters of spins in a large molecule (see below). In this scheme, we can evolve each subsystem separately from one another using any preferred method (Magnus, Floquet, path sum, Zassenhaus for short times, etc.), only to then combine these isolated evolutions exactly via a path sum to generate the true system evolution.
While thorough exploitation of the scale-invariance property is beyond the scope of this work, we demonstrate below how it can be used to tackle many-body Hamiltonians, with an emphasis on examples from NMR, i.e., 42 spins coupled by the homonuclear dipolar interaction and spin diffusion under MAS.

V. LARGE MOLECULE IN NMR
We now turn to the general problem of determining the temporal dynamics of spin diffusion as effected by the time-dependent high-field dipolar Hamiltonian for N homonuclear spins: where the interaction amplitude ω i j (t ) is time dependent due to the MAS rotation (see Appendix D for more details). We consider a cationic tin oxo cluster [(MeSn) 12 O 14 (OH) 6 ] 2+ [45] exhibiting N = 42 protons belonging to hydroxyl and methyl groups (see Fig. 8). This structure is idealized and exhibits the main characteristics of already synthesized clusters (distances, angles, crystal packing). The methyl groups are supposed fixed as is the case at low temperature, although this is no requirement of the path-sum method and methyl rotations can be tackled. A single orientation of the molecule toward the principal magnetic field B 0 is considered. Path sum yields analytical expressions for the entries of the evolution operator because the computational complexity of the calculations can be made to be only linear in the system size N depending on the initial state. We stress that this is due primarily to the peculiar structure of the high-field dipolar Hamiltonian, which allows for a particularly efficient usage of the scale invariance and graphical nature of path sums.
In particular, we do not claim to have solved the general many-body problem: there will be Hamiltonians for which this procedure cannot circumvent the exponential explosion of the state space. The methodology we employed is presented below, after the results. In Fig. 8 and in the Supplemental Material, Movie 2, [46], ω r is fixed at 2π × 10 kHz and the initial up spin is located on a hydroxyl proton, denoted H1. During the first 0.15 ms time period (or 1.5 rotor period), an oscillation is observed between two close hydroxyl protons H1 and H2, followed by a partial transfer to the closest methyl group (t 0.15 ms), in particular proton H3. Inside the methyl entity, the frequency of exchange is much faster as the three protons are subjected to much stronger dipolar couplings. In Figs. 9(a) and 9(b) and the Supplemental Material (Movies 1, 3, 4, 5, and 6 for ω r /2π = 5, 20, 40, 60, and 120 kHz [46], the return probability to spin 1 is expressed as a function of ω r and can be described analytically. These results provide an exact justification to recently proposed approximations in the context of the 1 H line dependence under ultrafast MAS [47]. Finally, in Fig. 9(c), strong offsets (roughly 30 ppm at 1.5 GHz, currently the highest magnetic field available for high-resolution solid-state NMR purposes) were added to all protons H i , except the two hydroxyl protons H 1,2 (see inset of Fig. 8 for identification). As the chemical shift offset corresponds simply to I z,i operators, the solution of the spin-diffusion problem remains analytical by using path sum. For strong offsets, spin diffusion is quenched. All of these results are in perfect agreement with experimental observations related to spin diffusion in NMR.

State-space reduction techniques
Simulating many-body quantum systems on classical computers is doomed to be an impossible task, barring the use of approximations. A general class of such approximations, called state-space reduction techniques, bypass the exponential computational hurdle by considering only the most relevant corners of the quantum state space that the system is likely to explore. But, path sum is, first and foremost, a mathematical technique for analytically solving systems of coupled linear differential equations with nonconstant coefficients. This holds regardless of what this system means and how it came about. Therefore, path sum can be used in conjunction with all state-space reduction techniques, as these intervene earlier in selecting the system to be considered.
In this work, which focuses on path sum, we achieve the desired reduction by choosing the initial density matrix ρ(0) to be a pure state with a small number k of up or down spins. Indeed, since the high-field Hamiltonian of Eq. (18) conserves this number at all times, the discrete graph structure G t encoding the quantum state space for path sum consists of exactly N disconnected components, of sizes N k ∼ N k when k N. Hence, the computational cost of finding the evolution operator using a path sum here is O(N k ), i.e., linear in N for a single initial up spin. This procedure is different from approximate state-space truncations approaches [15,16,48] since here the Hamiltonian rigorously enforces the state-space partition. As a result, our calculations retain quantum correlations of up to N spins. More general initial density matrices ρ(0) may be approximated with polynomial cost on expanding them over pure states with k N. In the sector of the quantum space with a single up spin, the difficulty is thus solely due to the time-dependent nature of the Hamiltonian. The evolution operator is then strictly analytical for static experiments and analytically soluble using path sums for MAS experiments. Physically, the time-dependent high-field dipolar Hamiltonian of Eq. (18) implements a continuous-time quantum random walk of the spin on the molecule. This interpretation remains true in the presence of more than one initial up spin, with the caveat that further interactions happen when quantum walkers meet.

Dynamics at the molecular scale
As stated above, the sector of the quantum state space that needs to be considered for an initial pure state with a single up spin is of dimension N. This reduces the problem of calculating the evolution operator to (analytically) solving an N × N system of coupled linear differential equations with nonconstant coefficients. Since, in principle, all pairs of spins interact directly, this system is full. Consequently, if no further partition of the Hamiltonian is used, the graph G t on which path sum is to be implemented is the complete graph on N vertices, which entails a huge (yet finite) number of terms in the path-sum continued fraction. The vast majority of these give negligible contributions to the overall dynamics, however, because of the scales of the interactions involved, one may therefore build up the path-sum continued fraction by brute force, progressively including longer cycles until convergence of the solution is obtained.
An alternative, physically motivated approach appealing once more to scale invariance nonetheless appears preferable as it yields further insights in the temporal dynamics. First, remark that at least one further nontrivial partition of the Hamiltonian is quite natural in the case of the cationic tin oxo cluster, that which puts together all spins belonging to the same methyl or 3 hydroxyls groups. Mathematically,  3 . Edges and self-loops correspond to intergroup and intragroup interactions, respectively. The adjacency matrix of this graph is the 14 × 14 Hamiltonian with 3 × 3 matrix-valued entries evoked in the text. Thus, the shape of the graph is essentially that of the molecule at the methyl and group of 3 hydroxyls level. It comprises two disconnected pathways for spin diffusion corresponding to the opposite sides of the cationic tin oxo cluster which become connected for higher cutoff values > 42.
this is equivalent to seeing the Hamiltonian as a 14 × 14 matrix with matrix-valued entries, each of size 3 × 3. Then, there is a path-sum continued fraction expressing any 3 × 3 block of the global evolution operator U(t , t ) in terms of the "small" Hamiltonians of the corresponding proton groups.
At this point, the path-sum continued fraction is already quite manageable without further approximations, but we can gain additional (analytical) insights into the spin dynamics by removing intergroup interactions weaker than a chosen cutoff value I B,B / , with I B,B the maximum inter-group interaction. Here, B indices mean "block." The value of is itself controlled by convergence of the overall solution. This procedure sends some off-diagonals blocks of the Hamiltonian to 0, givingG t a nontrivial topology which reveals the molecular structure at the methyl and 3 hydroxyls scale, as experienced by the spin excitation during diffusion. See Fig. 10 for an illustrative example, with ω r = 2π × 60 kHz and = 40. The corresponding path-sum continued fraction takes on the topology of the molecule and establishes mathematically the main pathways taken by the spin excitation: where, e.g., H (OH) 3 Me 3 * 3 * H Me 3 Me 2 * 2 * H Me 2 (OH) 3 is the weight of the triangle (OH) 3 Fig. 10(b)]. In these expressions, Id * = 1 * Id 3×3 , the j are given by and j designates the isolated evolution of the jth methyl group, i.e., These results illustrate again the "descending ladder principle" evoked in Fig. 1. Here, all inverses are * inverses and U (OH) 3  The reader may notice that the shape taken by the continued fraction for U (OH) 3 is immediately related to that of the graphG t [ Fig. 10(b)], with each term of the fraction being the weight of a fundamental cycle of the graph. This close, transparent, association between the mathematical form of the solution and the physical problem allows for physically motivated and better controlled approximations. For example, setting 5 to zero so that 4 ≡ 4 in the above solution is immediately understood to mean that one removes the possibility for the spin to diffuse to the remote methyl groups Me5 and Me6 before coming back to the initial group of 3 hydroxyls, an excellent approximation [see Fig. 10(a), red points to be compared to the red dashed line].
Finally, we remark that our choice of partition is not mathematically necessary. For example, larger blocks may be employed equally well or one may form blocks with protons scattered throughout the molecule. In principle, path sum's scale invariance guarantees that any choice, if properly implemented, leads to the same solution. In practice, however, there is a tradeoff between the size of the manipulated blocks and the complexity of the path-sum continued fraction. We do not know in general how to choose the best partition according to this tradeof,f but it seems that physically motivated partitions are a good starting point.

VI. CONCLUSION
In this contribution, we have demonstrated an approach to the problem of finding compact and exact expressions for the evolution operators of quantum dynamical systems driven by time-varying Hamiltonians. As illustrated in Fig. 1, path-sum calculations always involve a "descending ladder" of progressively simpler quantities yielding the exact solution after a finite number of steps. This is in strong contrast with traditional perturbation techniques (Magnus expansion, Floquet theory) which, when carrying out analytically, invariably lead to infinite series and an "ascending ladder" of increasingly intricate quantities, such as Magnus series' nested commutators. Most importantly, the solutions provided by path sums are always analytically accessible, e.g., through Neumann expansions.
As a fundamental and illustrative example, we used path sum to solve the Bloch-Siegert problem, related to the action of the counter-rotating component of the radio-frequency field, at any order. We analytically studied the spin diffusion effected by the homonuclear dipolar coupling Hamiltonian of NMR acting on a large molecule, starting from a pure-state initial density matrix. In general, on many-body systems, we are facing two kinds of "explosive" computational problems: (i) one, quantum in nature, related to the exponential size of the quantum state space; and (ii) one, graph theoretical in nature, related to the time required to construct the path-sum continued fraction, in particular if G t is large and not sparse. Issue (ii) can be managed with partitions and path sum's scale invariance and is further tackled with the implementation of a Lanczos path-sum algorithm [49]. This algorithm naturally exploits matrix sparsity, benefits from path sum's "descending ladder" principle, and was designed with a numerical outlook. It can, in principle, get excellent approximations after only a few iterations, equivalent to truncating a path-sum continued fraction in sufficient depth to reach the desired accuracy. This algorithm is best understood as an extension to time-ordered exponentials of modern numerical procedures for the computation of ordinary matrix exponentials. The first issue (i) is fundamental to quantum mechanics and its management inherently depends on the problem at hand. Here, path sum has the advantage that it works in conjunction with any state-space reduction technique. For the homonuclear dipolar coupling Hamiltonian, we bypassed the problem upon choosing certain initial pure states. The scale invariance of path sum offers further flexibility, as it allows to separately evolve chosen subsystems only to then combine all such evolutions in a globally exact way.
These results call for a discussion on the nature of the solutions sought after by physicists and mathematicians alike. A general assumption seems to be that an acceptable and interesting analytical solution to a problem has been found if and only if it can be presented with a finite number of symbols and preexisting functions. We think this is a restrictive if misleading expectation. For example, a Bessel or a Heun function solution would be considered "satifactory" when both are actually algebraically transcendant, known and understood from the equations they solve and from explicit series expansions involving simpler objects. It seems that at least in some cases our perception of mathematical objects may be biased by facts as simple as their having a name, yet the sine integral function Si(x) = sin(x)/x dx is no more undisputedly analytical than exp[sin(x)/x)]dx. We think that one cannot and should not ask a general purpose analytical method for solving systems of coupled linear differential equations with variable coefficients any more than what there is to be found: (i) finding, in a finite number of steps, an explicit differential or integral equation involving only one unknown function to be determined; and (ii) providing an unconditionally convergent means of expanding the solution as a series of some kind, be it Taylor, Neumann, accelerated Neumann, or other. We may add the requirement that (iii) all calculations should be feasible analytically, i.e., without giving numerical values to all parameters involved. Should one of these criteria fail to be met, a purely numerical strategy would surely be more interesting. But, if all of these demands are indeed satisfied, we may analyze the situation in greater depth and details than possible with numerical computations. This is exemplified by the CDT analysis provided here, the analytical formula for σ x revealing slight deviations from the expected zeros of the Bessel J 0 function.
With these understandings in place, we think that path sum opens an entire new field of research and is now open for the NMR and wider physics communities.
Note added. Recently, it was suggested to us that path sum may be related to Haydock's recursion method for calculating electronic states [51]. While both approaches share the same outlook of recursively resumming Feynman diagrams via path resummations, Haydock's method relies on fundamentally commutative mathematics, in particular determinants, which do not extend to the general setting required by ordered exponentials and scale invariance. Instead, it is possible that lifting Haydock's approach to noncommutativity using Gelfand's quasideterminants [52] would lead to path sum.

ACKNOWLEDGMENTS
C. Bonhomme thanks Dr. F. Ribot for communicating the molecular data pertaining to the cationic tin oxo cluster. P.-L.G. is supported by the Agence Nationale de la Recherche Grant No. ANR-19-CE40-0006. P.-L.G. is also grateful for the financial support from the Royal Commission for the Exhibition of 1851.

APPENDIX A: REMARKS ON THE STATE OF THE ART
While reviewing the state of the art in the course of this work, it appeared to us that a very vast corpus of research had accumulated on quantum dynamics driven by time-varying Hamiltonians. A host of special solutions have been found and numerous betterments of existing techniques have been developed. Some of these are recent enough that we could not cover them in our introduction, such as the flow-equation approach to periodic Hamiltonians [50]. It seems that a proper review article on the subject is urgently needed to gather all results and remedy the pitfalls of our modest introduction.

APPENDIX B: BLOCH-SIEGERT DYNAMICS
In this Appendix, we detail the calculation process for the transition probability P ↑→↓ at order 3 of the Neumann expan-sion of the exact path-sum solution. We work on resonance ω 0 = ω as this yields more compact expressions and also because this situation is the most challenging mathematically. Indeed, precisely when ω = ω 0 the kernel K ↑ has terms that are linear in time and which slow down convergence of P (n) ↑→↓ (t ) to P ↑→↓ (t ) (see Sec. III B 2).
As explained in the main text, at order 3 we have [see Eq. (9) of the main text]. Here, G (3) ↑ is the third-order Neumann expansion of the path-sum solution, i.e., Taken together, these calculations give the transition probability at the third Neumann order as P