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

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 multi-dimensional 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 time-dependent 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–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.

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

$$\begin{equation}f[{x}_{0},\dots ,{x}_{q}]\equiv \sum\limits _{j=0}^{q}\frac{f({x}_{j})}{{\prod }_{k\ne j}({x}_{j}-{x}_{k})}\end{equation} \tag{ 1 }$$

with respect to the list of real-valued input variables [x₀, ..., 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₀ = x₁ = ... = x_q = x, the definition of divided differences reduces to:

$$\begin{equation}f[{x}_{0},\dots ,{x}_{q}]=\frac{{f}^{(q)}(x)}{q!},\end{equation} \tag{ 2 }$$

where f⁽ⁿ⁾(⋅) stands for the nth derivative of f(⋅). Divided differences can alternatively be defined via the recursion relations

$$\begin{equation}f[{x}_{i},\dots ,{x}_{i+j}]=\frac{f[{x}_{i+1},\dots ,{x}_{i+j}]-f[{x}_{i},\dots ,{x}_{i+j-1}]}{{x}_{i+j}-{x}_{i}},\end{equation} \tag{ 3 }$$

with i ∈ {0, ..., q − j}, j ∈ {1, ..., q} and the initial conditions

$$\begin{equation}f[{x}_{i}]=f({x}_{i}),\quad i\in \left\{0,\dots ,q\right\}\enspace \forall \enspace i.\end{equation} \tag{ 4 }$$

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. (f + g)[x₀, ..., x_q ] = f[x₀, ..., x_q ] + g[x₀, ..., x_q ] and (λ ⋅ f)[x₀, ..., x_q ] = λ ⋅ f[x₀, ..., x_q ]. (ii) Divided differences obey the Leibniz rule $(\enspace f\cdot g)[{x}_{0},\dots ,{x}_{q}]={\sum }_{k=0}^{q}\enspace f[{x}_{0},\dots ,{x}_{k}]\cdot g[{x}_{k},\dots ,{x}_{q}]$. (iii) Invariance under permutation of inputs, namely, f[x₀, ..., x_q ] = f[x_σ(0), ..., x_σ(q)] where σ: {0, ..., q} → {0, ..., q} is a permutation. (iv) Any real-valued function f(⋅) obeys the mean-value theorem f[x₀, ..., x_q ] = f^(q)(ξ)/q! where f^(q)(⋅) stands for the qth derivative of f(⋅) and ξ is a number in the interval [min[x₀, ..., x_q ], max[x₀, ..., x_q ]].

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

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}t[{x}_{0},\dots ,{x}_{q}]}=\sum\limits _{n=0}^{\infty }\frac{{(-it)}^{n}{[{x}_{0},\dots ,{x}_{q}]}^{n}}{n!}.\end{equation} \tag{ 5 }$$

Moreover, it is easy to verify that

$$\begin{equation}{[{x}_{0},\dots ,{x}_{q}]}^{q+n}=\begin{cases}0\qquad n< 0\quad \\ 1\qquad n=0\quad \\ {\sum }_{\sum {k}_{j}=n}{\prod }_{j=0}^{q}{x}_{j}^{{k}_{j}}\quad \qquad n > 0\quad \end{cases}.\end{equation} \tag{ 6 }$$

One may therefore write:

$$\begin{align}{\mathrm{e}}^{-\mathrm{i}t[{x}_{0},\dots ,{x}_{q}]}& =\sum\limits _{n=0}^{\infty }\frac{{(-it)}^{n}{[{x}_{0},\dots ,{x}_{q}]}^{n}}{n!}\\ & =\sum\limits _{n=q}^{\infty }\frac{{(-it)}^{n}{[{x}_{0},\dots ,{x}_{q}]}^{n}}{n!}\\ & =\sum\limits _{n=0}^{\infty }\frac{{(-it)}^{q+n}{[{x}_{0},\dots ,{x}_{q}]}^{q+n}}{(q+n)!}.\end{align}$$

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

$$\begin{equation}H(t)={H}_{0}+V(t),\end{equation} \tag{ 7 }$$

namely as a sum of a static (or 'free') Hamiltonian H₀ (that is assumed to have a known spectrum) and an additional time-dependent Hamiltonian, V(t), that is treated as a perturbation to H₀. The time-evolution operator in the interaction-picture is given by ${U}_{I}(t)=\mathcal{T}\enspace \mathrm{exp}[-\mathrm{i}{\int }_{0}^{t}{H}_{I}(t)\mathrm{d}t]$, a shorthand for the Dyson series

$$\begin{equation}{U}_{I}(t)=\sum\limits _{q=0}^{\infty }{(-i)}^{q} {\int }_{0}^{t} \mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}} \mathrm{d}{t}_{1}\enspace {H}_{I}({t}_{q})\dots {H}_{I}({t}_{1}),\end{equation} \tag{ 8 }$$

where ${H}_{I}(t)={\mathrm{e}}^{\mathrm{i}{H}_{0}t}V(t){\mathrm{e}}^{-\mathrm{i}{H}_{0}t}$, and hereafter we use the q = 0 term to symbolize the identity operator.

The operator U_I (t) evolves the interaction-picture wave-function |ψ_I (t)⟩ which is related to the Schrödinger-picture wave-function via $\vert {\psi }_{I}(t)\rangle ={\mathrm{e}}^{\mathrm{i}{H}_{0}t}\vert \psi (t)\rangle$ (in our units, ℏ = 1). Similarly, the Schrödinger-picture time-evolution operator U(t) is related to the interaction-picture operator via $U(t)={\mathrm{e}}^{-\mathrm{i}{H}_{0}t}{U}_{I}(t)$. In what follows we present an equivalent form for the Dyson series, equation (8), by systematically evaluating the integrals in the sum, writing V(t) as a sum of exponentials in t.

3.1. Generalized permutation operator representation of the perturbation Hamiltonian

We begin by denoting the eigenstates and eigenenergies of the free Hamiltonian H₀ by $\mathcal{B}=\left\{\vert z\rangle \right\}$ and $\mathcal{E}=\left\{{E}_{z}\right\}$, respectively, such that H₀|z⟩ = E_z |z⟩. (For simplicity we assume a discrete countable set of eigenstates and eigenenergies). We will refer to $\mathcal{B}$ as the 'computational basis'. Next, we write the perturbation Hamiltonian V(t) as a sum of generalized permutation operators Π_i [15]:

$$\begin{equation}V(t)=\sum\limits _{i=0}^{M}{{\Pi}}_{i}(t)=\sum\limits _{i=0}^{M}{D}_{i}(t){P}_{i},\end{equation} \tag{ 9 }$$

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 $\vert {z}^{\prime }\rangle \in \mathcal{B}$ depending on i and z. The i = 0 permutation operator will be reserved to the identity operator, that is, ${P}_{0}=\mathbb{1}$. Equipped with these notations, the action of a generalized permutation operator Π_i on a basis state |z⟩ is given by D_i P_i |z⟩ = d_i (z')|z'⟩, 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 equation (9) as an exponential sum in t, that is,

$$\begin{equation}{D}_{i}(t)=\sum\limits _{k=1}^{{K}_{i}}{\mathrm{e}}^{\mathrm{i}{{\Lambda}}_{i}^{(k)}t}{D}_{i}^{(k)},\end{equation} \tag{ 10 }$$

where both ${{\Lambda}}_{i}^{(k)}$ and ${D}_{i}^{(k)}$ are (generally complex-valued) 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 references [16–18].) Thus, V(t) can be written as

$$\begin{equation}V(t)=\sum\limits _{i=0}^{M}\sum\limits _{k=1}^{K}{\mathrm{e}}^{\mathrm{i}{{\Lambda}}_{i}^{(k)}t}{D}_{i}^{(k)}{P}_{i},\end{equation} \tag{ 11 }$$

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 interaction-picture Hamiltonian H_I (t) as

$$\begin{equation}{H}_{I}(t)=\sum\limits _{i=1}^{M}\sum\limits _{k=1}^{K}{\mathrm{e}}^{\mathrm{i}{H}_{0}^{(i,k)}t}{D}_{i}^{(k)}{P}_{i}\enspace {\mathrm{e}}^{-\mathrm{i}{H}_{0}t},\end{equation} \tag{ 12 }$$

where we have defined ${H}_{0}^{(i,k)}\equiv {H}_{0}+{{\Lambda}}_{i}^{(k)}$. 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

$$\begin{align}{U}_{I}(t)& =\sum\limits _{q=0}^{\infty }\sum\limits _{z}\sum\limits _{{\mathbf{i}}_{q},{\mathbf{k}}_{q}}\enspace {d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}\enspace {\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]}{P}_{{\mathbf{i}}_{q}}\vert z\rangle \langle z\vert \\ & =\sum\limits _{q=0}^{\infty }\sum\limits _{z}\sum\limits _{{\mathbf{i}}_{q}}\enspace {\alpha }_{z}^{({\mathbf{i}}_{q})}{P}_{{\mathbf{i}}_{q}}\vert z\rangle \langle z\vert ,\end{align} \tag{ 13 }$$

where we have defined the (time-dependent) coefficients

$$\begin{equation}{\alpha }_{z}^{({\mathbf{i}}_{q})}=\sum\limits _{{\mathbf{k}}_{q}}\enspace {d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}\enspace {\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]}.\end{equation} \tag{ 14 }$$

In addition, we have denoted i_q = (i₁, i₂, ..., i_q ) and k_q = (k₁, k₂, ..., k_q ) as multi-indices, ${P}_{{\mathbf{i}}_{q}}={P}_{{i}_{q}}\dots {P}_{{i}_{1}}$, ${d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}={\prod }_{j=1}^{q}{d}_{z}^{({\mathbf{i}}_{j},{k}_{j})}$ with ${d}_{z}^{({\mathbf{i}}_{j},{k}_{j})}=\langle {z}_{{\mathbf{i}}_{j}}\vert {D}_{{i}_{j}}^{({k}_{j})}\vert {z}_{{\mathbf{i}}_{j}}\rangle$ and $\vert {z}_{{\mathbf{i}}_{j}}\rangle ={P}_{{\mathbf{i}}_{j}}\vert z\rangle$ for every j = 0, ..., q (remembering that the value of ${z}_{{\mathbf{i}}_{j}}$ depends on that of z). The inputs for the divided-difference exponential are given by

$$\begin{equation}{x}_{j}={E}_{z}^{({\mathbf{i}}_{j})}-{E}_{z}^{({\mathbf{i}}_{q})}-\sum\limits _{\ell =j+1}^{q}{\lambda }_{z}^{({\mathbf{i}}_{\ell },{k}_{\ell })}\end{equation} \tag{ 15 }$$

for j = 0, ..., q − 1, where ${E}_{z}^{({\mathbf{i}}_{j})}=\langle {z}_{{\mathbf{i}}_{j}}\vert {H}_{0}\vert {z}_{{\mathbf{i}}_{j}}\rangle$ and ${\lambda }_{z}^{({\mathbf{i}}_{j},{k}_{j})}=\langle {z}_{{\mathbf{i}}_{j}}\vert {{\Lambda}}_{{i}_{j}}^{({k}_{j})}\vert {z}_{{\mathbf{i}}_{j}}\rangle$ (we have omitted the dependence of x₀, ..., x_q−1 on z, i_q and on k_q for a lighter notation). Importantly, the complex-valued ${\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]}$ can be calculated efficiently, with computational complexity proportional to q² (see reference [19] for additional details). Equation (13) is the main technical result of our paper. The equation is a reformulation of the Dyson series, equation (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 time-evolution operator U(t):

$$\begin{equation}U(t)={\text{e}}^{-\text{i}{H}_{0}t}\enspace {U}_{I}(t)=\sum\limits _{q=0}^{\infty }\sum\limits _{z}\sum\limits _{{\mathbf{i}}_{q}}\enspace {\beta }_{z}^{({\mathbf{i}}_{q})}{P}_{{\mathbf{i}}_{q}}\vert z\rangle \langle z\vert ,\end{equation} \tag{ 16 }$$

where

$$\begin{equation}{\beta }_{z}^{({\mathbf{i}}_{q})}={\alpha }_{z}^{({\mathbf{i}}_{q})}\enspace {\text{e}}^{-\text{i}t{E}_{z}^{({\mathbf{i}}_{q})}} = \sum\limits _{{\mathbf{k}}_{q}}\enspace {d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}\enspace {\mathrm{e}}^{-\mathrm{i}t[{y}_{0},{y}_{1},\dots ,{y}_{q}]},\end{equation} \tag{ 17 }$$

with ${y}_{j}={E}_{z}^{({\mathbf{i}}_{j})}-{\sum }_{\ell =j+1}^{q}{\lambda }_{z}^{({\mathbf{i}}_{\ell },{k}_{\ell })}$ for j = 0, ..., q. We arrive at equation (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.

To showcase the immediate applicability of our approach, in what follows we illustrate, by examining a few examples, the utility of equation (16) in a variety of settings, specifically in scenarios where the time-ordered 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₀ due to the perturbation effect of V(t), cf equation (7). Using equation (16) we find

$$\begin{equation}A({z}_{\text{in}}\to {z}_{\text{fin}},t)=\sum\limits _{q=0}^{\infty } \sum\limits _{\enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace {\mathbf{i}}_{q}:\langle {z}_{\text{fin}}\vert {P}_{{\mathbf{i}}_{q}}\vert {z}_{\text{in}}\rangle =1} {\beta }_{{z}_{\text{in}}}^{({\mathbf{i}}_{q})},\end{equation} \tag{ 18 }$$

where the second sum is over all multi-indices i_q that obey ${P}_{{\mathbf{i}}_{q}}\vert {z}_{\text{in}}\rangle =\vert {z}_{\text{fin}}\rangle$.

Consider, for example, for the canonical, yet non trivial, case of the non-Hermitian Hamiltonian H = H₀ + γF e^iEt where γ is a perturbation parameter and F a permutation matrix such that F|z_in⟩ = |z_fin⟩ and ${F}^{2}=\mathbb{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. references [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 qth order contribution for the transition amplitude is given very neatly by

$$\begin{equation}{A}^{(q)}({z}_{\text{in}}\to {z}_{\text{fin}},t)={\gamma }^{q}\enspace {\mathrm{e}}^{-\mathrm{i}t{E}_{{z}_{\text{in}}}}\enspace {\mathrm{e}}^{-\mathrm{i}t[qE,(q-1)E,\dots ,E,0]},\end{equation} \tag{ 19 }$$

for odd values of q, and can be calculated with O(q²) 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 = ω₀ σ_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 figure 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.

Zoom In Zoom Out Reset image size

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,

$$\begin{equation}H(t)=\frac{{\hat{p}}^{2}}{2m}+\frac{1}{2}m{\omega }^{2}{\hat{x}}^{2}+\frac{1}{4}{\Gamma}\enspace \mathrm{cos}({\Omega}t){\hat{x}}^{4}.\end{equation} \tag{ 20 }$$

Identifying the static component as ${H}_{0}={\hat{p}}^{2}/2m+m{\omega }^{2}{\hat{x}}^{2}/2$, we rewrite the Hamiltonian in terms of the harmonic annihilation–creation operators (recall that ℏ = 1) as

$$\begin{equation}H(t)=\omega \left({\hat{a}}^{{\dagger}}\hat{a}+\frac{1}{2}\right)+\underset{V(t)}{\underbrace{\frac{1}{4}{\Gamma}\enspace \mathrm{cos}({\Omega}t){\left(\frac{1}{2m\omega }\right)}^{2}{({\hat{a}}^{{\dagger}}+\hat{a})}^{4}}}.\end{equation} \tag{ 21 }$$

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 ${\left({\hat{a}}^{{\dagger}}+\hat{a}\right)}^{4}$ as ${\left({\hat{a}}^{{\dagger}}+\hat{a}\right)}^{4}={\sum }_{i\in \left\{0,\pm 2,\pm 4\right\}}{D}_{i}{P}_{i}$ where P_i = ∑_n |n + i⟩⟨n| and the matrix elements of the diagonal operators D_i are given by

$$\begin{align}{({D}_{4})}_{nn}& =\sqrt{n(n-1)(n-2)(n-3)}\enspace ,\\ {({D}_{2})}_{nn}& =\sqrt{n(n-1)}(4n-2)\enspace ,\\ {({D}_{0})}_{nn}& =3(2{n}^{2}+2n+1)\enspace ,\\ {({D}_{-2})}_{nn}& =\sqrt{(n+1)(n+2)}(4n+6)\enspace ,\\ {({D}_{-4})}_{nn}& =\sqrt{(n+1)(n+2)(n+3)(n+4)}\enspace .\\ & \phantom{\rule{0ex}{0ex}}\end{align} \tag{ 22 }$$

In terms of the above operators, the Hamiltonian can be written as

$$\begin{equation}H(t)={H}_{0}+\gamma \sum\limits _{i}\sum\limits _{k=\pm 1}{\mathrm{e}}^{-\mathrm{i}k{\Omega}t}{D}_{i}{P}_{i},\end{equation} \tag{ 23 }$$

where ${H}_{0}=\omega {\sum }_{n}(n+\frac{1}{2})\vert n\rangle \langle n\vert$ and $\gamma =\frac{1}{8}{\Gamma}{\left(\frac{1}{2m\omega }\right)}^{2}$. Having cast the Hamiltonian in the appropriate form, it is straightforward to calculate the ${\beta }_{n}^{({\mathbf{i}}_{q})}$ terms in the divided-differences expansion of U(t), as per equation (17). For example for q = 1 we obtain

$$\begin{align}{\beta }_{n}^{({i}_{1})}& =\gamma \enspace {\mathrm{e}}^{-\mathrm{i}t{E}_{n+{i}_{1}}}{({D}_{{i}_{1}})}_{nn}({\mathrm{e}}^{-\mathrm{i}t[{\Omega}-{i}_{1}\omega ,0]}+{\mathrm{e}}^{-\mathrm{i}t[-({\Omega}+{i}_{1}\omega ),0]}),\end{align} \tag{ 24 }$$

where i₁ ∈ {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 $\vert \psi (t)\rangle =U(t)\vert n\rangle ={\sum }_{{\mathbf{i}}_{q}}{\beta }_{n}^{({\mathbf{i}}_{q})}{P}_{{\mathbf{i}}_{q}}\vert n\rangle$ (see appendix D for additional information and further analysis).

To conclude, in this work, we devised an alternative approach to time-dependent perturbation theory in quantum mechanics that allows one to readily 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–30]. However, our new formalism allowed us to write an integral-free perturbation expansion for the time-evolution 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 time-evolution of quantum states evolving under time-independent [12] and time-dependent [31] Hamiltonians. There, it was shown that the divided-differences 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 work is supported by the US Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES) under Award No. DE-SC0020280.

The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.

Here, we derive the expression

$$\begin{equation}{U}_{I}(t)=\sum\limits _{q=0}^{\infty }\sum\limits _{z}\sum\limits _{{\mathbf{i}}_{q},{\mathbf{k}}_{q}}\enspace {d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}\enspace {\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]}{P}_{{\mathbf{i}}_{q}}\vert z\rangle \langle z\vert .\end{equation} \tag{ A1 }$$

Using the form of H_I (t)

$$\begin{equation}{H}_{I}(t)=\sum\limits _{i=1}^{M}\sum\limits _{k=1}^{K}{\mathrm{e}}^{\mathrm{i}{H}_{0}^{(i,k)}t}{D}_{i}^{(k)}{P}_{i}\enspace {\mathrm{e}}^{-\mathrm{i}{H}_{0}t},\end{equation} \tag{ A2 }$$

we can write the Dyson series for U_I (t) as

$$\begin{align}{U}_{I}(t)& =\sum\limits _{q=0}^{\infty }{(-i)}^{q}{\int }_{0}^{t}\mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}}\mathrm{d}{t}_{1}\left(\sum\limits _{{i}_{q}=1}^{M}\sum\limits _{{k}_{q}=1}^{K}{\mathrm{e}}^{\mathrm{i}{H}_{0}^{({i}_{q},{k}_{q})}{t}_{q}}{D}_{{i}_{q}}^{({k}_{q})}{P}_{{i}_{q}}{\mathrm{e}}^{-i{H}_{0}{t}_{q}}\right)\dots \\ & \quad \times \dots \left(\sum\limits _{{i}_{1}=1}^{M}\sum\limits _{{k}_{1}=1}^{K}{\mathrm{e}}^{\mathrm{i}{H}_{0}^{({i}_{1},{k}_{1})}{t}_{1}}{D}_{{i}_{1}}^{({k}_{1})}{P}_{{i}_{1}}{\mathrm{e}}^{-\mathrm{i}{H}_{0}{t}_{1}}\right)\\ & =\sum\limits _{q=0}^{\infty }\sum\limits _{{\mathbf{i}}_{q}}\sum\limits _{{\mathbf{k}}_{q}}{(-i)}^{q}{\int }_{0}^{t}\mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}}\mathrm{d}{t}_{1}\left({\mathrm{e}}^{\mathrm{i}{H}_{0}^{({i}_{q},{k}_{q})}{t}_{q}}{D}_{{i}_{q}}^{({k}_{q})}{P}_{{i}_{q}}\enspace {\mathrm{e}}^{-\mathrm{i}{H}_{0}{t}_{q}}\dots {\mathrm{e}}^{\mathrm{i}{H}_{0}^{({i}_{1},{k}_{1})}{t}_{1}}{D}_{{i}_{1}}^{({k}_{1})}{P}_{{i}_{1}}\enspace {\mathrm{e}}^{-\mathrm{i}{H}_{0}{t}_{1}}\right),\end{align} \tag{ A3 }$$

where i_q = (i₁, i₂, ..., i_q ) and k_q = (k₁, k₂, ..., k_q ) are multi-indices. Acting with U_I (t) on an arbitrary computational basis state |z⟩, we get

$$\begin{align}{U}_{I}(t)\vert z\rangle & =\sum\limits _{q=0}^{\infty }\sum\limits _{{\mathbf{i}}_{q}}\sum\limits _{{\mathbf{k}}_{q}}{(-i)}^{q}{\int }_{0}^{t}\mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}}\mathrm{d}{t}_{1}\enspace {\mathrm{e}}^{\mathrm{i}{E}_{z}^{({\mathbf{i}}_{q},{k}_{q})}{t}_{q}}{d}_{z}^{({\mathbf{i}}_{q},{k}_{q})}\enspace {\mathrm{e}}^{-\mathrm{i}{E}_{z}^{({\mathbf{i}}_{q-1})}{t}_{q}}\dots {\mathrm{e}}^{\mathrm{i}{E}_{z}^{({\mathbf{i}}_{1},{k}_{1})}{t}_{1}}{d}_{z}^{({\mathbf{i}}_{1},{k}_{1})}\enspace {\mathrm{e}}^{-\mathrm{i}{E}_{z}{t}_{1}}{P}_{{\mathbf{i}}_{q}}\vert z\rangle \\ & =\sum\limits _{q=0}^{\infty }\sum\limits _{{\mathbf{i}}_{q}}\sum\limits _{{\mathbf{k}}_{q}}{d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}\left({(-i)}^{q}{\int }_{0}^{t}\mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}}\mathrm{d}{t}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}({E}_{z}^{({\mathbf{i}}_{q-1})}-{E}_{z}^{({\mathbf{i}}_{q},{k}_{q})}){t}_{q}}\dots {\mathrm{e}}^{-\mathrm{i}({E}_{z}-{E}_{z}^{({\mathbf{i}}_{1},{k}_{1})}){t}_{1}}\right){P}_{{\mathbf{i}}_{q}}\vert z\rangle ,\end{align} \tag{ A4 }$$

where we have defined ${P}_{{\mathbf{i}}_{q}}={P}_{{i}_{q}}\cdots {P}_{{i}_{1}}$ and $\vert {z}_{{\mathbf{i}}_{j}}\rangle ={P}_{{\mathbf{i}}_{j}}\vert z\rangle$ for every j = 0, ..., q (remembering that ${z}_{{\mathbf{i}}_{j}}$ depends on z). Moreover, we define ${E}_{z}^{({\mathbf{i}}_{j})}=\langle {z}_{{\mathbf{i}}_{j}}\vert {H}_{0}\vert {z}_{{\mathbf{i}}_{j}}\rangle$ and similarly

$$\begin{equation}{E}_{z}^{({\mathbf{i}}_{j},{k}_{j})}=\langle {z}_{{\mathbf{i}}_{j}}\vert {H}_{0}^{({i}_{j},{k}_{j})}\vert {z}_{{\mathbf{i}}_{j}}\rangle ={E}_{z}^{({\mathbf{i}}_{j})}+{\lambda }_{z}^{({\mathbf{i}}_{j},{k}_{j})},\end{equation} \tag{ A5 }$$

with ${\lambda }_{z}^{({\mathbf{i}}_{j},{k}_{j})}=\langle {z}_{{\mathbf{i}}_{j}}\vert {{\Lambda}}_{{i}_{j}}^{({k}_{j})}\vert {z}_{{\mathbf{i}}_{j}}\rangle$. In addition, we have denoted ${d}_{z}^{({\mathbf{i}}_{j},{k}_{j})}=\langle {z}_{{\mathbf{i}}_{j}}\vert {D}_{{i}_{j}}^{({k}_{j})}\vert {z}_{{\mathbf{i}}_{j}}\rangle$ and ${d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}={\prod }_{j=1}^{q}{d}_{z}^{({\mathbf{i}}_{j},{k}_{j})}$. Following identity 1 in appendix B, the integrals appearing inside the parentheses in equation (A4) evaluate explicitly to the divided-differences exponential, that is,

$$\begin{equation}{(-i)}^{q}{\int }_{0}^{t}\mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}}\mathrm{d}{t}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}({E}_{z}^{({\mathbf{i}}_{q-1})}-{E}_{z}^{({\mathbf{i}}_{q},{k}_{q})}){t}_{q}}\dots {\mathrm{e}}^{-\mathrm{i}({E}_{z}-{E}_{z}^{({\mathbf{i}}_{1},{k}_{1})}){t}_{1}}={\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]},\end{equation} \tag{ A6 }$$

with

$$\begin{equation}{x}_{j}=\sum\limits _{\ell =j+1}^{q}{E}_{z}^{({\mathbf{i}}_{\ell -1})}-{E}_{z}^{({\mathbf{i}}_{\ell },{k}_{\ell })}={E}_{z}^{({\mathbf{i}}_{j})}-{E}_{z}^{({\mathbf{i}}_{q})}-\sum\limits _{\ell =j+1}^{q}{\lambda }_{z}^{({\mathbf{i}}_{\ell },{k}_{\ell })}\quad \text{for}\enspace j=0,\dots ,q-1\end{equation} \tag{ A7 }$$

(we have omitted the dependence of x₀, ..., x_q−1 on z, i_q and on k_q for a lighter notation). Using equation (A6), we can write the action of U_I (t) on an arbitrary computational basis state in an integral-free form as

$$\begin{equation}{U}_{I}(t)\vert z\rangle =\sum\limits _{q=0}^{\infty }\sum\limits _{{\mathbf{i}}_{q}}\sum\limits _{{\mathbf{k}}_{q}}{d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}\enspace {\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]}{P}_{{\mathbf{i}}_{q}}\vert z\rangle ,\end{equation} \tag{ A8 }$$

with the inputs to the divided-differences exponential given in equation (A7). We thus arrive at:

$$\begin{equation}{U}_{I}(t)=\sum\limits _{q=0}^{\infty }\sum\limits _{z}\sum\limits _{{\mathbf{i}}_{q}}\sum\limits _{{\mathbf{k}}_{q}}{d}_{z}^{({\mathbf{i}}_{q},{\mathbf{k}}_{q})}\enspace {\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]}{P}_{{\mathbf{i}}_{q}}\vert z\rangle \langle z\vert ,\end{equation} \tag{ A9 }$$

as claimed.

Identity 1. 

$$\begin{equation}{(-i)}^{q}{\int }_{0}^{t}\mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}}\mathrm{d}{t}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}({\gamma }_{q}{t}_{q}+\cdots +{\gamma }_{1}{t}_{1})}={\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]},\end{equation} \tag{ B1 }$$

where ${x}_{j}={\sum }_{k=j+1}^{q}{\gamma }_{k}$. This identity is a variant of what is often known as the Hermite–Gennochi formula [8].

Proof. We start with the left-hand side of equation (B1). Making a change of variables t_j = ts_j we get

$$\begin{equation}{(-i)}^{q}{\int }_{0}^{t}\mathrm{d}{t}_{q}\dots {\int }_{0}^{{t}_{2}}\mathrm{d}{t}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}({\gamma }_{q}{t}_{q}+\cdots +{\gamma }_{1}{t}_{1})}={(-it)}^{q}{\int }_{0}^{1}\mathrm{d}{s}_{q}\dots {\int }_{0}^{{s}_{2}}\mathrm{d}{s}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}t({\gamma }_{q}{s}_{q}+\cdots +{\gamma }_{1}{s}_{1})}.\end{equation} \tag{ B2 }$$

Next we prove by induction that

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q-1},0]}={(-it)}^{q}{\int }_{0}^{1}\mathrm{d}{s}_{q}\dots {\int }_{0}^{{s}_{2}}\mathrm{d}{s}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}t({\gamma }_{q}{s}_{q}+\cdots +{\gamma }_{1}{s}_{1})},\end{equation} \tag{ B3 }$$

where ${x}_{j}={\sum }_{k=j+1}^{q}{\gamma }_{k}$.

Base step (proving for q = 1): starting from the left-hand-side of equation (B3) we have

$$\begin{equation}(-it){\int }_{0}^{1}\mathrm{d}{s}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}t{\gamma }_{1}{s}_{1}}=\frac{(-it)}{{\gamma }_{1}}{\int }_{0}^{{\gamma }_{1}}\mathrm{d}\xi \enspace {\mathrm{e}}^{-\mathrm{i}t\xi }=\frac{1}{{\gamma }_{1}}{\int }_{0}^{{\gamma }_{1}}\mathrm{d}\xi \frac{\mathrm{d}}{\mathrm{d}\xi }\enspace {\mathrm{e}}^{-\mathrm{i}t\xi }=\frac{{\mathrm{e}}^{-\mathrm{i}t{\gamma }_{1}}-1}{{\gamma }_{1}}={\mathrm{e}}^{-\mathrm{i}t[0,{\gamma }_{1}]}={\mathrm{e}}^{-\mathrm{i}t[{\gamma }_{1},0]}\end{equation} \tag{ B4 }$$

as required. In the first equality of equation (B4), we changed the integration variable ξ = s₁ γ₁, so that ξ ranges from 0 to γ₁, and γ₁ ds₁ = dξ.

Hypothesis step: next, we assume the validity of equation (B3).

Induction step: based on this assumption, we now prove that

$$\begin{equation}{(-it)}^{q+1}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{0}^{{s}_{2}}\mathrm{d}{s}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}t({\gamma }_{q+1}{s}_{q+1}+\cdots +{\gamma }_{1}{s}_{1})}={\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},\dots ,{x}_{q},0]}\end{equation} \tag{ B5 }$$

with ${x}_{j}={\sum }_{k=j+1}^{q+1}{\gamma }_{k}$. Changing the integration variable s₁ to ξ = γ_q+1 s_q+1 + ⋯ + γ₁ s₁, so that ξ ranges from γ_q+1 s_q+1 + ⋯ + γ₂ s₂ to γ_q+1 s_q+1 + ⋯ + (γ₂ + γ₁)s₂, and γ₁ ds₁ = dξ we get

$$\begin{align}& {(-it)}^{q+1}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{0}^{{s}_{2}}\mathrm{d}{s}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}t\left({\gamma }_{q+1}{s}_{q+1}+\cdots +{\gamma }_{1}{s}_{1}\right)}\\ & \qquad ={(-it)}^{q+1}\frac{1}{{\gamma }_{1}}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{{\gamma }_{q+1}{s}_{q+1}+\cdots +{\gamma }_{2}{s}_{2}}^{{\gamma }_{q+1}{s}_{q+1}+\cdots +({\gamma }_{2}+{\gamma }_{1}){s}_{2}}\mathrm{d}\xi \enspace {\mathrm{e}}^{-\mathrm{i}t\xi }\\ & \qquad ={(-it)}^{q}\frac{1}{{\gamma }_{1}}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{{\gamma }_{q+1}{s}_{q+1}+\cdots +{\gamma }_{2}{s}_{2}}^{{\gamma }_{q+1}{s}_{q+1}+\cdots +({\gamma }_{2}+{\gamma }_{1}){s}_{2}}\mathrm{d}\xi \frac{\mathrm{d}}{\mathrm{d}\xi }\enspace {\mathrm{e}}^{-\mathrm{i}t\xi }\\ & \qquad ={(-it)}^{q}\frac{1}{{\gamma }_{1}}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{0}^{{s}_{3}}\mathrm{d}{s}_{2}\left({\mathrm{e}}^{-\mathrm{i}t\left({\gamma }_{q+1}{s}_{q+1}+\cdots +({\gamma }_{2}+{\gamma }_{1}){s}_{2}\right)}-{\mathrm{e}}^{-\mathrm{i}t\left({\gamma }_{q+1}{s}_{q+1}+\cdots +{\gamma }_{2}{s}_{2}\right)}\right)\\ & \qquad =\frac{1}{{\gamma }_{1}}\left(\underset{({\ast})}{\underbrace{{(-it)}^{q}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{0}^{{s}_{3}}\mathrm{d}{s}_{2}\enspace {\mathrm{e}}^{-\mathrm{i}t\left({\gamma }_{q+1}{s}_{q+1}+\cdots +({\gamma }_{2}+{\gamma }_{1}){s}_{2}\right)}}}\right.\\ & \left.\quad \qquad \enspace \enspace \enspace -\underset{({\ast}{\ast})}{\underbrace{{(-it)}^{q}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{0}^{{s}_{3}}\mathrm{d}{s}_{2}\enspace {\mathrm{e}}^{-\mathrm{i}t\left({\gamma }_{q+1}{s}_{q+1}+\cdots +{\gamma }_{2}{s}_{2}\right)}}}\right),\end{align} \tag{ B6 }$$

where from the induction assumption,

$$\begin{equation}({\ast})={\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{2},\dots ,{x}_{q},0]},\quad \text{and}\quad ({\ast}{\ast})={\mathrm{e}}^{-\mathrm{i}t[{x}_{1},{x}_{2},\dots ,{x}_{q},0]},\end{equation} \tag{ B7 }$$

with ${x}_{j}={\sum }_{k=j+1}^{q+1}{\gamma }_{k}$. Therefore, we obtain,

$$\begin{align}& {(-it)}^{q+1}{\int }_{0}^{1}\mathrm{d}{s}_{q+1}\dots {\int }_{0}^{{s}_{2}}\mathrm{d}{s}_{1}\enspace {\mathrm{e}}^{-\mathrm{i}t\left({\gamma }_{q+1}{s}_{q+1}+\cdots +{\gamma }_{1}{s}_{1}\right)}\\ & \qquad =\frac{{\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{2},\dots ,{x}_{q},0]}-{\mathrm{e}}^{-\mathrm{i}t[{x}_{1},{x}_{2},\dots ,{x}_{q},0]}}{{\gamma }_{1}}=\frac{{\mathrm{e}}^{-\mathrm{i}t[{x}_{2},\dots ,{x}_{q},0,{x}_{0}]}-{\mathrm{e}}^{-\mathrm{i}t[{x}_{2},\dots ,{x}_{q},0,{x}_{1}]}}{{x}_{0}-{x}_{1}}\\ & \qquad ={\mathrm{e}}^{-\mathrm{i}t[{x}_{2},\dots ,{x}_{q},0,{x}_{0},{x}_{1}]}={\mathrm{e}}^{-\mathrm{i}t[{x}_{0},{x}_{1},{x}_{2},\dots ,{x}_{q},0]}.\end{align} \tag{ B8 }$$

This concludes the proof.

Identity 2. Given an arbitrary multi-set of inputs {x₀, ..., x_q },

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}t[{x}_{0},\dots ,{x}_{q}]}={\mathrm{e}}^{-\mathrm{i}tx}\enspace {\mathrm{e}}^{-\mathrm{i}t[{{\Delta}}_{0},\dots ,{{\Delta}}_{q}]},\end{equation} \tag{ B9 }$$

where x is an arbitrary constant and Δ_j = x_j − x.

Proof. We prove equation (B9). By definition [8],

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}t[{x}_{0},\dots ,{x}_{q}]}=\sum\limits _{j}\frac{{\mathrm{e}}^{-\mathrm{i}t{x}_{j}}}{{\prod }_{k\ne j}({x}_{j}-{x}_{k})}\end{equation} \tag{ B10 }$$

(assuming for now that all inputs are distinct). It follows then that

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}t[{x}_{0},\dots ,{x}_{q}]}=\sum\limits _{j}\frac{{\mathrm{e}}^{-\mathrm{i}t({{\Delta}}_{j}+x)}}{{\prod }_{k\ne j}({{\Delta}}_{j}-{{\Delta}}_{k})}={\mathrm{e}}^{-\mathrm{i}tx}\sum\limits _{j}\frac{{\mathrm{e}}^{-\mathrm{i}t{{\Delta}}_{j}}}{{\prod }_{k\ne j}({{\Delta}}_{j}-{{\Delta}}_{k})}={\mathrm{e}}^{-\mathrm{i}tx}\enspace {\mathrm{e}}^{-\mathrm{i}t[{{\Delta}}_{0},\dots ,{{\Delta}}_{q}]}.\end{equation} \tag{ B11 }$$

This result holds for arbitrarily close inputs and can be easily generalized to the case where inputs have repeated values, as claimed.

Consider a time-dependent qubit system whose Hamiltonian has the form H = ω₀ σ_z + gσ_x  cos ωt. We identify H₀ = ω₀ σ_z (hence the computational basis is the Pauli-Z eigenbasis {|z⟩: z = 0, 1}) and V(t) = gσ_x  cos ωt. The time-dependent component V(t) contains a one permutation operators (M = 1), P₁ = σ_x , whose associated diagonal operators, ${D}_{1}=g\enspace \mathrm{cos}\enspace \omega t\enspace \mathbb{1}=(g/2)({\mathrm{e}}^{\mathrm{i}\omega t}+{\mathrm{e}}^{-\mathrm{i}\omega t})\mathbb{1}$ has two exponential terms (K = 2) such that ${D}_{1}={\mathrm{e}}^{\mathrm{i}\omega t}{D}_{1}^{(1)}+{\mathrm{e}}^{-\mathrm{i}\omega t}{D}_{1}^{(-1)}$ with ${D}_{1}^{(1)}={D}_{1}^{(-1)}=(g/2)\mathbb{1}$. For reasons that would become clearer later we index the latter operators by k = ±1 (instead of k = 1, 2).

To write the time evolution operator, we first note that since there is only one off-diagonal operator in the Hamiltonian, there is also one sequence of off-diagonal operators per expansion order q, explicitly:

$$\begin{equation}{P}_{{\mathbf{i}}_{q}}={\sigma }_{x}^{q}=\begin{cases}{\sigma }_{x},\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{o}\mathrm{d}\mathrm{d}\enspace q\quad \\ \mathbb{1},\quad \enspace \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\enspace q.\quad \end{cases}\end{equation} \tag{ C1 }$$

Moreover ${\beta }_{z}^{(q)}$, defined in equation (12) in the main text, is given by

$$\begin{equation}{\beta }_{z}^{(q)}={\left(\frac{g}{2}\right)}^{q}\sum\limits _{{\mathbf{k}}_{q}}{\mathrm{e}}^{-\mathrm{i}t[{y}_{0},\dots ,{y}_{q}]},\end{equation} \tag{ C2 }$$

with ${y}_{j}={E}_{z}^{(j)}-\omega {\Sigma}{\mathbf{k}}_{j+1:q}$ for j = 0, ..., q, where ${E}_{z}^{(j)}=\langle z\vert {\sigma }_{x}^{j}{D}_{0}{\sigma }_{x}^{j}\vert z\rangle ={(-1)}^{z+j}{\omega }_{0}$ and Σk_j+1:q is a shorthand notation to ${\sum }_{\ell =j+1}^{q}{k}_{\ell }$. The last two equations define the time evolution operator of the system for any order expansion q, specifically,

$$\begin{equation}U(t)=\sum\limits _{z=0,1}\left(\sum\limits _{q=0}^{\infty }{\beta }_{z}^{(2q)}(t)\right)\vert z\rangle \langle z\vert +\sum\limits _{z=0,1}\left(\sum\limits _{q=0}^{\infty }{\beta }_{z}^{(2q+1)}(t)\right)\vert \bar{z}\rangle \langle z\vert \end{equation} \tag{ C3 }$$

where $\vert \bar{z}\rangle =\vert 1-z\rangle$.

Below we provide additional calculations for the time-dependent anharmonic oscillator

$$\begin{equation}H(t)=\omega \left({\hat{a}}^{{\dagger}}\hat{a}+\frac{1}{2}\right)+\frac{1}{4}{\Gamma}\enspace \mathrm{cos}({\Omega}t){\left(\frac{1}{2m\omega }\right)}^{2}{({\hat{a}}^{{\dagger}}+\hat{a})}^{4}\end{equation} \tag{ D1 }$$

discussed in the main text. Having computed the coefficients, ${\beta }_{n}^{({\mathbf{i}}_{q})}$ (see equation (19) in the main text), we may plot, as we do in figure 2(a), the population of every mode at various times. In the figure, the populations are given at t = 0.04 for ω = 1, Ω = 2, γ = 0.02 with the initial state |n = 4⟩ and an expansion cutoff of Q = 5. To ascertain the accuracy of the Dyson series truncated at different cutoff orders Q, we may also contrast the divided-differences expansion with exact-numerical results obtained via direct integration of the Schrödinger equation. In figure 2(b) we plot the infidelity 1 − |⟨ψ(t)|ψ_Q (t)⟩|² between the state |ψ(t)⟩, the solution of the time-dependent Schrödinger equation of this system, and the state as obtained by using divided-differences expansion for U(t), |ψ_Q (t)⟩ = U_Q (t)|4⟩, that is, the state obtained by evolving |4⟩ under U(t) given by

$$\begin{equation}U(t)={\text{e}}^{-\text{i}{H}_{0}t}\enspace {U}_{I}(t)=\sum\limits _{q=0}^{\infty }\sum\limits _{z}\sum\limits _{{\mathbf{i}}_{q}}{\beta }_{z}^{({\mathbf{i}}_{q})}{P}_{{\mathbf{i}}_{q}}\vert z\rangle \langle z\vert ,\end{equation} \tag{ D2 }$$

as a function of evolution time t.

Zoom In Zoom Out Reset image size