Pierre-Louis Giscard
email: giscard@univ-littoral.fr

Christian Bonhomme
email: christian.bonhomme@upmc.fr

Exact dynamics of quantum systems driven by time-varying Hamiltonians: solution for the Bloch-Siegert Hamiltonian and applications to NMR

Comprehending the dynamical behaviour of quantum systems driven by time-varying Hamiltonians is particularly difficult. Systems with as little as two energy levels are not yet fully understood as the usual methods including diagonalisation of the Hamiltonian do not work in this setting. In fact, since the inception of Magnus' expansion in 1954, no fundamentally novel mathematical approach capable of solving the quantum equations of motion with a time-varying Hamiltonian has been devised. We report here of an entirely different non-perturbative approach, termed path-sum, which is always guaranteed to converge, yields the exact analytical solution in a finite number of steps for finite systems and is invariant under scale transformations of the quantum state space. Path-sum can be combined with any state-space reduction technique and can exactly reconstruct the dynamics of a many-body quantum system from the separate, isolated, evolutions of any chosen collection of its sub-systems. As examples of application, we solve analytically for the dynamics of all two-level systems as well as of a many-body Hamiltonian with a particular emphasis on NMR (Nuclear Magnetic Resonance) applications: Bloch-Siegert effect, coherent destruction of tunneling and N -spin systems involving the dipolar Hamiltonian and spin diffusion.

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. ( = 1)

-i H(t ) U(t , t) = d dt U(t , t), (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) [START_REF] Slichter | Principles of Magnetic Resonance[END_REF][START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF][START_REF] Mehring | Object-Oriented Magnetic Resonance[END_REF].

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 diagonalisation of H. Rather, U is formally described by the action of a time-ordering operator on the Dyson series representation of the quantum evolution [START_REF] Dyson | [END_REF], 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 down explicitly so that practically, only low orders of the expansion are useable. 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 timedependent 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 polarised 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 non-avoidable 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 litterature background.

In this contribution, the path-sum method is applied for the very first time 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 non-constant coefficients. So far, no physical applications of these works has been presented. Consequently, they remained unnoticed outside of a specialised mathematical community, and their applicability to long-standing questions pertaining to quantum systems driven by timedependent Hamiltonians went completely unrecognised. 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 the present work is the first fundamentally new 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 emphasise that throughout this work, the time t is and remains a continuous variable, in particular path-sum does not rely on time-discretisation. As a corollary, path-sum yields a non-perturbative formulation of U(t , t), as will be illustrated below with the Bloch-Siegert effect. Further properties of path-sums ensures its scalability to multi-spin 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 sub-systems. In its general form, path-sum is best understood as a method to exactly and analytically solve systems of coupled linear differential equations with non-constant coefficients.

This article 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 time-dependent 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 lin-early polarised 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 scaletransformations 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 Fig 2(a,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. 1). In case another cycle c 2 is accessible to a walker while walking along c 1 , then the sum of all walks will take on the form 1/ 1 -c 1 /(1 -c 2 ) . If instead, both c 1 and c 2 are immediately accessible from the starting point α, the sum of all walks will be 1/(1 -c 1 -c 2 ). Finally, if two cycles c 2 and c 3 with different starting points are both accessible while walking along c 1 , then the sum of all walks will be 1

/ 1 -c 1 /(1 -c 2 ) × 1/(1 -c 3 ) .
There is a unique way to combine these constructions to describe the sum of all walks with chosen starting and ending points on any graph. For example, the walks from α to itself on the graph of Fig. 1 formally sum up to

w walk: α→α w = 1 1 -c 1 1 1-c2 1 1-c3 1 1-c 4 .
This procedure yields any s j |U(t)|s i as branched continued fractions comprising only the weights of the simple cycles and paths of the graph. See Fig 2(c,d,e). Because the graph is finite, there are finitely many such cycles and paths and the fraction is finite in both depth and breadth. It is thus 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 non-commutative convolution-like product

f * g (t , t) := t t f (t , τ )g(τ, t) dτ. (2)
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 super-exponentially convergent Neumann expansion [18] 1 * + n>0 f * n .

III. TWO-LEVEL SYSTEMS WITH TIME-DEPENDENT HAMILTONIANS

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

H(t) = h ↑ (t) h ↑↓ (t) h ↓↑ (t) h ↓ (t) . (3)
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 transcendantal 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. [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF] we get:

U(t , t) ↑↑ = t t G ↑ (τ, t)dτ, U(t , t) ↓↓ = t t G ↓ (τ, t)dτ, (4)
U(t , t) ↓↑ = -i t t τ0 t τ0 τ1 δ(τ 0 -τ 2 ) -ih ↓ (τ 0 ) e -i τ 0 τ 2 h ↓ (τ3)dτ3 × × h ↓↑ (τ 2 )G ↑ (τ 1 , t) dτ 2 dτ 1 dτ 0 , U(t , t) ↑↓ = -i t t τ0 t τ0 τ1 δ(τ 0 -τ 2 ) -ih ↑ (τ 0 ) e -i τ 0 τ 2 h ↑ (τ3)dτ3 × × h ↑↓ (τ 2 )G ↓ (τ 1 , t) dτ 2 dτ 1 dτ 0 .
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 ↑ G ↑ (t , t) = δ(t , t)

+ t t K ↑ (t , τ ) G ↑ (τ, t) dτ, (5)
and similarly for G ↓ . The kernel K ↑ of the above equation Step by step evaluation of GK 2 ,11(t , t) (dashed rectangle) showing the finite character of the continued fraction. The sum is performed on the simple cycles (C) of length 1 and 2 (these are indicated in red-other edges are indicated in dashed grey lines). At each step of the continued fraction, a vertex is removed (grey cross) and we work on subgraphs of less and less complexity. The calculation of entry U21 from GK 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. For 2 × 2 Hamiltonians there are only 2 steps on the ladder, i.e. all continued fractions stop at depth 2. For all 3 × 3 and 4 × 4 Hamiltonians the continued fractions stop at depth 3 and 4, respectively.

is

K ↑ (t , t) = -ih ↑ (t ) (6)
- t t t τ1 h ↑↓ (t ) δ(τ 2 -τ 1 ) -ih ↓ (τ 2 ) e -i τ 2 τ 1 h ↓ (τ3)dτ3 × × h ↓↑ (τ 1 )dτ 2 dτ 1 .
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 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 super-exponentially and uniformly [18]. Alternatively, Volterra equations can also be solved numerically [START_REF] Hackbusch | Integral Equations: Theory and Numerical Treatment[END_REF].

B. Bloch-Siegert dynamics

Background

The Bloch-Siegert Hamiltonian, here denoted H BS (t), is possibly the simplest model to exhibit non-trivial phys-ical effects due to time-dependencies in the driving radiofrequency 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 [START_REF] Zhang | [END_REF]. The Hamiltonian reads

H BS (t) = ω 0 /2 2β cos(ωt) 2β cos(ωt) -ω 0 /2 . (7)
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 ultra-weak 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 over-come 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 multi-photon 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 non-perturbative 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 2) has no self-loops and Eqs (4-5) thus give

H BS (t) = 2β cos(ωt) cos(ω 0 t)σ x -2β cos(ωt) sin(ω 0 t)σ y . (8) Since in the rotating frame, h ↑ (t) = h ↓ (t) = 0, the graph K 2 of Fig. (
U(t) ↑↑ = t 0 G ↑ (τ, 0)dτ, U(t) ↓↓ = t 0 G ↓ (τ, 0)dτ, U(t) ↓↑ = (9) -2iβ t 0 τ1 0 cos(ωτ 1 )e iω0τ1 G ↑ (τ 0 , 0)dτ 0 dτ 1 , U(t) ↑↓ = -2iβ t 0 τ1 0 cos(ωτ 1 )e -iω0τ1 G ↓ (τ 0 , 0)dτ 0 dτ 1 , while G ↑ (t , t) = (1 * -K ↑ ) * -1 , G ↓ = 1 * -K ↓ ) * -1 with K ↑ (t , t) = 4β 2 ω 2 -ω 2 0 cos(ωt ) k ↑ (t)e -iω0(t -t) -k ↑ (t ) , K ↓ (t , t) = iβ 2 ω 2 -ω 2 0 1 + e -2iωt k ↓ (t ) -k ↓ (t)e i(ω+ω0)(t -t) ,
where k ↑ (t) = iω 0 cos(ωt) + ω sin(ωt) and k ↓ (t) = e 2iωt (ω + ω 0 ) -(ω -ω 0 ). 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

K ↑ (t , t) = β 2 ω ie 2iωt -ie 2iωt -2ω(t -t) e -iωt cos(ωt ),
and

K ↓ (t , t) = β 2 ω -i + ie 2i(t -t)ω + 2ω(t -t )e 2iωt e -iωt cos(ωt ).
The peculiar mathematical nature of the resonant limit ω 0 → ω is responsible for the apparence 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 [18]. This observation holds for all N × N time-dependent Hamiltonians treated by path-sum.

G ↑ = δ + k K * k ↑ = δ + K ↑ + K ↑ * K ↑ + • • • , which is unconditionally convergent
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 G (n) ↑ (t, 0) = t 0 K ↑ (t, τ )G (n-1) ↑ (τ, 0)dτ.
Equivalently, it is sufficient to integrate the last term of the series at order n -1, namely

K * n-1 ↑ , to get G (n) ↑ : G (n) ↑ (t, 0) = t 0 K ↑ (t, τ )K * n-1 ↑ (τ, 0)dτ + G (n-1) (t, 0).
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 [START_REF]The Mathematica notebook generating these calculations is[END_REF]. Vastly faster computations are achieved upon assigning parameter values before performing the integrals. This calculations give (here displaying the first two orders on resonance

ω 0 = ω), G ↑ (t, 0) = δ(t) - β 2 ω e -iωt cos(ωt) -ie 2iωt + 2ωt + i + β 4 24ω 3 e -3iωt cos(ωt) 3ie 6iωt + 6e 4iωt (-2iω 2 t 2 + 2ωt + i) + e 2iωt (8ω 3 t 3 + 12iω 2 t 2 + 12ωt -15i) -12ωt + 6i + • • •
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,[START_REF] Angelo | [END_REF]. 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-like series

P ↑→↓ (t) = ∞ k=0 sin(2kωt)S 2k (β, t) + cos(2kωt)C 2k (β, t), (10)
P ↑↦↓ (t) (a) (b) (c) (d) (e) (f) (g) (h) (i)
↑→↓ (t) (dotted blue line) and P (13) ↑→↓ (t) (solid blue line), see Appendix B. As seen here, each of these formulas are equally valid throughout the coupling regimes, from weak to ultra-strong. Whenever longer times are desired, higher orders of the path-sum solution are readily available analytically. Also shown in figures (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.

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.

Visualizing the solution

Calculating G ↑ up to a finite order n as indicated earlier G ↑ ≡ G and t 4n+2-2k , respectively. Finally, we found analytically that at all orders n ≥ 0, P (n) ↑→↓ (0) = 0 as expected, although this is non-trivial to check. For t large enough P We plot on Fig. 3 the transition probabilities P ↑→↓ (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 ultra-strong 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 ↑→↓ (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.

are readily and analytically available, enabling precise evaluation of P ↑ →↓ (t) up to any desired target time. Recall that, as discussed above, P

↑→↓ (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 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 non-perturbative. 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. [START_REF] Dyson | [END_REF] we show the off-resonance ω 0 = ω dynamics of the analytical transition probability P ↑→↓ (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 novel insights into the Bloch-Siegert dynamics. For example, we can calculate the spin-flip duration t sf , 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 shows that C 0 (β, t) is the dominant contribution to t sf 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 sf from C 0 leads to

t sf = (11) 1 2 √ 2 12 β 2 - 15 ω 2 + √ 3 β 4 ω 2 91β 8 -88β 6 ω 2 + 16β 4 ω 4 = 1 β 1 2 3 + √ 3 - β 8ω 2 1 2 129 + 67 √ 3 - β 3 128ω 4 1 2 16131 + 5545 √ 3 + O β 4 .
This is remarkably close to the results obtain from numerical calculations, see Fig. (5). Mathematically, Eq. ( 11) assumes β/ω < 2 1 91 11 -√ 30 0.49. Beyond this point the above estimate yields a complex number as the real solution switches to another root of the derivative of C 0 . 11) (solid red line) and fully numerical results (blue dots). Discontinuities in the numerical results are due to small oscillations of P ↑→↓ (t), which make t sf undergo discrete jumps as one wins over the others. These are well captured analytically by a more advanced analysis including the S 2k and C 2k>0 functions. In insets: two examples of time evolution of P ↑→↓ (t), the straight red lines are the predictions of Eq. ( 11).

Even better analytical formulas for t sf with domains of validity that go much further into the stronger coupling regimes 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 ultra-weak coupling 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 G (0) ↑ = δ(t , t) and

P (0) ↑→↓ (t) = β 2 t ω sin(2ωt) + β 2 2ω 2 + β 2 t 2 - β 2 2ω 2 cos(2ωt). ( 12
) 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 counterrotating terms, and never trully 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. [START_REF] Mehring | Object-Oriented Magnetic Resonance[END_REF]. This competition also explains why P ↑→↓ (t) does not always peak at 1, as it results from a complicated superposition ↑→↓ (t) of Eq. ( 12) (solid blue line), second order Floquet theory [11] (solid green line and green points) and the numerical solution (dashed black line).

of oscilatory 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 [START_REF] Giscard | Journal of Integral Equations and Applications[END_REF]. Note that this procedure is always available when expanding path-sum 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 0th 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 P (acc,0)

↑→↑ (t) = cos 2β ω sin(ωt) + e -1 2 itω0 -1 (13)
+ t 0 iω 0 e -1 2 iτ ω0 sin 2 β ω sin(ωτ ) -sin(ωt) dτ 2 ,
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 Fig. (7 a,b,c). The remaining integral in P (acc,0) ↑→↑ (t) has no closed form but can be evaluated explicitely 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 (t) as given by Eq. ( 13) (solid blue line), the numerical solution (dashed black line), and its predicted time-average Eq. ( 14) (solid red straight line, indistinguishable from the numerically computed time-average). Middle line: transition probability P ψ -→ψ + (t) for a system starting in the |ψ-state at t = 0 with: (d) 4β/ω = 2.404..., first zero of J0(4β/ω); (e) 4β/ = 11.79..., fourth zero of J0(4β/ω); and (f ) 4β/ω = 27.49..., ninth zero of J0(4β/ω). Note the changes of scales. Shown here are the formula of Eq. ( 15) (solid blue line) and the numerical solution (dashed black line), these two being completely indistinguishable. Bottom line: far off-resonance ω = 100ω0 expectation value of σx for a system starting in the | ↑ state at t = 0 with: (g) 4β/ω = 2.404..., first zero of J0(4β/ω); (h) 4β/ω = 11.79..., fourth zero of J0(4β/ω); and (i) 4β/ω = 27.49..., ninth zero of J0(4β/ω). Note the changes of scales in σx . Shown here are the formula of Eq. ( 16) (solid blue line) and the numerical solution (dashed black line), these two being indistinguishable.

P (acc,0) ↑→↑ (t) t = 1 2 1 + J 0 4β ω , (14)
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)

P (acc,0) ψ-→ψ+ (t) =
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 non-perturbative 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 [START_REF] Thorwart | [END_REF], σ x presents anomalous fluctuations on CDT resonances, a fact that was interpreted as a hallmark of and resulting from a crossing Floquet states. 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 , σ x (acc,0) = ω 0 t 0 sin 4β ω sin(ωτ ) dτ. (16)
This simple expression fits once again absolutely flawlessly with the numerically computed expectation σ x , see Fig. (7 d,e,f). Now evaluating the integral remaining in Eq. ( 16) via Bessel functions shows that the time average of σ x is

σ x (acc,0) t = 2ω 0 ω ∞ n=0 J 2n+1 4β ω 1 2n + 1 ,
whose extrema are reached whenever

1 - π 2 H H H 1 4β ω = 0, (17)
with H H H 1 (.) the first Struve function. Remarquably, 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 zeroes 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 formulae 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 §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.

IV. FEW-TO MANY-BODY HAMILTONIANS A. Few-body, N>2 -level Hamiltonians

The path-sum approach is by no mean limited to twolevel 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. 2 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 spaces h :

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 subgraph Gt 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 timedependent weight P Vj .H(t) .P Vi . 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 Vj .U(t , t) .P Vi 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 homonuculear dipolar inter- action 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:

H II = i,j 1 2 ω ij (t) 3I iz I jz -I i . I j , (18)
where the interaction amplitude ω ij (t) is time-dependent due to the MAS rotation, see 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 idealised and exhibits the main characteristics of already synthesised 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 towards 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 Movie 2 (See Supplemental Material at [URL will be inserted by publisher] for this movie), ω r is fixed at 2π × 10 kHz and the initial up-spin is located on a hydroxyl proton, denoted H1. During the first 0.15ms 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.15ms), 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 Fig. 9(a,b) and Movies 1, 3, 4, 5 and 6 for ω r /2π = 5, 20, 40, 60 and 120 kHz ((See Supplemental Material at [URL will be inserted by publisher] for these movies), 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 ultra-fast MAS [46]. 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.

(ms) |〈 ↑ z,i ψ (t)〉 2 (a) (b) (c)
Figure 9. Return probability for the spin excitation on the initial hydroxyl proton H1 as a function of ωR (one plot point every 20Hz): (a) after a fixed time of t = 0.05ms, a situation exhibiting numerous peaks for small ωR values that are not all resolved on this picture; and (b) after two rotor periods t = 2 × (2π/ωR). (c) Probability of finding the spin excitation on hydroxyl protons H1 (solid red line) or H2 (dashed blue line) as a function of time for ωR = 10kHz and with a very strong offset of roughly 30 ppm at 1.5 GHz on all protons except H1 and H2. The total probability of being either on H1 or H2 (dotted black line) never goes below 0.94 over 3 rotor periods.

A. Setting up the path-sum: methodology

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 pathsum is, first and foremost, a mathematical technique for analytically solving systems of coupled linear differential equations with non-constant 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 the present 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,47], since here the Hamiltonian rigorously enforces the statespace 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 timedependent 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 non-constant 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 pathsum 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, 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 inter-group interactions weaker than a chosen cut-off 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, giving Gt a non-trivial 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:

U (OH) 3 = t t Id * + iH (OH) 3 + H (OH) 3 Me1 * Γ 1 * H Me1(OH) 3 + H (OH) 3 Me2 * Σ 2 * H Me2(OH) 3 + H (OH) 3 Me3 * Σ 3 * H Me3(OH) 3 -i H (OH) 3 Me2 * Γ 2 * H Me2Me3 * Σ 3 * H Me3(OH) 3 -i H (OH) 3 Me3 * Σ 3 * H Me3Me2 * Σ 2 * H Me2(OH) 3 * -1 (τ, t)dτ,
where e.g. [START_REF] Mehring | Object-Oriented Magnetic Resonance[END_REF] on Gt (Fig. 10 b). In these expressions, Id * = 1 * Id 3×3 , the Σ j are given by

H (OH) 3 Me3 * Σ 3 * H Me3Me2 * Σ 2 * H Me2(OH) 3 is the weight of the triangle (OH) 3 → Me 2 → Me 3 → (OH)
Σ 2 = 1 Id * + iH Me2 + H Me2Me3 * Σ 3 * H Me3Me2 , Σ 3 = 1 Id * + iH Me3 + H Me3Me4 * Σ 4 * H Me4Me3 , Σ 4 = 1 Id * + iH Me4 + H Me4Me5 * Σ 5 * H Me5Me4 , Σ 5 = 1 Id * + iH Me5 + H Me5Me6 * Γ 6 * H Me6Me5 ,
and Γ j designates the isolated evolution of the jth methyl group, i.e.

Γ j = 1 Id * + iH Mej .
These results illustrate again the "descending ladder principle" evoked in Fig. 2. Here, all inverses are *inverses and U (OH) 3 is the 3 × 3 block of the global evolution operator giving the probability amplitudes over a group of 3 hydroxyls. H Mex and H (OH) 3 are the Hamiltonians of isolated methyl and of a group of 3 hydroxyls, respectively. Similarly, H MeiMej is the interaction between neighbouring methyls and H Mei(OH) 3 the interaction between a methyl and a group of 3 hydroxyls. The reader may notice that the shape taken by the continued fraction for U (OH) 3 is immediately related to that of the graph Gt (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 trade-off 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 tradeoff but it seems that physically motivated partitions are a good starting point.

VI. CONCLUSION

In this contribution, we have demonstrated an entirely novel 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 Figure 2, 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 carying 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 Edges and self-loops correspond to inter-and intra-group 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 cut-off values Λ > 42.

tackled with the implementation of a Lanczos path-sum algorithm [48]. 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 pathsum 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/interesting analytical solution to a problem has been found if and only if it can be presented with a finite number of symbols and pre-existing functions. We think this is a restrictive if missleading 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 per-ception of mathematical objects may be biased by facts as simple as their having a name, yet the sine integral function Si(x) = sin(x)/xdx 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 mean 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 criterion fail to be met, a purely numerical strategy would surely be more interesting. But if all of these demands are indeed satisfied, we may analyse 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 zeroes of the Bessel J 0 function. With these understandings in place, we think that path-sum opens an entire new field of research is now open for the NMR and wider physics communities. All calculations were performed analytically on Math-ematica. The notebook generating these results, as well

Figure 1 .

 1 Figure 1. Graph illustrating the use of path-sum.

Figure 2 .

 2 Figure 2. The path-sum continued fraction for the exact calculation of the entries of U(t , t) is always of finite depth and breadth. (a) Illustrative example of a 2 × 2 time-dependent Hamiltonian H(t), involved for instance in spin dynamics. (b) Dynamical graph Gt = K2 with adjacency matrix H(t). Circles correspond to self-loops (diagonal terms of H(t)) and directed edges to off-diagonal terms. The associated weights are the entries of H(t). (c) Evolution operator U(t , t) as seen by path-sum, with integrals of the GK 2 ,ii quantities. (d)Step by step evaluation of GK 2 ,11(t , t) (dashed rectangle) showing the finite character of the continued fraction. The sum is performed on the simple cycles (C) of length 1 and 2 (these are indicated in red-other edges are indicated in dashed grey lines). At each step of the continued fraction, a vertex is removed (grey cross) and we work on subgraphs of less and less complexity. The calculation of entry U21 from GK 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. For 2 × 2 Hamiltonians there are only 2 steps on the ladder, i.e. all continued fractions stop at depth 2. For all 3 × 3 and 4 × 4 Hamiltonians the continued fractions stop at depth 3 and 4, respectively.

Figure 3 .

 3 Figure 3. Bloch-Siegert dynamics: resonant ω0 = ω transition probability P ↑ →↓ (t) as a function of time in the weak to ultra-strong coupling regimes: with (a) β/ω = 0.5, (b) β/ω = 0.7, (c) β/ω = 0.9, (d) β/ω = 1.2, (e) β/ω = 1.6, (f ) β/ω = 2, (g) β/ω = 3.5, (h) β/ω = 5 and (i) β/ω = 15. Shown here are the numerical solution (dashed black line) and the fully analytical formulas for the Neumann expansions of the exact path-sum solution P

( 3 )

 3 ↑→↓ (t) (dotted red line), P

  (n) ↑ , yields an expression P

  (n) ↑→↓ (t) which includes all terms of Eq. (10) up to sin (4n + 2)ωt and cos (4n + 2)ωt , while S (n) 2k≤4n+2 and C (n) 2k≤4n+2 are polynomials in β and t including up to β 4n+2 and t 4n+3-2k

  ↑→↓ (t) may diverge: truncated path-sums are not necessarily unitary.

( 7 )

 7 ↑→↓ (t) and P

Figure 4 .

 4 Figure 4. Bloch-Siegert dynamics: off-resonance ω0 = ω transition probability P ↑ →↓ (t) as a function of time in the weak to strong coupling regimes: with (a) ω0 = 2ω, β/ω = 0.7, (b) ω0 = 2ω, β/ω = 1.6, (c) ω0 = 2ω, β/ω = 3.5, (d) ω0 = 8ω, β/ω = 0.7, (e) ω0 = 8ω, β/ω = 1.6, and (f ) ω0 = 8ω, β/ω = 3.5. 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 )

 4 ↑→↓ (t) obtained from the fourth order Neumann expansion. Irrespectively of the coupling strength, at any fixed finite order n, P

Figure 5 .

 5 Figure 5. Bloch-Siegert dynamics: Resonant ω0 = ω spinflip duration t sf , at which P ↑→↓ (t) first peaks at or close to 1, as a function of the coupling strength β. Shown here are the exact formula of Eq. (11) (solid red line) and fully numerical results (blue dots). Discontinuities in the numerical results are due to small oscillations of P ↑→↓ (t), which make t sf undergo discrete jumps as one wins over the others. These are well captured analytically by a more advanced analysis including the S 2k and C 2k>0 functions. In insets: two examples of time evolution of P ↑→↓ (t), the straight red lines are the predictions of Eq. (11).

Figure 6 .

 6 Figure 6. Bloch-Siegert dynamics: resonant ω0 = ω transition probability P ↑→↓ (t) in the ultra-weak coupling regime β/ω = 0.05 for short times. Shown here are the order 0 formula P (0)

Figure 7 .

 7 Figure 7. Coherent destruction of tunneling: Top line, return probability P ↑→↑ (t) in the ultra-strong coupling regime β/ω0 = 30 for (a) ω = 4ω0; (b) ω = 20ω0; and (c) ω = 100ω0. Shown here are P (acc,0) ↑→↑

  ) -sin(ωτ ) dτ 2 .

Figure 8 .

 8 Figure 8. Analytical spin-diffusion on a cationic tin oxo-cluster with N = 42 protons (shown in inset) submitted to the time-dependent high-field dipolar Hamiltonian under MAS (rotor angular velocity ωr = 2π × 10kHz). The figure shows the time evolution of the probability | ψ(t)| ↑z,i | 2 of finding a spin-up along z on proton i for three protons: a hydroxyl proton H1 (solid red line), on which the excitation starts; a nearby hydroxyl proton H2 (dashed blue line); and a methyl proton H3 (dot-dashed black line).

Figure 10 .

 10 Figure 10. Building the path-sum on the cationic tin oxo-cluster. (a) Probability of return of the spin excitation on a hydroxyl proton H1 shown on Fig.8as a function of time for ωr = 2π × 60 kHz over 3 rotor periods: (i) solution with no cut-off (solid black line, identical with Λ > 100), (ii) analytical approximation with low interaction cut-off Λ = 40 (dashed red line) and (iii) further approximation obtained upon setting Σ5 to zero (red points). (b) Discrete structure Gt of the quantum state space as seen by path-sum when Λ = 40 and corresponding to the equations given in the text for U (OH) 3 . Edges and self-loops correspond to inter-and intra-group 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 cut-off values Λ > 42.

ACKNOWLEDGMENTS

C. Bonhomme thanks Dr. F. Ribot for communicating the molecular data pertaining to the cationic tin oxocluster. P.-L. Giscard is supported by the Agence Nationale de la Recherche grant ANR-19-CE40-0006. P.-L. Giscard is also grateful for the financial support from the Royal Commission for the Exhibition of 1851 over the period 2015-2018, during which time the present research was started.

Appendix A: Remarks on the state-of-the-art While reviewing the state-of-the-art in the course of the present 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 [49]. 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.

Following publication of the preprint of the present article, it was suggested to us that path-sum may be related to Haydock's recursion method for calculating electronic states [50]. 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 non-commutativity using Gelfand's quasi-determinants [START_REF]Gel'fand and V. Retakh[END_REF] would lead to path-sum.

Appendix B: Bloch Siegert dynamics

In this appendix, we detail the calculation process for the transition probability P ↑→↓ at order 3 of the Neumann expansion 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

As explained in the main text, at order 3 we have P

↑ (τ 0 , 0)dτ 0 dτ 1 , see Eq. ( 9) of the main text. Here G

↑ is the third order Neumann expansion of the path-sum solution, i.e.

G

(3)

Taken together, these calculations give the transition probability at the third Neumann order as Suppose that we are given a function or matrix of two times K(t , t) = K 1 (t , t) + K 2 (t , t) such that in some sense K 1 is much larger than K 2 . Suppose further that we are interested in the solution of the linear Volterra integral equation of the second kind G(t , t) = δ+K * G = δ(t -t) + t t K(t , τ )G(τ, t)dτ , as will always be the case when expanding the exact path-sum solution to quantum dynamical problems at any scale.

Instead of expanding G as usual, G(t , t) = δ+ n K * n , one can exploit the fact that K 1 is dominant over K 2 to accelerate convergence of the Neumann expansion by expressing G in terms of the solutions G i of the "individual" Volterra equations

, see [START_REF] Giscard | Journal of Integral Equations and Applications[END_REF] for details. Since T * 0 = δ(t -t), the 0th order term of the accelerated expansion is then simply the * -product of the solutions of the individual Volterra equations:

This is particularly well suited to physical situations where a certain parameter dominates over the others: not only because the so-obtained expression for G is greatly improved, but also because in general the individual G i are known exactly. Furthermore, this acceleration procedure continues to hold for any number of kernels K i [START_REF] Giscard | Journal of Integral Equations and Applications[END_REF].

Taking the Bloch-Siegert Hamiltonian of Eq. ( 7) as an example, let us use path-sum's scale invariance to work in the trivial situation where we have single subsystem, namely the entire system itself. Then, we get that

with G the solution of the matrix-valued linear integral Volterra equation of the second kind with matrix kernel

The ultra-strong coupling regime β/ω 0 1 thus corresponds to the situation described above as K 1 dominates K 2 . Since furthermore both G i are immediately accessible as

we get G (acc,0) easily and integrating it with respect to t yields Higher orders of the accelerated expansion of the pathsum solution are also available although they are not necessary given the machine-precision accuracy with respect to numerical solutions already reached by order 0. The integrals in U (acc,0) (t) have no closed form but can be determined exactly via standard expansions over Bessel functions since e.g.

The modulus squared of U (acc,0) (t) 11 gives Eq. ( 13) of the main text, while calculating other quantities such as P ψ-→ψ+ (t) and σ x from U (acc,0) (t) is now a simple task, giving e.g.

In the regime ω 0 ω, both cos(ω 0 t/2) and sin(ω 0 t/4) are essentially equal to their initial t = 0 values, leading to Eq. ( 16). We consider the time-dependent high-field dipolar Hamiltonian of Eq. ( 18) presented in the main text, with interaction terms under MAS

where r ij is the distance between protons i and j and [START_REF] Slichter | Principles of Magnetic Resonance[END_REF] ξ ij (t) := 2 √ 2 sin(ψ ij ) cos(ψ ij ) sin(φ ij + ω r t)

+ sin(ψ ij ) 2 cos(2φ ij + 2ω r t).

In this expression, ψ ij is the angle between ij and the z-axis and φ ij is the angle between ij and the x-axis for a coordinate system fixed to the sample. Finally, ω r is the angular velocity of the rotor.

The raw molecular data pertaining to the cationic tin oxo-cluster is available online on the webpage http://www-lmpa.univ-littoral.fr/ ~plgiscard/ and included here as a dataset.