Exponential Perturbative Expansions and Coordinate Transformations

Abstract

We propose a unified approach for different exponential perturbation techniques used in the treatment of time-dependent quantum mechanical problems, namely the Magnus expansion, the Floquet–Magnus expansion for periodic systems, the quantum averaging technique, and the Lie–Deprit perturbative algorithms. Even the standard perturbation theory fits in this framework. The approach is based on carrying out an appropriate change of coordinates (or picture) in each case, and it can be formulated for any time-dependent linear system of ordinary differential equations. All of the procedures (except the standard perturbation theory) lead to approximate solutions preserving by construction unitarity when applied to the time-dependent Schrödinger equation.

1. Introduction

Linear differential equations of the form

$$\frac{dU}{dt}=A\left(t\right)U,U\left(0\right)=I$$

(1)

are ubiquitous in many branches of physics, chemistry and mathematics. Here, U is a real or complex $d\times d$ matrix, and $A\left(t\right)$ is a sufficiently smooth matrix to ensure the existence of solutions. Perhaps the most important example corresponds to the Schrödinger equation for the evolution operator in quantum systems with a time-dependent Hamiltonian $H\left(t\right)$, in which case $A\left(t\right)=-iH\left(t\right)/\hslash$. Particular cases include spin dynamics in magnetic resonance (Nuclear Magnetic Resonance—NMR, Electronic Paramagnetic Resonance—EPR, Dynamic Nuclear Polarization—DNP, etc.) [1,2,3], electron-atom collisions in atomic physics, pressure broadening of rotational spectra in molecular physics, control of chemical reactions with driving induced by laser beams, etc. When the time-dependence of the Hamiltonian is periodic, as occurs, for instance, in periodically driven quantum systems, atomic quantum gases in periodically driven optical lattices, etc. [4,5], the Floquet theorem [6] relates $H\left(t\right)$ with a constant Hamiltonian. More specifically, it implies that the evolution operator is factorized as $U\left(t\right)=P\left(t\right)\exp \left(tF\right)$, with $P\left(t\right)$ a periodic time-dependent matrix and F a constant matrix. This theorem has been widely used in problems of solid state physics, in NMR, in the quantum simulation of systems with time-independent Hamiltonians by periodically driven quantum systems, etc. [7]. The Average Hamiltonian Theory is also closely related with this result, and the effective Hamiltonian is an important tool in the description of the system [8,9].

In general, Equation (1) cannot be solved in closed form, and so different approaches have been proposed along the years to obtain approximations, both analytic and numerical. Among the former, we can mention the standard perturbation theory, the average Hamiltonian theory, and the Magnus expansion [5,8,10]. Concerning the second approach, different numerical algorithms have been applied to obtain solutions on specific time intervals [11,12,13].

In this work, we will concentrate on different techniques providing analytic approximations to the solution of (1) that also share some of its most salient qualitative features. In particular, if (1) represents the Schrödinger equation, it is well known that $U\left(t\right)$ is unitary for all t, and this guarantees that its elements represent probabilities of transition between the different states of the system. However, it happens that not every approximate scheme (either analytic or numerical) renders unitary matrices and, thus, the physical description they provide is unreliable, especially for large integration times.

The Magnus expansion [14] presents the remarkable feature that it allows one to express the solution as the exponential of a series, $U\left(t\right)=\exp \left(\mathsf{\Omega }\right(t\left)\right)$, so that, even when the series $\mathsf{\Omega }\left(t\right)={\sum }_{n\ge 1}{\mathsf{\Omega }}_n\left(t\right)$ is truncated, the corresponding approximation is still unitary when applied to the Schrödinger equation. More generally, if (1) is defined in a Lie group $\mathcal{G}$, then it provides an approximate solution also belonging to $\mathcal{G}$. Moreover, it has also been used to construct efficient numerical integration algorithms also preserving this feature [12,15].

When $A\left(t\right)$ periodically depends on time with period T, the Magnus expansion does not explicitly provide the structure of the solution ensured by the Floquet theorem, i.e., $U\left(t\right)=P\left(t\right){\mathrm{e}}^{tF}$. with $P\left(t\right)$ periodic and F a constant matrix. Nevertheless, it can be conveniently generalized to cover also this situation, in such a way that the matrix $P\left(t\right)$ is expressed as the exponential of a series of periodic terms. The resulting approach (the so-called Floquet–Magnus expansion [10]) has been used during the last years in a variety of physical problems [4,5,7,16].

Very often, the coefficient matrix in (1) is of the form $A\left(t\right)=A_0+\epsilon A_1\left(t\right)$, where $A_0$ is constant, $A_1(t+T)=A_1\left(t\right)$, and $\epsilon >0$ is a (small) parameter. In other words, one is dealing with a time-dependent perturbation of a solvable problem defined by $A_0$. In that case, several perturbative procedures exist to construct $U\left(t\right)$ as a power series in $\epsilon$, either directly (by applying standard perturbation techniques in the interaction picture defined by $A_0$) or taking into account the structure ensured by the Floquet theorem and constructing both matrices $P\left(t\right)$ and F as power series [9]. Of course, if $P\left(t\right)=\exp \left(\mathsf{\Lambda }\right(t\left)\right)$ and $\mathsf{\Lambda }\left(t\right)$ is also constructed as a power series in $\epsilon$, then the qualitative properties of the solution are inherited by the approximations (in particular, they are unitary in quantum evolution problems) [17].

An alternative manner of viewing the Floquet theorem is to interpret the matrix $P\left(t\right)$, provided that it is invertible, as a time-periodic transformation to a new set of coordinates where the new evolution equation has a constant coefficient matrix F, so that $\exp \left(tF\right)$ is the exact solution in the new coordinates. In the language of Quantum Mechanics, this corresponds to a change of picture. This interpretation leads to the important mathematical notion of a reducible system: according to Lyapunov, the general system (1) is called reducible if there exists a matrix $P\left(t\right)$ which together with $det\left(P^{-1}\left(t\right)\right)$ is bounded on $0\le t<+\infty$, such that the system obtained from (1) by applying the linear transformation defined by $P\left(t\right)$ has constant coefficients [18]. In this sense, if $A\left(t\right)$ is a periodic matrix, then (1) is reducible by means of a periodic matrix. The situation is not so clear; however, when $A\left(t\right)$ is quasi-periodic or almost periodic: although several general results exist ensuring reducibility [19,20,21], there are also examples of irreducible systems [22].

In this work, we pursue and generalize this interpretation to show that all of the above mentioned exponential perturbative treatments can be considered as particular instances of a generic change of variables applied to the original differential equation. The idea of making a coordinate change to analyze a problem arises of course in many application settings, ranging from canonical (or symplectic) transformations in classical mechanics (based either on generating functions expressed in terms of mixed, old and new, variables [23], or on the Lie-algebraic setting [24,25]) to changes of picture and unitary transformations in quantum mechanics. What we intend here is to show that several widely used perturbative expansions in quantum mechanics can be indeed derived from the same basic principle using different variations of a unique ansatz based on a generic linear transformations of coordinates. We believe that this interpretation sheds new light into the different expansions and, moreover, allows one to elaborate a unique procedure for analyzing a given problem and compare in an easy way how they behave in practice.

It is important to remark that all the procedures considered here (with the exception of course of the standard perturbation theory) preserve by construction the unitarity of the solution when (1) refers to the Schrödinger equation. More generally, the approximations obtained evolve in the same matrix Lie group as the exact solution of the differential Equation (1).

2. Coordinate Transformations and Linear Systems

To begin with, let us consider the most general situation of a linear differential equation

$$\dot{x}\equiv \frac{dx}{dt}=A\left(t\right)x,x\left(0\right)=x_0\in {\mathbb{C}}^d,$$

(2)

with $A\left(t\right)$ a $d\times d$ matrix whose entries are integrable functions of t. Notice that $U\left(t\right)$ in (1) can be considered the fundamental matrix of (2). We analyze a change of variables $x⟼X$ transforming the original system (2) into

$$\frac{dX}{dt}=F\left(t\right)X,X\left(0\right)=X_0,$$

(3)

where the matrix F adopts some desirable form. Because we are interested in preserving qualitative properties of (2), we impose an additional requirement for the transformation, namely, it has to be of the form

$$x\left(t\right)=\exp \left(\mathsf{\Omega }\right(t\left)\right)X\left(t\right),\mathrm{with}\mathsf{\Omega }\left(0\right)=0,$$

(4)

so that $x\left(0\right)=X\left(0\right)$. Thus, in particular, if (2) is the Schrödinger equation with Hamiltonian $H\left(t\right)=i\hslash A\left(t\right)$, $\mathsf{\Omega }\left(t\right)$ is skew-Hermitian.

It is clear that the transformation (4) is completely determined once the generator $\mathsf{\Omega }\left(t\right)$ is obtained. An evolution equation for $\mathsf{\Omega }$ is obtained by introducing (4) in (2) and also taking into account (3) as

$$\frac{d}{dt}\exp \left(\mathsf{\Omega }\right)=A\left(t\right)\exp \left(\mathsf{\Omega }\right)-\exp \left(\mathsf{\Omega }\right)F\left(t\right).$$

(5)

The derivative of a matrix exponential can be written as [26]

$$\frac{d}{dt}\exp \left(\mathsf{\Omega }\left(t\right)\right)=d{\exp }_{\mathsf{\Omega }\left(t\right)}\left(\dot{\mathsf{\Omega }}\left(t\right)\right)\exp \left(\mathsf{\Omega }\right),$$

(6)

where the symbol $d{\exp }_{\mathsf{\Omega }}\left(C\right)$ stands for the (everywhere convergent) power series

$$d{\exp }_{\mathsf{\Omega }}\left(C\right)=\sum _{k=0}^{\infty }\frac{1}{(k+1)!}{\mathrm{ad}}_{\mathsf{\Omega }}^k\left(C\right)\equiv \frac{\exp \left({\mathrm{ad}}_{\mathsf{\Omega }}\right)-I}{{\mathrm{ad}}_{\mathsf{\Omega }}}\left(C\right).$$

(7)

Here ${\mathrm{ad}}_{\mathsf{\Omega }}^{0}C=C$, ${\mathrm{ad}}_{\mathsf{\Omega }}^{k}C=[\mathsf{\Omega },{\mathrm{ad}}_{\mathsf{\Omega }}^{k-1}C]$, and $[\mathsf{\Omega },C]$ denotes the usual commutator. By inserting (6) into (5), one gets

$$d{\exp }_{\mathsf{\Omega }}\left(\dot{\mathsf{\Omega }}\right)=A-\exp \left(\mathsf{\Omega }\right)F\exp (-\mathsf{\Omega })=A-\exp \left({\mathrm{ad}}_{\mathsf{\Omega }}\right)F,$$

(8)

where

$$\exp \left({\mathrm{ad}}_{\mathsf{\Omega }}\right)F=\sum _{n\ge 0}\frac{1}{n!}{\mathrm{ad}}_{\mathsf{\Omega }}^{n}F={\mathrm{e}}^{\mathsf{\Omega }}F{\mathrm{e}}^{-\mathsf{\Omega }}.$$

If we invert the $d{\exp }_{\mathsf{\Omega }}$ operator given by (7), we get from (8) the formal identity

$$\dot{\mathsf{\Omega }}=\frac{x}{{\mathrm{e}}^x-1}(A-{\mathrm{e}}^{x}F),\mathrm{with}x\equiv {\mathrm{ad}}_{\mathsf{\Omega }}.$$

Now, taking into account that

$$\frac{x{\mathrm{e}}^x}{{\mathrm{e}}^x-1}=\frac{x}{{\mathrm{e}}^x-1}+x,$$

it is clear that

$$\dot{\mathsf{\Omega }}=\frac{x}{{\mathrm{e}}^x-1}(A-F)-xF.$$

(9)

A more convenient way of expressing this identity is obtained by recalling that

$$d{\exp }_{\mathsf{\Omega }}^{-1}\left(C\right)=\frac{{\mathrm{ad}}_{\mathsf{\Omega }}}{\exp \left({\mathrm{ad}}_{\mathsf{\Omega }}\right)-I}\left(C\right)=\sum _{k=0}^{\infty }\frac{B_k}{k!}{\mathrm{ad}}_{\mathsf{\Omega }}^k\left(C\right),$$

where $B_k$ are the Bernoulli numbers, so that (9) reads

$$\dot{\mathsf{\Omega }}=\sum _{k=0}^{\infty }\frac{B_k}{k!}{\mathrm{ad}}_{\mathsf{\Omega }}^k(A-F)-{\mathrm{ad}}_{\mathsf{\Omega }}F.$$

(10)

With more generality, we can assume that both $A\left(t\right)$ and $F\left(t\right)$ are power series of some appropriate parameter $\epsilon$,

$$A\left(t\right)=A_0\left(t\right)+\sum _{n\ge 1}{\epsilon }^n{A}_n\left(t\right),F\left(t\right)=F_0\left(t\right)+\sum _{n\ge 1}{\epsilon }^n{F}_n\left(t\right),$$

(11)

and, thus, the generator $\mathsf{\Omega }$ will be also a power series,

$$\mathsf{\Omega }\left(t\right)=\sum _{n\ge 1}{\epsilon }^n{\mathsf{\Omega }}_n\left(t\right),{\mathsf{\Omega }}_n\left(0\right)=0.$$

(12)

The successive ${\mathsf{\Omega }}_n\left(t\right)$ can be determined by inserting (11) and (12) into Equation (10) and equating terms of the same power in $\epsilon$. Thus, one obtains the following recursive procedure:

$$\begin{array}{cc}\hfill & n=0:F_0=A_0\hfill \\ \hfill & n=1:{\dot{\mathsf{\Omega }}}_1+[{\mathsf{\Omega }}_1,A_0]=W_1^{\left(0\right)}\hfill \\ \hfill & n\ge 2:{\dot{\mathsf{\Omega }}}_n+[{\mathsf{\Omega }}_n,A_0]=W_n^{\left(0\right)}+G_n-V_n,\hfill \end{array}$$

(13)

where

$$\begin{array}{cc}\hfill & W_n^{\left(0\right)}=A_n-F_n,n\ge 1\hfill \\ \hfill & W_n^{\left(k\right)}=\sum _{m=1}^{n-k}[{\mathsf{\Omega }}_m,W_{n-m}^{(k-1)}],n\ge 1,1\le k\le n-1\hfill \\ \hfill & G_n=\sum _{k=1}^{n-1}\frac{B_k}{k!}{W}_n^{\left(k\right)},V_n=\sum _{k=1}^{n-1}[{\mathsf{\Omega }}_k,F_{n-k}],n\ge 2.\hfill \end{array}$$

(14)

Of course, although the change of variables $x⟼X$ is completely general, it is only worth to be considered if Equation (3) is simpler to solve than (2). In the following we analyze several ways of choosing F fulfilling this basic requirement, and how, in this way, we are able to recover different exponential perturbative expansions.

3. Magnus Expansion

The simplest choice one can imagine of is taking $F=0$ or, in other words, one is looking for a linear transformation rendering the original system (2) into

$$\frac{dX}{dt}=0,$$

(15)

with trivial solution $X\left(t\right)=X\left(0\right)=x_0$. A sufficient condition for the reducibility of Equation (2) to (15) is [18]

$${\int }_0^{+\infty }{\parallel A\left(t\right)\parallel }_{F}dt<+\infty ,$$

where ${\parallel A\left(t\right)\parallel }_F=\sqrt{{\sum }_{i,j=1}^d{\left|a_{ij}\right|}^2}$. If this is the case, from (4),

$$x\left(t\right)=\exp \left(\mathsf{\Omega }\left(t\right)\right)X\left(t\right)=\exp \left(\mathsf{\Omega }\left(t\right)\right)x_0,$$

where $\mathsf{\Omega }\left(t\right)$ is determined from (10) with $F=0$, i.e.,

$$\dot{\mathsf{\Omega }}=\sum _{k=0}^{\infty }\frac{B_k}{k!}{\mathrm{ad}}_{\mathsf{\Omega }}^{k}A.$$

(16)

This, of course, corresponds to the well known Magnus expansion for the solution $x\left(t\right)$ of (2) [14,26]. The terms ${\mathsf{\Omega }}_n\left(t\right)$ are then determined by the recursion (14) by taking $F_n=0$. If we take $A_0=0$ and $A_1\left(t\right)=A\left(t\right)$ in (11), then we get the familiar recursive procedure [26]

$$\begin{array}{cc}\hfill & W_1^{\left(0\right)}=A,W_n^{\left(0\right)}=0,n\ge 2\hfill \\ \hfill & W_n^{\left(k\right)}=\sum _{m=1}^{n-k}[{\mathsf{\Omega }}_m,W_{n-m}^{(k-1)}],n\ge 1,1\le k\le n-1\hfill \\ \hfill & G_n=\sum _{k=1}^{n-1}\frac{B_k}{k!}{W}_n^{\left(k\right)},n\ge 2\hfill \\ \hfill & {\dot{\mathsf{\Omega }}}_1=A_1,{\dot{\mathsf{\Omega }}}_n=G_n,n\ge 2,\hfill \end{array}$$

(17)

whence the successive terms ${\mathsf{\Omega }}_n$ are obtained by integration. An explicit expression for ${\mathsf{\Omega }}_n\left(t\right)$ involving only independent nested commutators can be obtained by working out this recurrence and using the class of bases proposed by Dragt and Forest [27] for the Lie algebra generated by the operators $A_1\left(t_1\right),\dots A_1\left(t_n\right)$. Specifically, one has [28]

$$\begin{array}{cc}\hfill & {\mathsf{\Omega }}_n\left(t\right)=\frac{1}{n}\sum _{\sigma \in S_{n-1}}{(-1)}^{d_{\sigma }}\frac{1}{\left(\genfrac{}{}{0pt}{}{n-1}{d_{\sigma }}\right)}{\int }_0^{t}d{t}_1{\int }_0^{t_1}d{t}_2\cdots {\int }_0^{t_{n-1}}d{t}_n\hfill \\ \hfill & [A\left(t_{\sigma \left(1\right)}\right),[A\left(t_{\sigma \left(2\right)}\right)\cdots [A\left(t_{\sigma (n-1)}\right),A\left(t_n\right)]\cdots ]],\hfill \end{array}$$

(18)

where $\sigma$ is a permutation of $\{1,2,\dots ,n-1\}$ and $d_{\sigma }$ corresponds to the number descents of $\sigma$. We recall that $\sigma$ has a descent in i if $\sigma \left(i\right)>\sigma (i+1)$, $i=1,\dots ,n-2$. Notice that the argument in the last term is fixed to $t_n$, and one considers all of the permutations in $A\left(t_1\right),A\left(t_2\right),\dots ,A\left(t_{n-1}\right)$. Moreover, the series (12) converges in this case in the interval $t\in [0,t_{\mathrm{f}}]$, such that

$${\int }_0^{t_{\mathrm{f}}}{\parallel A\left(s\right)\parallel }_{2}ds<\pi$$

(19)

and the sum $\mathsf{\Omega }\left(t\right)$ satisfies $\exp \left(\mathsf{\Omega }\right(t\left)\right)=U\left(t\right)$ [26]. Here, ${\parallel ·\parallel }_2$ denotes the spectral norm, characterized as ${\parallel A\parallel }_2=max\{\sqrt{\lambda }:\lambda \mathrm{is}\mathrm{an}\mathrm{eigenvalue}\mathrm{of}{A}^{†}A\}$.

The Magnus expansion has a long history as a tool to approximate solutions in a wide spectrum of fields in Physics and Chemistry, from atomic and molecular physics to Nuclear Magnetic Resonance to Quantum Electrodynamics (see [26] and references therein). Additionally, in computational mathematics, it has been used to construct efficient algorithms for the numerical integration of differential equations within the widest field of geometric numerical integration [12,13,15]. Recently, it has also been used in the treatment of dissipative driven two-level systems [29].

4. Floquet–Magnus Expansion

The Magnus expansion can be, in principle applied for any particular time dependence in $A\left(t\right)$, as long as the integrals in ${\mathsf{\Omega }}_n\left(t\right)$ can be computed or conveniently approximated. When $A\left(t\right)$ in (2) is periodic in t with period T; however, other changes of variables may be more suitable. According to Floquet’s theorem, the original system is reducible to a system with a constant coefficient matrix F, whose eigenvalues (the so-called characteristic exponents) determine the asymptotic stability of the solution $x\left(t\right)$. In addition, the linear transformation is periodic with the same period T [6].

In our general framework, then, it makes sense to determine a change of variables $x=\exp \left(\mathsf{\Omega }\right(t\left)\right)X\left(t\right)$ in such a way that F in (3) is constant, so that $X\left(t\right)=\exp \left(tF\right)x_0$ and $\mathsf{\Omega }\left(t\right)$ is periodic. Proceeding as before, if we take $A_0=0$ and $A_1\left(t\right)=A\left(t\right)$ in (11), the procedure (13)–(14) simplifies to

$$\begin{array}{cc}\hfill & W_1^{\left(0\right)}=A-F_1,W_n^{\left(0\right)}=-F_n,n\ge 2\hfill \\ \hfill & W_n^{\left(k\right)}=\sum _{m=1}^{n-k}[{\mathsf{\Omega }}_m,W_{n-m}^{(k-1)}],n\ge 1,1\le k\le n-1\hfill \\ \hfill & G_n=\sum _{k=1}^{n-1}\frac{B_k}{k!}{W}_n^{\left(k\right)},V_n=\sum _{k=1}^{n-1}[{\mathsf{\Omega }}_k,F_{n-k}],n\ge 2\hfill \end{array}$$

(20)

and

$$\begin{array}{cc}\hfill & n=1:{\dot{\mathsf{\Omega }}}_1=A-F_1\hfill \\ \hfill & n\ge 2:{\dot{\mathsf{\Omega }}}_n=-F_n+G_n-V_n\equiv {\mathcal{F}}_n-F_n\hfill \end{array}$$

(21)

Notice that, in general, ${\mathcal{F}}_n=G_n-V_n$ depends on the previously computed ${\mathsf{\Omega }}_k$, $F_k$, $k=1,\dots ,n-1$, so that equations (21) can be solved recursively, as follows. First, we determine $F_1$ and $F_n$ by taking the average of A and ${\mathcal{F}}_n$, respectively, over one period T,

$$F_1=\frac{1}{T}{\int }_0^{T}A\left(s\right)ds,F_n=\langle {\mathcal{F}}_n\rangle \equiv \frac{1}{T}{\int }_0^T(G_n\left(s\right)-V_n\left(s\right))ds,$$

and then compute the integrals

$${\mathsf{\Omega }}_1\left(t\right)={\int }_0^{t}A\left(s\right)ds-t{F}_1,{\mathsf{\Omega }}_n\left(t\right)={\int }_0^t(G_n\left(s\right)-V_n\left(s\right))ds-t{F}_n,$$

respectively, thus ensuring that ${\mathsf{\Omega }}_n$ is periodic with the same period T. This results in the well known Floquet–Magnus expansion for the solution of (2),

$$x\left(t\right)=\exp \left(\mathsf{\Omega }\left(t\right)\right)\exp \left(tF\right)x_0,$$

originally introduced in [10] and subsequently applied in different areas [4,5,7,16]. In the context of periodic quantum mechanical problems, $H_{\mathrm{ef}}\equiv i\hslash F$ is called the effective Hamiltonian of the problem. This expansion presents the great advantage that, in addition to preserving unitarity as the Magnus expansion, also allows one to determine the stability of the system by analyzing the eigenvalues of F.

As shown in [10], the resulting series for $\mathsf{\Omega }\left(t\right)$ is absolutely convergent at least for $t\in [0,t_{\mathrm{f}}]$, such that

$${\int }_0^{t_{\mathrm{f}}}{\parallel A\left(s\right)\parallel }_{2}ds<0.20925\dots$$

The procedure can be easily generalized to quasi-periodic problems [30]. We recall that $A\left(t\right)$ is quasi-periodic in t with frequencies $({\omega }_1,\dots ,{\omega }_r)$ if $A\left(t\right)=\tilde{A}({\theta }_1,\dots ,{\theta }_r)$, where $\tilde{A}$ is $2\pi$-periodic with respect to ${\theta }_1,\dots ,{\theta }_r$ and ${\theta }_j={\omega }_{j}t$ for $j=1,\dots ,r$. In that case, we can write

$$A\left(t\right)=\sum _{k\in {\mathbb{Z}}^r}{a}_k{\mathrm{e}}^{i(k,\omega )t}$$

where $(k,\omega )\equiv k_1{\omega }_1+\cdots +k_r{\omega }_r$ and $\sum |a_k{|}^2<\infty$ [31].

In this case, ${\mathcal{F}}_n$ is also quasi-periodic (by induction) and we take

$$F_n=\langle {\mathcal{F}}_n\rangle \equiv \underset{T\to \infty }{lim}\frac{1}{T}{\int }_a^{a+T}{\mathcal{F}}_n\left(s\right)ds,$$

the limiting mean value of ${\mathcal{F}}_n\left(t\right)$, independent of the particular value of a. In consequence,

$${\dot{\mathsf{\Omega }}}_n=-F_n+{\mathcal{F}}_n\left(t\right)=\sum _{k\in {\mathbb{Z}}^r\backslash \left\{0\right\}}{f}_{j,k}{\mathrm{e}}^{i(k,\omega )t}$$

and

$${\mathsf{\Omega }}_n\left(t\right)={\int }_0^t{\mathcal{F}}_n\left(s\right)ds-t{F}_n$$

is also quasi-periodic with the same basic frequencies as ${\mathcal{F}}_n$.

5. Perturbed Problems

There are many problems that can be formulated as (2) with $A\left(t\right)$, depending on a perturbation parameter $\epsilon$,

$$A\left(t\right)=A_0+\sum _{n\ge 1}{\epsilon }^n{A}_n\left(t\right),$$

so that $A_0$ is a constant matrix and $A_n\left(t\right)$ are periodic or quasi-periodic functions of t. Here, again the issue of reducibility has received much attention along the years, and the problem consists in determining both the linear transformation $P\left(t\right)$ and the constant matrix F. In this section, we consider several possible ways to proceed, depending on the final matrix F one is aiming for.

5.1. Removing the Perturbation

One obvious approach is to take $F=A_0$. In other words, we try to construct a transformation rendering the original system into another one in which the perturbation is removed [17]. Therefore, the transformation fulfilling this requirement is

$$x\left(t\right)=\exp \left(\mathsf{\Omega }(t,\epsilon )\right)\exp \left(t{A}_0\right)x_0,\mathsf{\Omega }(t,\epsilon )=\sum _{n\ge 1}{\epsilon }^n{\mathsf{\Omega }}_n\left(t\right),$$

(22)

where we have written explicitly the dependence on $\epsilon$, the parameter of the perturbation. The recurrence for determining the generators ${\mathsf{\Omega }}_n$ is obtained from (13)–(14) as

$$\begin{array}{cc}\hfill & W_n^{\left(0\right)}=A_n,n\ge 1\hfill \\ \hfill & W_n^{\left(k\right)}=\sum _{m=1}^{n-k}[{\mathsf{\Omega }}_m,W_{n-m}^{(k-1)}],n\ge 1,1\le k\le n-1\hfill \\ \hfill & G_1=0,G_n=\sum _{k=1}^{n-1}\frac{B_k}{k!}{W}_n^{\left(k\right)},n\ge 2\hfill \\ \hfill & {\dot{\mathsf{\Omega }}}_n+[{\mathsf{\Omega }}_n,A_0]=A_n+G_n,n\ge 1.\hfill \end{array}$$

(23)

Alternatively, we can write

$${\dot{\mathsf{\Omega }}}_n={\mathrm{ad}}_{A_0}{\mathsf{\Omega }}_n+A_n+G_n$$

(24)

with solution verifying ${\mathsf{\Omega }}_j\left(0\right)=0$ given by

$${\mathsf{\Omega }}_n\left(t\right)={\mathrm{e}}^{t{\mathrm{ad}}_{A_0}}{\int }_0^t{\mathrm{e}}^{-s{\mathrm{ad}}_{A_0}}\left(A_n\left(s\right)+G_n\left(s\right)\right)ds.$$

As a matter of fact, this scheme can be related with the usual perturbation treatment of Equation (2). To keep the formalism simple, let us assume that $A\left(t\right)$ in (11) is $A\left(t\right)=A_0+\epsilon A_1\left(t\right)$ and work directly with the power series of the transformation,

$$P(t,\epsilon )=\exp \left(\mathsf{\Omega }(t,\epsilon )\right)=\sum _{n\ge 0}{\epsilon }^n{P}_n\left(t\right),\mathrm{with}{P}_0=I.$$

Subsequently, by inserting $x\left(t\right)=P(t,\epsilon ){\mathrm{e}}^{t{A}_0}{x}_0$ into Equation (2), one determines the differential equation satisfied by P as

$$\dot{P}+P{A}_0=AP,$$

whence the successive terms $P_n\left(t\right)$ verify

$${\dot{P}}_n+[P_n,A_0]=A_1{P}_{n-1},P_n\left(0\right)=0,n\ge 1.$$

The solution, as with Equation (24), is given by

$$P_n\left(t\right)={\mathrm{e}}^{t{A}_0}{g}_n\left(t\right){\mathrm{e}}^{-t{A}_0},$$

with

$$\begin{array}{cc}\hfill g_n\left(t\right)& ={\int }_0^{t}ds{\mathrm{e}}^{-s{A}_0}{A}_1\left(s\right)P_{n-1}\left(s\right){\mathrm{e}}^{s{A}_0}={\int }_0^{t}ds{\mathrm{e}}^{-s{A}_0}{A}_1\left(s\right){\mathrm{e}}^{s{A}_0}{g}_{n-1}\left(s\right)\hfill \\ \hfill & ={\int }_0^{t}d{s}_1{\int }_0^{s_1}d{s}_2\dots {\int }_0^{s_{n-1}}d{s}_n{A}_I\left(s_1\right)A_I\left(s_2\right)\dots A_I\left(s_n\right).\hfill \end{array}$$

(25)

Here, $g_0\equiv I$ and we have denoted $A_I\left(s\right)\equiv {\mathrm{e}}^{-s{A}_0}{A}_1\left(s\right){\mathrm{e}}^{s{A}_0}$. Therefore, the solution of (2) reads

$$\begin{array}{cc}\hfill x\left(t\right)& =P(t,\epsilon ){\mathrm{e}}^{t{A}_0}{x}_0=\left(I+\sum _{n\ge 1}{\epsilon }^n{P}_n\left(t\right)\right){\mathrm{e}}^{t{A}_0}{x}_0=\left({\mathrm{e}}^{t{A}_0}+\sum _{n\ge 1}{\epsilon }^n{\mathrm{e}}^{t{A}_0}{g}_n\left(t\right)\right)x_0\hfill \\ & ={\mathrm{e}}^{t{A}_0}\left(I+\sum _{n\ge 1}{\epsilon }^n{g}_n\left(t\right)\right)x_0\hfill \end{array}$$

and so, taking into account the explicit expression (25) for $g_n\left(t\right)$, one recovers the solution provided by the standard perturbation theory in the picture defined by $A_0$.

5.2. Lie–Deprit Perturbative Algorithm

The previous treatment has one important drawback (in addition to the lack of preservation of unitarity when one deals with the Schrödinger equation and the expansion is truncted), namely, when $A\left(t\right)$ is a periodic or quasi-periodic function of time, the secular terms are not removed and, as a result, the time interval of validity of the resulting approximations is very small indeed. For this reason it is worth considering in the general framework (13) and (14) the quasi-periodic case and look for a quasi-periodic transformation leading to a constant coefficient matrix

$$F\left(\epsilon \right)=F_0+\sum _{n\ge 1}{\epsilon }^n{F}_n,\mathrm{with}{F}_0=A_0.$$

(26)

This problem has been addressed a number of times in the literature (see, e.g., [21,32] and references therein). In particular, in [33] a perturbative algorithm is presented for constructing this transformation as the exponential of a quasi-periodic matrix. Here, we show how the results of [17,33] can be reproduced by the generic scheme proposed here in a more direct way, in the sense that the terms ${\mathsf{\Omega }}_j$ are determined at once from (13)–(14).

It is clear that the problem reduces to solve

$${\dot{\mathsf{\Omega }}}_n+[{\mathsf{\Omega }}_n,A_0]={\mathcal{F}}_n-F_n,{\mathsf{\Omega }}_n\left(0\right)=0$$

where now

$${\mathcal{F}}_1\equiv A_1,{\mathcal{F}}_n\equiv A_n+G_n-V_n,n\ge 2,$$

and the goal is to determine the constant term $F_n$ and construct ${\mathsf{\Omega }}_n\left(t\right)$ as a quasi-periodic function with the same basic frequencies as $A_n$ for $n\ge 1$.

As shown in Appendix A, the solution verifying all of these requirements is given by

$$F_n=\langle {\mathcal{F}}_n\left(t\right)\rangle -[A_0,M_n\left(0\right)],{\mathsf{\Omega }}_n\left(t\right)=-M_n\left(0\right)+{\mathrm{e}}^{t{\mathrm{ad}}_{A_0}}{M}_n\left(t\right),$$

where $M_n\left(t\right)$ is the antiderivative

$$M_n\left(t\right)=\int {\mathrm{e}}^{-t{\mathrm{ad}}_{A_0}}\left({\mathcal{F}}_n\left(t\right)-\langle {\mathcal{F}}_n\left(t\right)\rangle \right)dt.$$

Finally, the solution of (2) is written as

$$x\left(t\right)=\exp \left(\mathsf{\Omega }(t,\epsilon )\right)\exp \left(tF\left(\epsilon \right)\right)x_0,$$

where $\mathsf{\Omega }(t,\epsilon )$ and $F\left(\epsilon \right)$ are appropriate truncations of the corresponding series.

5.3. Quantum Averaging

In the context of time-dependent quantum mechanical systems, the averaging method has been used to construct quantum analogues of classical perturbative treatments [34,35]. Essentially, the basic approach in that setting is to transform the original Hamiltonian of the problem by a unitary transformation, so that the problem of finding the time evolution of the transformed Hamiltonian is easier than the original one. The idea has been also applied in the perturbative treatment of pulse-driven quantum problems [36,37,38].

This approach to quantum averaging can be fit into our general framework of Section 2 by taking $F_0=A_0$ and the terms $F_n\left(t\right)$ in (11) verifying in addition

$${\dot{F}}_n\left(t\right)+[F_n\left(t\right),A_0]=0,n\ge 1.$$

(27)

Clearly, the solution of (27) verifies $F_n\left(t\right)={\mathrm{e}}^{t{\mathrm{ad}}_{A_0}}{F}_n\left(0\right)={\mathrm{e}}^{t{A}_0}{F}_n\left(0\right){\mathrm{e}}^{-t{A}_0}$ and, thus,

$$\frac{d}{dt}\left({\mathrm{e}}^{-t{A}_0}X\left(t\right)\right)=\left(F\left(0\right)-A_0\right){\mathrm{e}}^{-t{A}_0}X\left(t\right),$$

which leads to $\exp (-t{A}_0)X\left(t\right)=\exp (t(F\left(0\right)-A_0)X\left(0\right)$ and finally

$$X\left(t\right)={\mathrm{e}}^{t{A}_0}{\mathrm{e}}^{t(F\left(0\right)-A_0)}X\left(0\right).$$

In other words, condition (27) guarantees that Equation (3) can be solved in a closed way and the dynamics of (2) is obtained once the generator $\mathsf{\Omega }\left(t\right)$ is determined, even when $F_n$ depends explicitly on time. As shown in [35], the corresponding solutions are

$$F_n\left(t\right)=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_0^T{B}_n(\sigma ,t)d\sigma ,{\mathsf{\Omega }}_n\left(t\right)=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_0^{T}d\tau {\int }_0^{\tau }d\sigma \left(B_n(\sigma ,t)-F_n\left(t\right)\right),$$

where

$$B_n(\sigma ,t)={\mathrm{e}}^{(t-\sigma )A_0}{\mathcal{F}}_n(t-\sigma ){\mathrm{e}}^{-(t-\sigma )A_0}$$

provided

$$\underset{T\to \infty }{lim}\frac{1}{T}\left({\mathcal{F}}_n\left(t\right)-{\mathrm{e}}^{-T{A}_0}{\mathcal{F}}_n(t-T){\mathrm{e}}^{T{A}_0}\right)=0.$$

Here ${\mathcal{F}}_1=A_1$ and ${\mathcal{F}}_n=A_n+G_n-V_n$ for $n\ge 2$.

6. Illustrative Examples

In previous sections, we have reviewed several perturbative expansions aimed to get analytic approximations for the solution of (1), and how they can be derived from a same basic principle, namely a transformation of coordinates. This allows one, in particular, to design a unique computational procedure to deal with a particular problem defined by $A\left(t\right)$ and take the most appropriate variant in each case. To illustrate the technique we next consider two examples describing the dynamics of simple time-dependent quantum systems, although the same treatment can be applied of course to other more involved problems.

6.1. The Three-Lambda System

The so-called driven three-lambda system describes an atomic three energy-level system with two ground states $|1\rangle$, $|2\rangle$ with the same energy $E_1$ that are coupled with an excited state $|3\rangle$ with energy $E_3$ via a time-dependent laser field. In the interaction picture, the Hamiltonian is given by $H\left(t\right)=f\left(t\right)|3\rangle \left(\langle 1|+\langle 2|\right)$ (+ Hermitian conjugate), or equivalently,

$$H\left(t\right)=\left(\begin{array}{ccc}0& 0& \overline{f}\left(t\right)\\ 0& 0& \overline{f}\left(t\right)\\ f\left(t\right)& f\left(t\right)& 0\end{array}\right),$$

(28)

where $f\left(t\right)$ is usually a periodic function. The corresponding Schrödinger equation, $i\hslash \dot{U}\left(t\right)=H\left(t\right)U\left(t\right)$, is then a particular case of Equation (2) with $A\left(t\right)=-iH\left(t\right)/\hslash$, and one is typically interested in obtaining the induced probability of transition between states $|1\rangle$ and $|2\rangle$, namely $P_{12}\left(t\right)={|\langle 1|U\left(t\right)|2\rangle |}^2$.

We follow [30] and take

$$f\left(t\right)=\beta {\mathrm{e}}^{i\omega t}$$

(29)

with $\beta =1$ and $\omega =\frac{10}{{(1+\frac{\sqrt{2}}{2})}^{\frac{1}{2}}}\approx 7.65$ as an example of periodic function. Starting from the initial condition $U\left(0\right)=I$, we compute the approximate solution until the final time $\omega t=400$, i.e., for 487 periods, with the Floquet–Magnus (FM) expansion (20)–(21) up to $n=3$ and $n=4$ and determine the corresponding transition probability. In this way, we get Figure 1, where the result achieved by the effective Hamiltonian $X\left(t\right)={\mathrm{e}}^{tF}$, with $F=-i{H}_{\mathrm{ef}}/\hslash$, is also depicted. We see that, for this value of $\omega$, considering more terms in the series for both $\mathsf{\Omega }\left(t\right)$ and F leads to better results and that working only with the effective Hamiltonian gives indeed a good approximation. The reference solution is computed with the DSolve function of Mathematica. For completeness, the value of the effective Hamiltonian up to $n=4$ reads

$$H_{\mathrm{ef}}=\left(\begin{array}{ccc}\frac{6-{\omega }^2}{{\omega }^4}& \frac{6-{\omega }^2}{{\omega }^4}& \frac{4}{{\omega }^3}\\ \frac{6-{\omega }^2}{{\omega }^4}& \frac{6-{\omega }^2}{{\omega }^4}& \frac{4}{{\omega }^3}\\ \frac{4}{{\omega }^3}& \frac{4}{{\omega }^3}& \frac{2\left({\omega }^2-6\right)}{{\omega }^4}\end{array}\right),$$

whereas, the transition probability obtained with $H_{\mathrm{ef}}$ up to $n=3$ is given by ($\hslash =1$)

$$P_{12}^{\left(3\right)}=\frac{{sin}^2\left(t\tilde{\omega }\right)}{({\omega }^2+8)}\left(-4cos\left(2t\tilde{\omega }\right)+{\omega }^2+4\right)$$

with $\tilde{\omega }\equiv \left(\sqrt{{\omega }^2+8}\right)/{\omega }^3$.

Next, we repeat the simulation with $\omega =6$ and $\omega =12$, by taking $n=3$ terms in both the Floquet–Magnus and the Magnus expansion (17). The results are shown in Figure 2 (top). Notice that the approximations get worse for smaller values of $\omega$. To get a more quantitative view, also in Figure 2 (bottom) we show the error in the transition probability (in logarithmic scale) for the same values of $\omega$, but now include $n=3$ and $n=7$ terms in the FM expansion, as well as the result obtained with just the effective Hamiltonian in this last case. As is clear from the figure, taking into account more terms in the FM expansion improves the approximations, and this improvement is more remarkable for larger values of $\omega$.

Here, the approximate solution is obtained as

$$U\left(t\right)=\exp \left(\tilde{\mathsf{\Omega }}\left(t\right)\right)\exp (-it{H}_{\mathrm{ef}}/\hslash )$$

and the two parts of U encode different behaviors of the system. Thus, $\exp \left(\tilde{\mathsf{\Omega }}\left(t\right)\right)$ describes fast fluctuations around the envelope evolution provided by $\exp (-it{H}_{\mathrm{ef}}/\hslash )$, but these fluctuations are not always negligible, especially if the driving frequency is not too large. In any case, it is worth noticing that taking only into consideration the effective Hamiltonian renders less accurate results than working with the full approximation.

Next, we take $f\left(t\right)$ in (28) to be a quasi-periodic function, namely

$$f\left(t\right)=\beta ({\mathrm{e}}^{i\omega t}+{\mathrm{e}}^{i\sqrt{2}\omega t}),$$

(30)

with $\beta =1$ and $\omega =12$, and compute approximations using the Floquet–Magnus and Magnus expansions with $n=2$ and $n=3$ terms. The corresponding errors in the transition probability are shown in Figure 3. We see that, in this case, including more terms in the expansions does not improve significantly the accuracy of the approximations. This is an indication of the lack of convergence of the expansions for these values of the parameters. In fact, it is straightforward to verify that ${\parallel A\left(t\right)\parallel }_2=\sqrt{8}\beta \left|cos\left(\frac{\sqrt{2}-1}{2}\omega t\right)\right|$, so that for the values of the parameters considered, the convergence of the Magnus expansion is ensured for $t<1.6117$ according to the estimate (19). To illustrate the behavior in this range of times, in Figure 4 we collect the results obtained by taking up to $n=6$ terms in the expansion. It is worth noticing the systematic decrease in the error when more terms are included in the expansion. A similar pattern is observed for the Floquet–Magnus expansion, but with much smaller times (the convergence is only assured for $t<0.074$).

6.2. Bloch–Siegert Dynamics

This is an example of a periodically driven two-level quantum system described by the time-dependent Hamiltonian

$$H\left(t\right)=\left(\begin{array}{cc}{\omega }_0/2& 2bcos\left(\omega t\right)\\ 2bcos\left(\omega t\right)& -{\omega }_0/2\end{array}\right).$$

(31)

The coupling parameter b is the amplitude of a driving radio-frequency field and $\omega$ is its frequency. We can write

$$H\left(t\right)=H_0+\epsilon H_1\left(t\right)$$

with $\epsilon =2b$ and

$$H_0=\frac{{\omega }_0}{2}{\sigma }_3,H_1\left(t\right)=cos\left(\omega t\right){\sigma }_1,$$

where ${\sigma }_i$ are Pauli matrices. With this simple example, we can illustrate the different treatments previously considered, both perturbative and non-perturbative, and also analyze how the results vary with the perturbation parameter $\epsilon$.

The use of the Magnus expansion for the approximate analytic treatment of time-dependent two-level systems has a long history. In particular, the work of Pechukas and Light [39] includes, for the first time, an specific analysis of the convergence of the procedure for this class of systems.

As in [40], we consider the resonant case $\omega ={\omega }_0=1$ and two different values of the perturbation, $\epsilon =0.2$ and $\epsilon =1$. One could use the Magnus expansion (17) directly to the matrix $A\left(t\right)=-iH\left(t\right)/\hslash$. However, in that case, the convergence of the series, as given by (19), is only ensured for times $t\in [0,t_{\mathrm{f}}]$, such that

$${\int }_0^{t_{\mathrm{f}}}\sqrt{\frac{1}{4}+{\epsilon }^2{cos}^2\left(t\right)}dt<\pi .$$

Thus, in particular, if $\epsilon =0.2$, then $t_{\mathrm{f}}\approx 6.056$, whereas, for $\epsilon =1$, it reduces to $t_{\mathrm{f}}\approx 3.608$. We cannot expect, then, accurate results for larger time intervals, and this is indeed what is observed in practice. Something analogous also happens for the Floquet–Magnus expansion, with the difference that the interval of convergence is even smaller. It makes sense, in consequence, to apply both procedures in the interaction picture, i.e., to the matrix $A_I\left(t\right)={\mathrm{e}}^{-t{A}_0}{A}_1\left(t\right){\mathrm{e}}^{t{A}_0}$ and this is what we have done in our computations by applying the procedure of Section 4. We also compare with the Lie–Deprit (LD) perturbation algorithm and the formalism removing the perturbation, i.e., such that $F=A_0$. As in the previous example, we compute $|U_{12}{\left(t\right)|}^2$, i.e., the probability of transition between the two levels, obtained with each procedure with $n=3$ terms. The corresponding results are shown in Figure 5 as a function of time, where we also include the exact result for comparison. Left diagrams are obtained with FM (in the interaction picture), whereas we depict in the right diagrams the results obtained by LD and also when F is taken as $A_0$. Notice that top and bottom panels correspond to $\epsilon =0.2$ and $\epsilon =1$, respectively. We have also depicted the result obtained with just the effective Hamiltonian in the case of the Floquet–Magnus expansion.

It can be observed that both the LD perturbative algorithm and the one obtained by removing the perturbation provide accurate results only when $\epsilon$ takes small values, in contrast with Floquet–Magnus in the interaction picture, with an enlarged validity domain in $\epsilon$. Notice that, also in this case, the effective Hamiltonian is not enough to get a precise description of the dynamics when $\epsilon$ increases. We should remark that all schemes render unitary approximations by construction.

For comparison, in Figure 6 we collect the results achieved by the standard Magnus expansion in the interaction picture for this same problem and values of the perturbation parameter: $\epsilon =0.2$ (left) and $\epsilon =1$ (right). Here again, increasing the value of $\epsilon$ leads to a loss of accuracy for long times. If the Magnus expansion is applied to the matrix $A\left(t\right)$, then only accurate results are obtained for times in the convergence domain.

To analyze how incorporating an increasing number of terms in the Floquet–Magnus expansion leads to an improvement of the approximation, even for large values of $\epsilon$, we show in Figure 7 the results for the transition probability for $\epsilon =1.5$ obtained with $n=3,5,7,9$ terms (left to right, top to bottom) as a function of time. Notice that, even for this large value of the perturbation parameter, we still have a correct description of the dynamics when a sufficiently large number of terms is taken into account.

To get a more quantitative description of the different approximations, in Figure 8 we depict the error in the transition probability committed by the Floquet–Magnus expansion (left) and Lie–Deprit (right) with $n=3,5,7,9$ terms each for $\epsilon =0.5$. Notice that taking into account more terms leads to a reduction in the error, and the LD approximation clearly deteriorates with time, whereas FM still provides excellent results, even for large time intervals.

7. Concluding Remarks

The subject of the perturbative treatment of linear systems of time-dependent differential equations has a long history in physics and mathematics. In physics, in particular, it plays a fundamental role in the study of the evolution of quantum systems. Among the different procedures, the so-called exponential perturbation algorithms have the relevant property of preserving the unitary character of the evolution. More in general, the approximations they provide belong to the same Lie group than the exact solution of the problem. Archetypical examples of exponential perturbation theories are the Magnus expansion since its inception in the 1950s and, more recently, the Floquet–Magnus expansion, several quantum averaging procedures and a generalization of the well known Hori–Deprit perturbation theory of classical mechanics. Each of these algorithms has been derived in an independent way and it is not always easy to establish connections and common points between them.

The present paper tries to bridge this gap by showing that all of them can be seen in fact as the result of linear changes of coordinates expressed as the exponential of a certain series whose terms can be obtained recursively. In addition to the Magnus and Floquet–Magnus expansions, other techniques can also be incorporated into our general framework, including the quantum averaging method and the Lie–Deprit perturbative algorithm. In the process, we have also considered a novel approach, namely an exponential transformation rendering the original system into another one in which the perturbation is removed. The resulting approximations preserve whatever qualitative features the exact solution may possess (unitarity, orthogonality, symplecticness, etc.). Even the standard perturbation theory in the interaction picture is recovered in this setting when the exponential defining the transformation is truncated.

With this same framework one might of course consider other exponential transformations, and this would automatically lead to new perturbation formalisms that might be particularly adapted to the problem at hand. For instance, we could choose $F\left(t\right)$ as a diagonal time-dependent matrix, so that equation (3) is easy to solve, and construct the corresponding generator $\mathsf{\Omega }\left(t\right)$. We believe that the proposed framework sheds new light into nature of the different expansions and allows one to create a unique procedure to analyze a given problem, obtaining all of the expansions by choosing appropriately the coefficient matrix in the new variables. Moreover, it also allows one to determine in an easy way what is the best algorithm for a given problem. This feature has been illustrated in this work by applying the procedure to two simple examples.

We have seen, in particular, that for problems with a periodic time dependence, the Floquet–Magnus expansion leads to more accurate results over longer time intervals, even for perturbed problems of the form $A\left(t\right)=A_0+\epsilon A_1\left(t\right)$ in the interaction picture, unless the parameter $\epsilon$ is exceedingly small, in which other purely perturbative procedures are competitive. It would be interesting to check whether this conclusion remains valid for more involved problems.

The interested reader may find useful the file available in [41] containing the Mathematica implementation of the techniques exposed here for the quantum mechanical problem describing the Bloch–Siegert dynamics.

In this paper we have always taken $X\left(0\right)=x\left(0\right)=x_0$, but it is clear that we can also take instead $X\left(0\right)$ as the image by the transformation of $x_0$. In that case, the solution is factorized as

$$x\left(t\right)=\exp \left(\mathsf{\Omega }\left(t\right)\right)\exp \left(tF\right)X\left(0\right)=\exp \left(\mathsf{\Omega }\left(t\right)\right)\exp \left(tF\right)\exp (-\mathsf{\Omega }\left(0\right))x_0,$$

where, obviously, the recurrences for determining $\mathsf{\Omega }\left(t\right)$ and F are slightly different from those collected here. The final results can also vary, as shown in [33].

Although only examples of quantum mechanics have been considered here, the formalism of this paper can be applied in fact to any explicitly time-dependent linear system of differential equations. This might include the treatment of open quantum systems described by a non-Hermitian Hamiltonian operator [42]. This will be the subject of future research.