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

TL;DR

This paper introduces a novel non-perturbative path-sum method for exactly solving the dynamics of quantum systems with time-varying Hamiltonians, demonstrated on two-level and many-body NMR systems.

Contribution

The paper presents the path-sum approach, a new mathematical technique that guarantees convergence and provides exact solutions for finite quantum systems with time-dependent Hamiltonians.

Findings

Exact solutions for two-level quantum systems.

Analytical results for NMR-related many-body Hamiltonians.

Demonstrated applicability to Bloch-Siegert effect and spin diffusion.

Abstract

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.

Equations

$$-i\mathsf{H}(t^{\prime })\mathsf{U}(t^{\prime },t)=\frac{d}{d{t}^{\prime }}\mathsf{U}(t^{\prime },t),$$

$$\sum_{w\,\text{walk:}\,\alpha\to\alpha}\!\!\!\!w\,=\frac{1}{1-c_{1}\frac{1}{1-{\color[rgb]{0.0,0.5,0.0}c_{2}}}\frac{1}{1-{\color[rgb]{0,0,1}c_{3}}\frac{1}{1-{\color[rgb]{1,0,0}c_{4}}}}}.$$

$$\big{(}f\ast g\big{)}(t^{\prime},t):=\int_{t}^{t^{\prime}}\!\!f(t^{\prime},\tau)g(\tau,t)\,d\tau.$$

$$\mathsf{H}(t)=\left(\begin{array}{cc}{h}_{\uparrow }(t)& h_{\uparrow \downarrow }(t)\\ h_{\downarrow \uparrow }(t)& h_{\downarrow }(t)\end{array}\right).$$

$$\mathsf{U}(t^{\prime },t{)}_{\uparrow \uparrow }={\int }_t^{t^{\prime }}{G}_{\uparrow }(\tau ,t)d\tau ,\mathsf{U}(t^{\prime },t{)}_{\downarrow \downarrow }={\int }_t^{t^{\prime }}{G}_{\downarrow }(\tau ,t)d\tau ,$$

$$\mathsf{U}(t^{\prime },t{)}_{\downarrow \uparrow }=$$

$$\displaystyle-i\int_{t}^{t^{\prime}}\!\!\!\!\int_{t}^{\tau_{0}}\!\!\!\int_{\tau_{1}}^{\tau_{0}}\!\!\!\Big{(}\delta(\tau_{0}-\tau_{2})-ih_{\downarrow}(\tau_{0})\,e^{-i\!\int_{\tau_{2}}^{\tau_{0}}h_{\downarrow}(\tau_{3})d\tau_{3}}\!\Big{)}\times$$

$$\times h_{\downarrow \uparrow }({\tau }_2)G_{\uparrow }({\tau }_1,t)d{\tau }_{2}d{\tau }_{1}d{\tau }_0,$$

$$\mathsf{U}(t^{\prime },t{)}_{\uparrow \downarrow }=$$

$$\displaystyle-i\int_{t}^{t^{\prime}}\!\!\!\!\int_{t}^{\tau_{0}}\!\!\!\int_{\tau_{1}}^{\tau_{0}}\!\!\!\Big{(}\delta(\tau_{0}-\tau_{2})-ih_{\uparrow}(\tau_{0})\,e^{-i\!\int_{\tau_{2}}^{\tau_{0}}h_{\uparrow}(\tau_{3})d\tau_{3}}\!\Big{)}\times$$

$$\times h_{\uparrow \downarrow }({\tau }_2)G_{\downarrow }({\tau }_1,t)d{\tau }_{2}d{\tau }_{1}d{\tau }_0.$$

$$G_{\uparrow }(t^{\prime },t)=\delta (t^{\prime },t)+{\int }_t^{t^{\prime }}{K}_{\uparrow }(t^{\prime },\tau )G_{\uparrow }(\tau ,t)d\tau ,$$

$$K_{\uparrow }(t^{\prime },t)$$

$$\displaystyle-\!\!\int_{t}^{t^{\prime}}\!\!\int_{\tau_{1}}^{t^{\prime}}\!\!h_{\uparrow\downarrow}(t^{\prime})\Big{(}\delta(\tau_{2}-\tau_{1})-ih_{\downarrow}(\tau_{2})\,e^{-i\int_{\tau_{1}}^{\tau_{2}}h_{\downarrow}(\tau_{3})d\tau_{3}}\Big{)}\times$$

$$\times h_{\downarrow \uparrow }({\tau }_1)d{\tau }_{2}d{\tau }_1.$$

$${\mathsf{H}}_{BS}(t)=\left(\begin{array}{cc}{\omega }_0/2& 2\beta \cos (\omega t)\\ 2\beta \cos (\omega t)& -{\omega }_0/2\end{array}\right).$$

$${\mathsf{H}}_{BS}(t)=2\beta \cos (\omega t)\cos ({\omega }_{0}t){\sigma }_x-2\beta \cos (\omega t)\sin ({\omega }_{0}t){\sigma }_y.$$

$$\mathsf{U}(t{)}_{\uparrow \uparrow }={\int }_0^t{G}_{\uparrow }(\tau ,0)d\tau ,\mathsf{U}(t{)}_{\downarrow \downarrow }={\int }_0^t{G}_{\downarrow }(\tau ,0)d\tau ,$$

$$\mathsf{U}(t{)}_{\downarrow \uparrow }=$$

$$-2i\beta {\int }_0^t{\int }_0^{{\tau }_1}\cos (\omega {\tau }_1)e^{i{\omega }_0{\tau }_1}{G}_{\uparrow }({\tau }_0,0)d{\tau }_{0}d{\tau }_1,$$

$$\mathsf{U}(t{)}_{\uparrow \downarrow }=$$

$$-2i\beta {\int }_0^t{\int }_0^{{\tau }_1}\cos (\omega {\tau }_1)e^{-i{\omega }_0{\tau }_1}{G}_{\downarrow }({\tau }_0,0)d{\tau }_{0}d{\tau }_1,$$

$$K_{\uparrow }(t^{\prime },t)$$

$$K_{\downarrow }(t^{\prime },t)$$

$$\frac{i{\beta }^2}{{\omega }^2-{\omega }_0^2}\left(1+e^{-2i\omega t^{\prime }}\right)\left(k_{\downarrow }(t^{\prime })-k_{\downarrow }(t)e^{i(\omega +{\omega }_0)(t^{\prime }-t)}\right),$$

$$K_{\uparrow }(t^{\prime },t)$$

$$\frac{{\beta }^2}{\omega }\left(i{e}^{2i\omega t^{\prime }}-i{e}^{2i\omega t}-2\omega (t^{\prime }-t)\right)e^{-i\omega t^{\prime }}\cos (\omega t^{\prime }),$$

$$K_{\downarrow }(t^{\prime },t)$$

$$\frac{{\beta }^2}{\omega }\left(-i+i{e}^{2i(t^{\prime }-t)\omega }+2\omega (t-t^{\prime })e^{2i\omega t^{\prime }}\right)e^{-i\omega t^{\prime }}\cos (\omega t^{\prime }).$$

$$G_{\uparrow }^{(n)}(t,0)={\int }_0^t{K}_{\uparrow }(t,\tau )G_{\uparrow }^{(n-1)}(\tau ,0)d\tau .$$

$$G_{\uparrow }^{(n)}(t,0)={\int }_0^t{K}_{\uparrow }(t,\tau )K_{\uparrow }^{\ast n-1}(\tau ,0)d\tau +G^{(n-1)}(t,0).$$

$$G_{\uparrow }(t,0)=\delta (t)$$

$$\displaystyle+\frac{\beta^{4}}{24\omega^{3}}e^{-3i\omega t}\cos(\omega t)\Big{(}3ie^{6i\omega t}$$

$$+6e^{4i\omega t}(-2i{\omega }^2{t}^2+2\omega t+i)$$

$$+e^{2i\omega t}(8{\omega }^3{t}^3+12i{\omega }^2{t}^2+12\omega t-15i)$$

$$\displaystyle\hskip 14.22636pt-12\omega t+6i\Big{)}$$

$$+\cdots$$

$$P_{\uparrow \to \downarrow }(t)=\sum _{k=0}^{\infty }\sin (2k\omega t)S_{2k}(\beta ,t)+\cos (2k\omega t)C_{2k}(\beta ,t),$$

$$t_{sf}$$

$$\frac{1}{2\sqrt{2}}\sqrt{\frac{12}{{\beta }^2}-\frac{15}{{\omega }^2}+\frac{\sqrt{3}}{{\beta }^4{\omega }^2}\sqrt{91{\beta }^8-88{\beta }^6{\omega }^2+16{\beta }^4{\omega }^4}}$$

$$=\frac{1}{\beta }\sqrt{\frac{1}{2}\left(3+\sqrt{3}\right)}-\frac{\beta }{8{\omega }^2}\sqrt{\frac{1}{2}\left(129+67\sqrt{3}\right)}$$

$$-\frac{{\beta }^3}{128{\omega }^4}\sqrt{\frac{1}{2}\left(16131+5545\sqrt{3}\right)}+O\left({\beta }^4\right).$$

$$P_{\uparrow \to \downarrow }^{(0)}(t)=\frac{{\beta }^{2}t}{\omega }\sin (2\omega t)+\frac{{\beta }^2}{2{\omega }^2}+{\beta }^2{t}^2-\frac{{\beta }^2}{2{\omega }^2}\cos (2\omega t).$$

Full text

Time-varying Hamiltonian $|$ Path-sum $|$ Analytical and numerical methods $|$ Bloch-Siegert effect $|$ Nuclear Magnetic Resonance $|$ $N$-spins systems and dipolar couplings

I Introduction

The unitary evolution operator $\mathsf{U}(t^{\prime},t)$ describing the time dynamics of quantum systems is defined as the unique solution of Schrödinger’s equation with quantum Hamiltonian $\mathsf{H}$, i.e. $(\hbar=1)$

[TABLE]

and such that $\mathsf{U}(t^{\prime}=t,\,t)=\mathsf{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) Slichter (1990); Ernst et al. (1987); Mehring and Weberruß (2001). As is typically the case in NMR, the Hamiltonian may be time-dependent and might furthermore not commute with itself at various times, $\mathsf{H}(t)\mathsf{H}(t^{\prime})-\mathsf{H}(t^{\prime})\mathsf{H}(t)\neq\mathbf{0}$ for $t^{\prime}\neq 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 $\mathsf{H}$. Rather, $\mathsf{U}$ is formally described by the action of a time-ordering operator on the Dyson series representation of the quantum evolution Dyson (1952), a formulation which does not permit concrete calculations to be carried out.

As a consequence, only approximate expressions of $\mathsf{U}(t^{\prime},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 Waugh et al. (1968). This relies exclusively on the Magnus expansion Magnus (1954) of $\mathsf{U}(t^{\prime},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 Magnus (1954) and Fel’dman Fel’dman (1984). 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 $\mathsf{U}(t,0)=\mathsf{P}(t)\exp(\mathsf{F}t)$, with $\mathsf{P}(t)$ a periodic time-dependent matrix and $\mathsf{F}$ a constant matrix, both of which are determined perturbatively when working analytically Blanes et al. (2009), or otherwise via numerical procedures Großmann and Hänggi (1992); Grifoni and Hänggi (1998). Floquet formalism was first used by Shirley Shirley (1965), 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 Bloch and Siegert (1940).

We also mention numerical methods: (i) Fer and Magnus-Floquet hybrids proposed recently as potential expansions for the evolution operator Takegoshi et al. (2015); Mananga (2018), (ii) Zassenhaus and Suzuki-Trotter propagator approximations Brüschweiler and Ernst (1997); Dumez et al. (2010); Mentink-Vigier et al. (2017). 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\times 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 $\mathsf{U}(t^{\prime},t)$. The rigorous underpinnings of this approach were laid out in Giscard et al. (2015, 2012) 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 time-dependent 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 $\mathsf{U}(t^{\prime},t)$ as the inverse of an operator with respect to certain $\ast$-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 $\mathsf{U}(t^{\prime},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 $\mathsf{U}(t^{\prime},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 Giscard et al. (2015) 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\times 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 linearly 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 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\leavevmode\nobreak\ :\leavevmode\nobreak\ |s_{1}\rangle\mapsto|s_{2}\rangle\mapsto|s_{3}\rangle\cdots$, each transition $|s_{i}\rangle\mapsto|s_{i+1}\rangle$ 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}\rangle\mapsto|s_{i+1}\rangle$ is dictated by the Hamiltonian as $\langle s_{i+1}|\mathsf{H}(t_{i})|s_{i}\rangle$.

A natural representation of such discrete histories is as walks on a graph. To see this, let $\mathcal{G}_{t}$ be the graph such that each vertex $v_{i}$ corresponds to one member $|s_{i}\rangle$ of an orthonormal basis for the entire state-space and give the directed edge $v_{i}\mapsto v_{j}$ the time-dependent weight $\langle s_{j}|\mathsf{H}(t)|s_{i}\rangle$. In this picture, a system history as defined earlier is a walk on $\mathcal{G}_{t}$ and $\mathsf{H}(t)$ is the adjacency matrix of $\mathcal{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 $\langle s_{j}|\mathsf{U}(t)|s_{i}\rangle$ of the evolution operator $\mathsf{U}(t)$ is given by the sum over all walks from $v_{i}$ to $v_{j}$ on $\mathcal{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 Giscard et al. (2012), in particular they can be uniquely factored into products of prime walks: the simple cycles ($\mathcal{C}$) and paths ($\mathcal{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 $\alpha$ to itself are of the form $c_{1}^{n}$, i.e. $c_{1}$ repeated $n$ times. Therefore the sum of all such walks is formally $\sum_{n}c_{1}^{n}=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/\big{(}1-c_{1}/(1-c_{2})\big{)}. If instead, both $c_{1}$ and $c_{2}$ are immediately accessible from the starting point $\alpha$, 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/\big{(}1-c_{1}/(1-c_{2})\times 1/(1-c_{3})\big{)}. 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 $\alpha$ to itself on the graph of Fig. 1 formally sum up to

[TABLE]

This procedure yields any $\langle s_{j}|\mathsf{U}(t)|s_{i}\rangle$ 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^{\prime},t)$ that multiply via a non-commutative convolution-like product

[TABLE]

This means that for general time-dependent Hamiltonians the continued fraction formulation for $\mathsf{U}(t)$ involves products and resolvent with respect to the $\ast$ multiplication and that the order of traversal of the edges along the cycles must be respected. The $\ast$-resolvents such as $(1_{\ast}-f)^{\ast-1}$ with $1_{\ast}\equiv\delta(t^{\prime}-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 $\ast$-resolvent can be always represented analytically by the super-exponentially convergent Neumann expansion Giscard et al. (2015) $1_{\ast}+\sum_{n>0}f^{\ast n}$.

III Two-level systems with time-dependent Hamiltonians

III.1 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 Kayanuma (1994); Blanes et al. (2009); Xie and Hai (2010); Irish et al. (2005); Ashhab et al. (2007); Saiko and Fedoruk (2008); Gangopadhyay et al. (2010); Rapedius (2015); Barnes and Das Sarma (2012); Yan et al. (2015); Schmidt (2018). 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 Barnes and Das Sarma (2012); Zeuch et al. (2018). The most general two-level Hamiltonian is

[TABLE]

In this expression we only require that $h_{\downarrow\uparrow}(t)$, $h_{\uparrow\downarrow}(t)$, $h_{\uparrow}(t)$ and $h_{\downarrow}(t)$ be bounded functions of time over the interval $[t,t^{\prime}]$ of interest. So far, no analytical expression has been found for the corresponding evolution operator $\mathsf{U}(t)$ defined as the unique solution of Eq. (1) with the Hamiltonian of Eq. (3). It is known that particular choices for $\mathsf{H}(t)$ lead to evolution operators that involve higher transcendantal functions Xie and Hai (2010); Braak (2011). Thus the best possible result for the general $\mathsf{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 time-dependent two-level systems. Following the process of Fig. (2) we get:

[TABLE]

while the ‘usual’ $\mathsf{U}(t)$ is actually $\mathsf{U}(t,0)$. The two-times Green’s functions $G_{\uparrow}:=(1_{\ast}-K_{\uparrow})^{\ast-1}$ and $G_{\downarrow}:=(1_{\ast}-K_{\downarrow})^{\ast-1}$ solve linear Volterra equations of the second kind, e.g. for $G_{\uparrow}$

[TABLE]

and similarly for $G_{\downarrow}$. The kernel $K_{\uparrow}$ of the above equation is

[TABLE]

while kernel $K_{\downarrow}$ entering $G_{\downarrow}$ 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 Xie and Hai (2010)–the solution is at least analytically available from its Neumann expansion; in the case Eq. (5) it is $G_{\uparrow}=1_{\ast}+\sum_{n>0}K_{\uparrow}^{\ast n}$. If every entry of the Hamiltonian is a bounded function of time, this series representation converges super-exponentially and uniformly Giscard et al. (2015). Alternatively, Volterra equations can also be solved numerically Hackbusch (1995).

III.2 Bloch-Siegert dynamics

III.2.1 Background

The Bloch-Siegert Hamiltonian, here denoted $\mathsf{H}_{BS}(t)$, is possibly the simplest model to exhibit non-trivial 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 Zhang et al. (2018). The Hamiltonian reads

[TABLE]

In these expressions, the coupling parameter $\beta$ 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 Shirley (1965). Beyond the rotating wave approximation–which omits the field’s counter-rotating terms and is limited to near resonant $\omega\sim\omega_{0}$ ultra-weak couplings $\beta/\omega\ll 1$–one of the most successful approaches used a combination of Floquet formalism and almost degenerate perturbation theory Aravind and Hirschfelder (1984). Still, this could only approximate the temporal dynamics in the vicinity of resonance in the weak coupling regime $\beta/\omega\lesssim 0.6$.

In the case of quantum systems driven by large amplitude fields $\beta/\omega>1$ to $\beta/\omega\gg 1$ Ashhab et al. (2007), 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 Lü and Zheng (2012); Yan et al. (2015); Zhang and Chen (2015); Yan et al. (2017). 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 Lü and Zheng (2012). 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 Laucht et al. (2016). 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.

III.2.2 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

[TABLE]

Since in the rotating frame, $h_{\uparrow}(t)=h_{\downarrow}(t)=0$, the graph $K_{2}$ of Fig. (2) has no self-loops and Eqs (4–5) thus give

[TABLE]

while $G_{\uparrow}(t^{\prime},t)=(1_{\ast}-K_{\uparrow})^{\ast-1}$, G_{\downarrow}=\big{(}1_{\ast}-K_{\downarrow})^{\ast-1} with

[TABLE]

where $k_{\uparrow}(t)=i\omega_{0}\cos(\omega t)+\omega\sin(\omega t)$ and $k_{\downarrow}(t)=e^{2i\omega t}(\omega+\omega_{0})-(\omega-\omega_{0})$. In spite of the apparent divergences in the resonant case $\omega_{0}\to\omega$, the kernels $K_{\uparrow}$ and $K_{\downarrow}$ are actually well defined in this limit where they simplify to

[TABLE]

The peculiar mathematical nature of the resonant limit $\omega_{0}\to\omega$ is responsible for the apparence of the term $2\omega(t^{\prime}-t)$ which is proportional to time in the kernel.

The quantity $G_{\uparrow}$ as obtained from $K_{\uparrow}$ has no closed form, rather it is a hitherto unknown higher special function. It is nonetheless analytically available thanks to the Neumann expansion $G_{\uparrow}=\delta+\sum_{k}K_{\uparrow}^{\ast k}=\delta+K_{\uparrow}+K_{\uparrow}\ast K_{\uparrow}+\cdots$, which is unconditionally convergent Giscard et al. (2015). This observation holds for all $N\times 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)}_{\uparrow}=\delta+\sum_{k=1}^{n}K^{\ast k}$ is simply given as

[TABLE]

Equivalently, it is sufficient to integrate the last term of the series at order $n-1$, namely $K_{\uparrow}^{\ast n-1}$, to get $G_{\uparrow}^{(n)}$:

[TABLE]

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 111The Mathematica notebook generating these calculations is available for download at http://www-lmpa.univ-littoral.fr/~plgiscard/. Vastly faster computations are achieved upon assigning parameter values before performing the integrals. This calculations give (here displaying the first two orders on resonance $\omega_{0}=\omega$),

[TABLE]

Of particular interest for qubit-driving experiments is the evolution of the transition probability $P_{\uparrow\mapsto\downarrow}(t):=|U_{\downarrow\uparrow}(t)|^{2}$ between the two levels Angelo and Wreszinski (2005); Lü and Zheng (2012); Zeuch et al. (2018). This quantity is usually found perturbatively using Floquet theory Shirley (1965) as Magnus series again suffer from divergences Maricq (1987). It is here easily accessible–$\mathsf{U}_{\downarrow\uparrow}$ being given by Eq. (9). We find that $P_{\uparrow\to\downarrow}(t)$ takes on the form of a Fourier-like series

[TABLE]

with $S_{2k}$ and $C_{2k}$ functions of $\beta$ and $t$, a representation of which is analytically available (see Appendix B). This form of $P_{\uparrow\to\downarrow}(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\times 2$ Hamiltonians, but solely of the present Hamiltonian with linearly polarized driving.

III.2.3 Visualizing the solution

Calculating $G_{\uparrow}$ up to a finite order $n$ as indicated earlier $G_{\uparrow}\equiv G_{\uparrow}^{(n)}$, yields an expression $P^{(n)}_{\uparrow\to\downarrow}(t)$ which includes all terms of Eq. (10) up to \sin\big{(}(4n+2)\omega t\big{)} and \cos\big{(}(4n+2)\omega t\big{)}, while $S^{(n)}_{2k\leq 4n+2}$ and $C^{(n)}_{2k\leq 4n+2}$ are polynomials in $\beta$ and $t$ including up to $\beta^{4n+2}$ and $t^{4n+3-2k}$ and $t^{4n+2-2k}$, respectively. Finally, we found analytically that at all orders $n\geq 0$, $P^{(n)}_{\uparrow\to\downarrow}(0)=0$ as expected, although this is non-trivial to check. For $t$ large enough $P^{(n)}_{\uparrow\to\downarrow}(t)$ may diverge: truncated path-sums are not necessarily unitary.

We plot on Fig. 3 the transition probabilities $P^{(3)}_{\uparrow\to\downarrow}(t)$, $P^{(7)}_{\uparrow\to\downarrow}(t)$ and $P^{(13)}_{\uparrow\to\downarrow}(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 $\omega_{0}=\omega$. Here this situation was chosen because: i) it is mathematically the most difficult to approach exactly owing to the peculiar form of $K_{\uparrow}$ 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_{\uparrow\mapsto\downarrow}(t)$ up to any desired target time. Recall that, as discussed above, $P^{(13)}_{\uparrow\to\downarrow}(t)$ is actually a single analytical formula involving all even frequencies sines and cosines up to $\sin(54\omega t)$ and $\cos(54\omega t)$ with coefficients up to $\beta^{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)}_{\uparrow\mapsto\downarrow}(t)$ is an equally good approximation to the exact transition probability in all parameter regimes, i.e. from $\beta/\omega_{0}\ll 1$ to $\beta/\omega_{0}\gg 1$ is a signature that the path-sum approach is non-perturbative. For the same reason, we observe that $P^{(n)}_{\uparrow\mapsto\downarrow}(t)$ captures roughly the same number of spin flips in time regardless of $\beta$: 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 $\beta/\omega\ll 1$ Shirley (1965), while the diverging Magnus series is limited to very short times.

In Fig. (4) we show the off-resonance $\omega_{0}\neq\omega$ dynamics of the analytical transition probability $P^{(4)}_{\uparrow\to\downarrow}(t)$ obtained from the fourth order Neumann expansion. Irrespectively of the coupling strength, at any fixed finite order $n$, $P^{(n)}_{\uparrow\to\downarrow}(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_{\uparrow}$. Once again, this is purely a feature of the Bloch-Siegert Hamiltonian and not of the path-sum approach.

III.2.4 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_{\uparrow\to\downarrow}(t)$ first peaks close to 1 when on resonance $\omega_{0}=\omega$. Analysis of Eq. (10) with e.g. the analytic expressions of Appendix B shows that $C_{0}(\beta,t)$ is the dominant contribution to $t_{sf}$ in the weak coupling regimes $\beta/\omega\lesssim 1/2$, while the $C_{2k>0}$ and $S_{2k}$ functions describe further oscillations smaller by a factor of at least $\beta^{2}$. Extracting $t_{sf}$ from $C_{0}$ leads to

[TABLE]

This is remarkably close to the results obtain from numerical calculations, see Fig. (5). Mathematically, Eq. (11) assumes $\beta/\omega<2\sqrt{\frac{1}{91}\left(11-\sqrt{30}\right)}\simeq 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}$.

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_{\uparrow\to\downarrow}(t)$ as $\beta/\omega$ is increased from the weak to strong regimes. For example, in the ultra-weak coupling regime $\beta/\omega\ll 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_{\uparrow}^{(0)}=\delta(t^{\prime},t)$ and

[TABLE]

This shows that the oscillations missed by earlier treatments have a linearly-growing amplitude at short times on the order of $\beta^{2}t$, originate purely from the counter-rotating terms, and never trully vanish as long as $\beta\neq 0$. The diverging parabola in $\beta^{2}t^{2}$ reflects the humble beginning of the Rabi oscillation, unsurprisingly missed by order 0. As $\beta$ 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_{\uparrow\to\downarrow}(t)$ does not always peak at 1, as it results from a complicated superposition 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) Grifoni and Hänggi (1998) in the strong coupling $\beta/\omega_{0}\gg 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 Giscard (2020). 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 $\mathsf{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 $|\!\uparrow\rangle$ state is found to be

[TABLE]

This formula becomes exact when either $\omega_{0}\to 0$ or $\beta\to 0$, as expected from the acceleration procedure. In general, it provides excellent approximations when $\beta/\omega_{0}$ is large, see Fig. (7 a, b, c).

The remaining integral in $P^{(acc,0)}_{\uparrow\to\uparrow}(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

[TABLE]

which is exactly $1/2$ on CDT resonances where $J_{0}(4\beta/\omega)=0$, consistent with the current understanding of CDT. To be more precise let us study CDT directly by considering the states $|\psi_{\pm}\rangle=\frac{1}{\sqrt{2}}(|\!\uparrow\rangle\pm|\!\downarrow\rangle)$. The probability of transition between these states, denoted $P_{\psi_{-}\to\psi_{+}}\!(t)$, is found in the situation where $\omega_{0}\ll(\beta/\omega)^{1/2}$, as (Appendix C)

[TABLE]

This expression flawlessly reproduces the numerical solution in its finest details, details which had hitherto not been captured with such accuracy Lü and Zheng (2012). Minimizing the time-average of this formula confirms that the CDT condition is exactly $J_{0}(4\beta/\omega)=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 $\tau$ 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 $\sigma_{x}$ for a system initially prepared in the $|\!\uparrow\rangle$ state. As observed by Thorwart et al. (2000), $\langle\sigma_{x}\rangle$ 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 $\omega_{0}\ll(\beta/\omega)^{1/2}$,

[TABLE]

This simple expression fits once again absolutely flawlessly with the numerically computed expectation $\langle\sigma_{x}\rangle$, see Fig. (7 d, e, f). Now evaluating the integral remaining in Eq. (16) via Bessel functions shows that the time average of $\langle\sigma_{x}\rangle$ is

[TABLE]

whose extrema are reached whenever

[TABLE]

with $\boldsymbol{H}_{1}(.)$ the first Struve function. Remarquably, the difference $\Delta_{n}$ between the location of the $n$th zero of $J_{0}(.)$ and of the $n$th zero of Eq. (17) tends asymptotically to 0 as $\Delta_{n}\sim 1/(2\pi n)$ for $n\gg 1$. This asymptotics develops quite quickly: while $\Delta_{1}\simeq 0.4$, already $\Delta_{2}\simeq 0.03$. The fact that the anomalous fluctuations in the expectation value of $\sigma_{x}$ peak at the zeroes of Eq. (17) rather than on CDT resonances is confirmed by the numerical simulations. This analysis indicates that while $\langle\sigma_{x}\rangle$ 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_{\psi_{-}\to\psi_{+}}(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.1. This same solution is valid for dissipative non-Hermitian operators Sergi and Zloshchastiev (2013), and will always be amenable to analytic Neumann and accelerated Neumann expansions, should it lack a closed form.

IV Few- to many-body Hamiltonians

IV.1 Few-body, $N\!\!>\!\!2\,$-level Hamiltonians

The path-sum approach is by no mean limited to two-level systems: e.g. solutions to all time-dependent $3\times 3$ and $4\times 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\gg 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.

IV.2 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 $\tilde{h}:\leavevmode\nobreak\ V_{1}\mapsto V_{2}\mapsto V_{3}\cdots$. 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 $\tilde{\mathcal{G}}_{t}$ of $\mathcal{G}_{t}$. To see this, take a complete family of orthogonal spaces, i.e. $\bigoplus_{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}\mapsto v_{j}$ the time-dependent weight $\mathsf{P}_{V_{j}}.\mathsf{H}(t)\,.\mathsf{P}_{V_{i}}$. Here $\mathsf{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 $\mathsf{P}_{V_{j}}.\mathsf{U}(t^{\prime},t)\,.\mathsf{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 homonuculear 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:

[TABLE]

where the interaction amplitude $\omega_{ij}(t)$ is time-dependent due to the MAS rotation, see D for more details. We consider a cationic tin oxo-cluster \big{[}(\text{MeSn})_{12}\text{O}_{14}(\text{OH})_{6}\big{]}^{2+} Vivas-Reyes et al. (2002) 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), $\omega_{r}$ is fixed at 2\pi\times 10\leavevmode\nobreak\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\gtrsim 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 Fig. 9(a,b) and Movies 1, 3, 4, 5 and 6 for $\omega_{r}/2\pi=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 $\omega_{r}$ and can be described analytically. These results provide an exact justification to recently proposed approximations in the context of the 1H line dependence under ultra-fast MAS Sternberg et al. (2018). 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 Hi, except the two hydroxyl protons H1,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.

V.1 Setting up the path-sum: methodology

V.1.1 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 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 $\rho(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 $\mathcal{G}_{t}$ encoding the quantum state space for path-sum consists of exactly $N$ disconnected components, of sizes $\binom{N}{k}\sim N^{k}$ when $k\ll N$. Hence, the computational cost of finding the evolution operator using a path-sum here is O\big{(}N^{k}), i.e. linear in $N$ for a single initial up-spin. This procedure is different from approximate state space truncations approaches Brüschweiler and Ernst (1997); Butler et al. (2009); Dumez et al. (2010), 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 $\rho(0)$ may be approximated with polynomial cost on expanding them over pure states with $k\ll 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.

V.1.2 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\times 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 $\mathcal{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 non-trivial 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\times 14$ matrix with matrix valued entries, each of size $3\times 3$. Then there is a path-sum continued fraction expressing any $3\times 3$ block of the global evolution operator $\mathsf{U}(t^{\prime},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^{\prime}}/\Lambda$, with $I_{B,B^{\prime}}$ the maximum inter-group interaction. Here $B$ indices mean ”block”. The value of $\Lambda$ is itself controlled by convergence of the overall solution. This procedure sends some off-diagonals blocks of the Hamiltonian to 0, giving $\tilde{\mathcal{G}}_{t}$ 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 $\omega_{r}=2\pi\times 60$ kHz and $\Lambda=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:

[TABLE]

where e.g. $\mathsf{H}_{\text{(OH)}_{3}\text{Me}_{3}}\ast\Sigma_{3}\ast\mathsf{H}_{\text{Me}_{3}\text{Me}_{2}}\ast\Sigma_{2}\ast\mathsf{H}_{\text{Me}_{2}\text{(OH)}_{3}}$ is the weight of the triangle $\text{(OH)}_{3}\mapsto\text{Me}_{2}\mapsto\text{Me}_{3}\mapsto\text{(OH)}_{3}$ on $\tilde{\mathcal{G}}_{t}$ (Fig. 10 b). In these expressions, $\mathsf{Id}_{\ast}=1_{\ast}\mathsf{Id}_{3\times 3}$, the $\Sigma_{j}$ are given by

[TABLE]

and $\Gamma_{j}$ designates the isolated evolution of the $j$th methyl group, i.e.

[TABLE]

These results illustrate again the “descending ladder principle” evoked in Fig. 2. Here, all inverses are $\ast$-inverses and $\mathsf{U}_{\text{(OH)}_{3}}$ is the $3\times 3$ block of the global evolution operator giving the probability amplitudes over a group of 3 hydroxyls. $\mathsf{H}_{\text{Me}_{x}}$ and $\mathsf{H}_{\text{(OH)}_{3}}$ are the Hamiltonians of isolated methyl and of a group of 3 hydroxyls, respectively. Similarly, $\mathsf{H}_{\text{Me}_{i}\text{Me}_{j}}$ is the interaction between neighbouring methyls and $\mathsf{H}_{\text{Me}_{i}\text{(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 $\mathsf{U}_{\text{(OH)}_{3}}$ is immediately related to that of the graph $\tilde{\mathcal{G}}_{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 $\Sigma_{5}$ to zero so that $\Sigma_{4}\equiv\Gamma_{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 trade-off 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 $\mathcal{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 Giscard and Pozza (2019). 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/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 perception of mathematical objects may be biased by facts as simple as their having a name, yet the sine integral function $\text{Si}(x)=\int\sin(x)/xdx$ is no more undisputedly analytical than $\int\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 $\langle\sigma_{x}\rangle$ 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.

Acknowledgements.

C. Bonhomme thanks Dr. F. Ribot for communicating the molecular data pertaining to the cationic tin oxo-cluster. 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 Vogl et al. (2019). 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 Haydock and Ashcroft (1991). 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 I. Gel’fand and V. Retakh (1997) would lead to path-sum.

Appendix B Bloch Siegert dynamics

In this appendix, we detail the calculation process for the transition probability $P_{\uparrow\to\downarrow}$ at order 3 of the Neumann expansion of the exact path-sum solution. We work on resonance $\omega_{0}=\omega$ as this yields more compact expressions and also because this situation is the most challenging mathematically. Indeed, precisely when $\omega=\omega_{0}$ the kernel $K_{\uparrow}$ has terms that are linear in time and which slow down convergence of $P^{(n)}_{\uparrow\to\downarrow}(t)$ to $P_{\uparrow\to\downarrow}(t)$ (see §III.2.2).

As explained in the main text, at order 3 we have $P^{(3)}_{\uparrow\to\downarrow}(t)=|U^{(3)}(t)_{\downarrow\uparrow}|^{2}$ with

[TABLE]

see Eq. (9) of the main text. Here $G^{(3)}_{\uparrow}$ is the third order Neumann expansion of the path-sum solution, i.e.

[TABLE]

Taken together, these calculations give the transition probability at the third Neumann order as

[TABLE]

in accordance with Eq. (10) of the main text. Here we have

[TABLE]

Now on to the $C^{(3)}_{2k}$ functions:

[TABLE]

All calculations were performed analytically on Mathematica. The notebook generating these results, as well as any desired higher order of the Neumann expansion of the exact path-sum solution is available for download at http://www-lmpa.univ-littoral.fr/~plgiscard/. Everytime the order is increased by one, e.g. from $P^{(3)}_{\uparrow\to\downarrow}(t)$ to $P^{(4)}_{\uparrow\to\downarrow}(t)$, each expression above gains new high order terms while four new functions also appear, e.g. $S^{(4)}_{16}(\beta,t)$, $S^{(4)}_{18}(\beta,t)$, $C^{(4)}_{16}(\beta,t)$ and $C^{(4)}_{18}(\beta,t)$ all enter $P^{(4)}_{\uparrow\to\downarrow}(t)$).

Appendix C Accelerated Neumann series

Suppose that we are given a function or matrix of two times $K(t^{\prime},t)=K_{1}(t^{\prime},t)+K_{2}(t^{\prime},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^{\prime},t)=\delta+K\ast G=\delta(t^{\prime}-t)+\int_{t}^{t^{\prime}}K(t^{\prime},\tau)G(\tau,t)d\tau$, 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^{\prime},t)=\delta+\sum_{n}K^{\ast 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 $G_{i}=\delta+K_{i}\ast G_{i}$. More specifically one gets

[TABLE]

where $T=\delta(t^{\prime}-t)-G_{1}\ast G_{2}+G_{1}\ast G_{2}\ast(K_{1}+K_{2})$, see Giscard (2020) for details. Since $T^{\ast 0}=\delta(t^{\prime}-t)$, the 0th order term of the accelerated expansion is then simply the $\ast$-product of the solutions of the individual Volterra equations:

[TABLE]

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}$ Giscard (2020).

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

[TABLE]

with $G$ the solution of the matrix-valued linear integral Volterra equation of the second kind with matrix kernel $\mathsf{K}=\mathsf{K}_{1}+\mathsf{K}_{2}$, where

[TABLE]

The ultra-strong coupling regime $\beta/\omega_{0}\gg 1$ thus corresponds to the situation described above as $K_{1}$ dominates $K_{2}$. Since furthermore both $\mathsf{G}_{i}$ are immediately accessible as

[TABLE]

we get $\mathsf{G}^{(acc,0)}$ easily and integrating it with respect to $t$ yields

[TABLE]

Higher orders of the accelerated expansion of the path-sum 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 $\mathsf{U}^{(acc,0)}(t)$ have no closed form but can be determined exactly via standard expansions over Bessel functions since e.g.

[TABLE]

The modulus squared of $\mathsf{U}^{(acc,0)}(t)_{11}$ gives Eq. (13) of the main text, while calculating other quantities such as $P_{\psi_{-}\to\psi_{+}}(t)$ and $\langle\sigma_{x}\rangle$ from $\mathsf{U}^{(acc,0)}(t)$ is now a simple task, giving e.g.

[TABLE]

In the regime $\omega_{0}\ll\omega$, both $\cos(\omega_{0}t/2)$ and $\sin(\omega_{0}t/4)$ are essentially equal to their initial $t=0$ values, leading to Eq. (16).

Appendix D Interaction terms in the high-field dipolar Hamiltonian

We consider the time-dependent high-field dipolar Hamiltonian of Eq. (18) presented in the main text, with interaction terms under MAS

[TABLE]

where $r_{ij}$ is the distance between protons $i$ and $j$ and Slichter (1990)

[TABLE]

In this expression, $\psi_{ij}$ is the angle between $\vec{ij}$ and the $z$-axis and $\phi_{ij}$ is the angle between $\vec{ij}$ and the $x$-axis for a coordinate system fixed to the sample. Finally, $\omega_{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.