An integral-free representation of the Dyson series using divided differences

The Dyson series is an infinite sum of multi-dimensional time-ordered integrals, which serves as a formal representation of the quantum time-evolution operator in the interaction-picture. Using the mathematical tool of divided differences, we introduce an alternative representation for the series that is entirely free from both time ordering and integrals. In this new formalism, the Dyson expansion is given as a sum of efficiently-computable divided differences of the exponential function, considerably simplifying the calculation of the Dyson expansion terms, while also allowing for time-dependent perturbation calculations to be performed directly in the Schrödinger-picture. We showcase the utility of this novel representation by studying a number of use cases. We also discuss several immediate applications.


I. INTRODUCTION
The Dyson series [1] is one of the fundamental results of quantum scattering theory.It is a perturbative expansion of the quantum time-evolution operator in the interaction-picture, wherein each summand is formally represented as a multidimensional integral over a time-ordered product of the interaction-picture Hamiltonian at different points in time.As such, the series serves as a pivotal tool for studying transition properties of timedependent quantum many-body systems [2].The Dyson series also has a close relation to Feynman diagrams in quantum field theory [3] and to the Magnus expansion in the theory of first-order homogeneous linear differential equation for linear operators [4].
Despite its fundamental role in quantum theory, one of the key challenges in using the Dyson series in practical applications remains the evaluation of the multi-dimensional integrals over products of time-ordered operators, making the calculation of the terms in the series an exceedingly complicated task [5].In this paper, we derive an analytical, closed-form expression for the summands of the Dyson series, by explicitly evaluating the Dyson integrals.We accomplish this by utilizing the machinery of 'divided differences' -a mathematical tool normally used for computing tables of logarithms and trigonometric functions and for calculating the coefficients in the interpolation polynomial in the Newton form [6][7][8][9][10][11][12].
The main technical contribution of this work is an alternative, yet equivalent, formulation of the Dyson series wherein the summands are given in terms of efficiently computable divided differences of the exponential function.We argue that our novel representation, which is devoid of integrals and time ordering, makes it a very useful tool in the study of perturbation effects in many-body quantum systems -a fundamental branch of quantum physics.
In particular, this representation allows us to write an explicit, integral-free time-dependent perturbation expansion for the time-evolution operator in the Schrödinger-picture enabling us to carry perturbation calculations without the need to switch to the interaction-picture.To showcase the utility of our formulation, we work out a number of examples, for which we explicitly calculate the first few terms of the series.We begin by a brief introduction of divided differences, followed by the derivation of the integral-free representation of the Dyson series.

II. DIVIDED DIFFERENCES -A BRIEF OVERVIEW
The divided differences [7,8], which will be a major ingredient in the derivation detailed in what follows, is a recursive division process.The divided differences of a function f (•) is defined as with respect to the list of real-valued input variables [x 0 , . . ., x q ].The above expression is ill-defined if some of the inputs have repeated values, in which case one must resort to the use of limits.For instance, in the case where x 0 = x 1 = . . .= x q = x, the definition of divided differences reduces to: where f (n) (•) stands for the n-th derivative of f (•).Divided differences can alternatively be defined via the recursion relations with i ∈ {0, . . ., q − j}, j ∈ {1, . . ., q} and the initial conditions As is evident from the above, the divided differences of any given function f (•) with q inputs can be calculated with q(q − 1)/2 basic operations.Divided differences obey the following properties: (i) Linearity: for any two functions f (•) and g(•) and a number λ, the following holds.
In addition, a function of divided differences can be defined in terms of its Taylor expansion.For the special case of f (x) = e −itx , which will be of special importance in this study, we have Moreover, it is easy to verify that ) One may therefore write: (−it) q+n [x 0 , . . ., x q ] q+n (q + n)! .

III. DIVIDED DIFFERENCES FORMULATION OF THE TIME-EVOLUTION OPERATOR
Usually, the go-to approach to tackling the dynamics of time-dependent systems is through the use of (time-dependent) perturbation theory in the interaction-picture [13,14] where the Hamiltonian is written as namely as a sum of a static (or 'free') Hamiltonian H 0 (that is assumed to have a known spectrum) and an additional time-dependent Hamiltonian, V (t), that is treated as a perturbation to H 0 .The timeevolution operator in the interaction-picture is given by U I (t) = T exp[−i  (8) where H I (t) = e iH0t V (t)e −iH0t , and hereafter we use the q = 0 term to symbolize the identity operator.
The operator U I (t) evolves the interactionpicture wave-function |ψ I (t)⟩ which is related to the Schrödinger-picture wave-function via |ψ I (t)⟩ = e iH0t |ψ(t)⟩ (in our units, ℏ = 1).Similarly, the Schrödinger-picture time-evolution operator U (t) is related to the interaction-picture operator via U (t) = e −iH0t U I (t).In what follows we present an equivalent form for the Dyson series, Eq. ( 8), by systematically evaluating the integrals in the sum, writing V (t) as a sum of exponentials in t.

A. Generalized permutation operator representation of the perturbation Hamiltonian
We begin by denoting the eigenstates and eigenenergies of the free Hamiltonian H 0 by B = {|z⟩} and (For simplicity we assume a discrete countable set of eigenstates and eigenenergies).We will refer to B as the 'computational basis'.Next, we write the perturbation Hamiltonian V (t) as a sum of generalized permutation operators Π i [15]: where every generalized permutation operator is further expressed as a product of a (time dependent) diagonal (in the computational basis) operator D i and a bona-fide permutation operator P i .Specifically, the action of D i and P i on a computational basis states is given by D i |z⟩ = d i (z)|z⟩, where d i (z) is in general a complex number, and P i |z⟩ = |z ′ ⟩ for some |z ′ ⟩ ∈ B depending on i and z.The i = 0 permutation operator will be reserved to the identity operator, that is, P 0 = 1.Equipped with these notations, the action of a generalized permutation operator Π i on a basis state |z⟩ is given by where z ′ depends on both the state z and the operator index i.We note that any Hamiltonian can be readily cast in the above form [11].
At this point, we write each diagonal operator, D i (t), in Eq. ( 9) as an exponential sum in t, that is, where both Λ (k) i and D (k) i are (generally complexvalued) diagonal matrices and K i denotes the number of terms in the decomposition of D i .(For more details as to how to carry out this decomposition efficiently, see Refs.[16][17][18].)Thus, V (t) can be written as where, for simplicity, hereafter we fix K i = K ∀i, though this assumption can be easily removed.With this expansion of V (t), we can write the interactionpicture Hamiltonian H I (t) as where we have defined i .Using this form of H I (t) allows us to explicitly evaluate the time-ordered integrals of the Dyson series (8) (see Appendix A), arriving at where we have defined the (time-dependent) coefficients In addition, we have denoted i q = (i 1 , i 2 , . . ., i q ) and k q = (k 1 , k 2 , . . ., k q ) as multi-indices, P iq = ij |z ij ⟩ and |z ij ⟩ = P ij |z⟩ for every j = 0, . . ., q (remembering that the value of z ij depends on that of z).The inputs for the divided-difference exponential are given by for j = 0, . . ., q − 1, where ij |z ij ⟩ (we have omitted the dependence of x 0 , . . ., x q−1 on z, i q and on k q for a lighter notation).Importantly, the complex-valued e −it[x0,x1,...,xq−1,0] can be calculated efficiently, with computational complexity proportional to q 2 (see Ref. [19] for additional details).Equation ( 13) is the main technical result of our paper.The equation is a reformulation of the Dyson series, Eq. ( 8), yet it includes only sums of efficiently-computable terms and is devoid of both integrals and time-ordering operators.
The current representation of U I (t) also allows us to obtain an integral-free time-dependent perturbation expansion of the Schrödinger-picture timeevolution operator U (t): where with for j = 0, . . ., q.We arrive at Eq. ( 17) by using Identity 2 in Appendix B.
Equation ( 16) allows us to explicitly carry out time-dependent perturbation calculations directly in the Schrödinger-picture without the need to invoke the interaction-picture.

IV. USE CASES
To showcase the immediate applicability of our approach, in what follows we illustrate, by examining a few examples, the utility of Eq. ( 16) in a variety of settings, specifically in scenarios where the timeordered integrals of the Dyson series are cumbersome to calculate or lead to unnecessary complicated expressions.
Example 1: Transition amplitudes and Fermi's golden rule.Using the divided-differences series formulation of the time-evolution operator one may arrive at a rather simple expression for transition amplitudes, A(z in → z fin , t) = ⟨z fin |U (t)|z in ⟩, between an initial eigenstate |z in ⟩ and a final eigenstate |z fin ⟩ of the free Hamiltonian H 0 due to the perturbation effect of V (t), cf.Eq. ( 7).Using Eq. ( 16) we find where the second sum is over all multi-indices i q that obey P iq |z in ⟩ = |z fin ⟩.Consider, for example, for the canonical, yet non trivial, case of the non-Hermitian Hamiltonian H = H 0 + γF e iEt where γ is a perturbation parameter and F a permutation matrix such that F |z in ⟩ = |z fin ⟩ and F 2 = 1.While the standard Dyson-series-based derivation of the transition amplitude beyond the first order expansion q = 1 (that is Fermi's golden rule) is cumbersome to derive and to calculate (see, e.g., Refs.[20,21]), the present formulation allows us to immediately write a succinct expression for transition amplitudes that are accurate to any order in γ.In particular, the q-th order contribution for the transition amplitude is given very neatly by for odd values of q, and can be calculated with O(q 2 ) basic operations (for even values of q, A (q) (z in → z fin , t) = 0).The case q = 1 corresponds to Fermi's golden rule [22].This simple expression showcases how the new formulation of the Dyson series can be used to explore properties of time-varying Hamiltonian systems that arise from perturbation orders that have not explored with prior methods.
Example 2: Time-oscillating two-level Hamiltonian system.To further illustrate the practical power of the integral-free description of the Dyson series and its ability to assist in uncovering new physics, we next consider the dynamics of two-level systems driven by highly oscillating time-dependent Hamiltonians of the form H = ω 0 σ z + gσ x cos ωt.Despite its apparent simplicity, calculating dynamical properties of this system still poses challenges to numerical integration methods whenever the driving oscillations are very fast, and it is an open problem that continues to be a very active area of research, for example in the study of evolution of qubits undergoing quantum gates [23] or in Bloch-Siegert systems, which play an important role in atomic physics [24] (It should be noted that in certain cases two-level systems admit analytical solutions in which case a numerical series expansion is not strictly necessary, see e.g.[25]).Our formulation on the other hand results in a compact description of the unitary operator (more information is given in Appendix C, which in turn allows for a numerically-exact calculation of various properties of the system, such as transition amplitudes, even in regions of large ω where direct integration of the Schördinger equation becomes challenging.An example is given in Fig. 1 wherein a few numerically-exact results are given for various scenarios, as calculated by summing over the appropriate divided-difference terms.The figure also depicts independently derived numerically-exact solutions obtained using 4th-order Runge-Kutta simulations [26] carried out independently.As is evident, there is an excellent agreement between the two methods.
Example 3: Time-oscillating infinite-dimensional Hamiltonian system.Next, we work out the expansion of the time-evolution operator for a particle moving under the influence of a harmonic potential perturbed by a periodic time-dependent anharmonic term, namely, Identifying the static component as H 0 = p2 /2m + mω 2 x2 /2, we rewrite the Hamiltonian in terms of the harmonic annihilation-creation operators (recall that ℏ = 1) as To express V (t) in terms of generalized permutation operators, using the harmonic oscillator eigenstates {|n⟩} as the computational basis states, we rewrite the anharmonic operator â † + â 4 as â † + â 4 = i∈{0,±2,±4} D i P i where P i = n |n + i⟩⟨n| and the matrix elements of the diagonal operators D i are given by In terms of the above operators, the Hamiltonian can be written as where Having cast the Hamiltonian in the appropriate form, it is straightforward to calculate the β (iq) n terms in the divided-differences expansion of U (t), as per Eq. ( 17).For example for q = 1 we obtain × (D i1 ) nn (e −it[Ω−i1ω,0] + e −it[−(Ω+i1ω),0] ) where i 1 ∈ {0, ±2, ±4}.Our formulation provides an easy way to compute the state of the system after time t.For an initial state |ψ(0)⟩ = |n⟩, we get the analytic closed-form expression |ψ(t)⟩ = U (t)|n⟩ = iq β n P iq |n⟩ (see Appendix D for additional information and further analysis).

V. SUMMARY AND DISCUSSION
To conclude, in this work, we devised an alternative approach to time-dependent perturbation theory in quantum mechanics that allows one to read-ily calculate expansion terms.We derived an expansion that is equivalent to the usual Dyson series but which includes only sums of closed-form analytical simple-to-calculate expressions rather than the usual Dyson multi-dimensional time-ordered integrals.The terms at every expansion order in our new formulation coincide with those of the standard Dyson series, except that the Dyson integrals at every order replaced with finite sums.Therefore, both series share the same convergence criteria [27][28][29][30].However, our new formalism allowed us to write an integral-free perturbation expansion for the timeevolution operator.We illustrated the utility of our approach by working out a number of use cases and calculated the series coefficients for a number of examples for which the usual Dyson series calculation is cumbersome, demonstrating the functionality and practicality of our approach.
Another area in which the divided-differences expansion can be applied, which we have not explored here, is quantum algorithms -algorithms designed to be executed on quantum computers.Specifically, the divided-differences expansion was recently shown to be a valuable tool in the derivation of quantum algorithm devised to simulate the timeevolution of quantum states evolving under timeindependent [12] and time-dependent [31] Hamiltonians.There, it was shown that the divideddifferences expansion allows for the time-evolution operator to be written as a sum of generalized permutation matrices, equivalently a linear combination of unitary (LCU) operators.As such, this representation of the time-evolution operator lends itself naturally to the quantum LCU lemma [32] which provides a prescription for efficiently simulating such operators on quantum circuits.We leave that for future work.
We believe that the formulation introduced here will prove itself to be a powerful tool in the study of time-dependent quantum many-body systems.This result holds for arbitrarily close inputs and can be easily generalized to the case where inputs have repeated values, as claimed.

t 0 H
I (t)dt], a shorthand for the Dyson series