## Abstract

We formulate a unitary perturbation theory for quantum mechanics inspired by the Lie-Deprit formulation of canonical transformations. The original Hamiltonian is converted into a solvable one by a transformation obtained through a Magnus expansion. This ensures unitarity at every order in a small parameter. A comparison with the standard perturbation theory is provided. We work out the scheme up to order ten with some simple examples.

## 1. Introduction

Historically, classical mechanics played a crucial role in the early formulation of quantum mechanics. Afterwards, the new ideas led to a quantum formalism that loosened contact with the classical language. But nevertheless, every now and then, mutual relationships between them have been profitably considered.

Fixing our attention on the treatment of the effect of a perturbation on a solvable physical system, matrix mechanics adapted, with daring modifications and reinterpretations, the classical lore of canonical transformations, generating functions, action and angle variables and the like. This is apparent, e.g. in the famous paper by Born *et al.* (1926). Nowadays, both time-independent and time-dependent quantum perturbation theories are presented in the Hilbert space formalism and make no mention of their classical counterparts.

Suppose one is interested in a system with Hamiltonian *H*(*t*,*ϵ*)=*H*_{0}+*H*′(*t*,*ϵ*), where *H*_{0} is time-independent and *H*′(*t*,*ϵ*) is a time-dependent perturbation, depending on a small parameter *ϵ* that controls the intensity of the perturbation in such a way that *H*′(*t*,*ϵ*=0)=0. The aim is then to determine the quantum time-evolution operator *U*(*t*,*t*_{0}), which takes the state of the system from time *t*_{0} to time *t*. One of its most salient features is its unitary character. The evolution operator satisfies the Schrödinger equation which, in the conventional approach, is solved by series expansion in *ϵ*. After truncation at a given order, it allows approximation of the different observable magnitudes of physical interest. The main drawback of the method is that any finite order approximation of *U*(*t*,*t*_{0}) ceases to be unitary.

By contrast, the procedure followed in the founding years of quantum mechanics and directly inherited from the classical perturbation theory was more geometric in character. The power series expansions were also the basic tool, but the goal was to construct a geometric (unitary) transformation that renders solvable the Hamiltonian *H*(*t*,*ϵ*) (e.g. Birtwistle 1928). The reader acquitted with the Hamilton–Jacobi equation will recognize the idea of finding a canonical transformation (namely, a change of phase space coordinates) that transforms a Hamiltonian *H* into a Hamiltonian *K* easily solvable in the new coordinates. The problem is then translated to the search of the corresponding generating function.

The original canonical transformations theory used in classical mechanics relies on the effective use of generating functions expressed in terms of mixed, old and new, variables (e.g. Arnold *et al.* 2006). It means that, once the equations of motion have been integrated, one has to solve implicit functional equations so as to express the results in terms of either the old or the new variables. In general, this is difficult to do.

This classical canonical theory, however, fitted in the late 1960s into a Lie-algebraic setting, mainly in connection with problems of celestial mechanics. In that respect, the paper (Deprit 1969) *Canonical transformations depending on a small parameter* played a seminal role by introducing what it is known nowadays as the *Lie-Deprit perturbation method*. Essentially, a canonical transformation is generated by a new ‘Hamiltonian’ *W* in such a way that the shift by ‘time’ *ϵ* along the trajectories of *W* produces the required transformation. In this way, no mixed variables appear in the formalism. The algorithm can be readily implemented in a computer, thus allowing carrying out of perturbative analysis up to high order.

Since then, techniques and procedures of classical perturbation theory have been used in a quantum mechanical context from time to time (e.g. Kummer 1971; Eckhardt 1986; Scherer 1994; Daems *et al.* 2004). In fact, this is quite natural, given the existing parallelism between perturbation theory in classical and quantum systems when both are formulated in terms of geometric transformations (canonical maps and unitary transformations, respectively). The common formal setting is indeed provided by the language of Lie algebras. From this perspective, the perturbation problem defined by *H*_{0}+*H*′(*t*,*ϵ*) can be solved by changing through a (canonical or unitary) transformation to a Hamiltonian *K* whose dynamics is easier to solve, at least up to a certain order in the *ϵ* parameter. When the perturbation expansion converges, the original problem is completely solved. Even when the convergence problem is not settled, the procedure very often provides very accurate results.

The purpose of the present paper consists in adapting and developing an implementation of the Lie-Deprit method that gives results amenable to analyse time-dependent problems in quantum mechanics. For this reason, we have adapted the title of the 1969 Deprit paper. The procedure leads to operator linear differential equations, which we solve by use of the Magnus expansion. This will ensure unitarity at any order of truncation in *ϵ*. This is an important feature of the algorithm proposed here. Thus, it is guaranteed that the constructed approximation shares several important qualitative features with the exact solution. In addition, in this formalism, several approximation schemes are possible, better suited to the mathematical nature of the problem at hand. We will explain how to take advantage of this. The connection with the standard time-dependent perturbation theory will also be provided. Here, we will proceed at a purely formal level, keeping the technical requisites to a minimum.

The plan of the paper is the following: in §2, we briefly summarize the standard time-dependent perturbation theory in quantum mechanics in the interaction picture defined by *H*_{0}. In §3, we develop the new unitary formalism, paying special attention to the construction of the transformation and the derivation of the basic equation. The procedure is highly flexible in the sense that several possible choices of the new Hamiltonian *K* are possible. Also the time-independent case fits naturally into the algorithm. These issues are analysed in §4, whereas §5 is devoted to a pair of examples to illustrate the main features of the scheme. Finally, §6 contains an additional discussion and further comparisons with the so-called quantum averaging approach.

## 2. Time-dependent perturbation theory in quantum mechanics

Let us consider the dynamics governed by the Hamiltonian
2.1where we have assumed that *H*′(*t*,*ϵ*) is analytic in *ϵ*. In most situations, the power series expansion of *H*′(*t*,*ϵ*) collapses to the first term: *H*′(*t*,*ϵ*)=*ϵH*_{1}(*t*).

The corresponding time-evolution operator *U*(*t*,*t*_{0}) is unitary, i.e. *U*(*t*,*t*_{0})*U*^{†}(*t*,*t*_{0})=*U*^{†}(*t*,*t*_{0})*U*(*t*,*t*_{0})=*I*, where the dagger means adjoint. This property guarantees, in physical terms, probability conservation. Although for simplicity we have not shown it explicitly, *U*(*t*,*t*_{0}) depends also on *ϵ*. In turn, it obeys the Schrödinger equation
2.2

When the Hamiltonian is independent of time, i.e. *H*(*t*)=*H*_{0}, the solution of (2.2) reads . Otherwise, the solution is by no means so simple, unless and *H*(*t*) do commute.

For the Hamiltonian (2.1), the problem can be approximately solved by the conventional time-dependent perturbation theory, yielding the solution of (2.2) in power series of *ϵ*. In this setting, it is very important to isolate the piece of the problem to be treated approximately and allow the solvable part of the Hamiltonian to be exactly integrated. A way to deal with this question consists in transforming the problem into the interaction picture, i.e. we factorize
2.3where the unknown operator *U*_{I} obeys the Schrödinger equation
2.4with the new Hamiltonian given by
2.5In this way, the interaction picture defined in (2.3) provides a unitary transformation that allows the mathematical integration of the solvable piece *H*_{0} in *H*(*t*,*ϵ*) and focuses the approximations just on *H*′(*t*,*ϵ*). Explicitly, one deals with the expansion
2.6where *U*_{n}(*t*,*t*_{0}) stands for the contribution of order *n*.

Substitution of (2.6) in (2.4) and then equating terms of the same power *ϵ*^{n} yields the successive terms *U*_{n}. In particular, the first two read
2.7Unfortunately, neither the approximant *U*_{I}≃*I*+*U*_{1} nor any of higher finite order are unitary. This, of course, has undesirable consequences in practice. For instance, the computed transition probability between different quantum states may exceed unity. Only when the infinite series is summed up, the unitary character is restored.

## 3. Unitary time-dependent perturbation theory

### (a) The formalism

We will start by rephrasing, in modern quantum mechanics language, the main idea of perturbation theory such as it was first formulated in matrix mechanics. It will help us to introduce the central issue of the present article.

Suppose we have a *hard* time-independent quantum problem defined by the Hamiltonian *H*(*ϵ*)=*H*_{0}+*ϵH*_{1}+*ϵ*^{2}*H*_{2}+⋯ , and assume the mathematical problem defined by *H*_{0} has been solved. This is equivalent to finding a basis of the Hilbert space such that *H*_{0} is diagonal. Then the idea is to look for a unitary transformation *T*(*ϵ*) (or ‘canonical transformation’ in the language of matrix mechanics) such that the new Hamiltonian *K*=*T*^{†}*HT* is diagonal. The transformation *T* can be obtained recursively as a power series in the perturbation parameter *ϵ*. Notice that the starting point of the method is rather geometric in nature, in clear contrast with time-dependent perturbation theory that is algebraic from the very beginning. Of course, an appropriate choice of *K* is essential.

In the following, we consider the more general, time-dependent problem defined by (2.1). Assume that the dynamics corresponding to *H*_{0} has been solved, i.e. we have determined
Let us now define
3.1where *T*(*t*,*ϵ*) is a unitary transformation we want to construct in order to solve our problem.

The Schrödinger equation for *U*_{K}(*t*,*t*_{0}) reads
3.2where the new Hamiltonian is
3.3

Our goal is to determine *T*(*t*,*ϵ*) in such a way that (3.2) is easier to solve than the original equation (2.2). Notice that the factorization (3.1) may be considered as a generalization of (2.3). Very often we will take *T*(*t*_{0},*ϵ*)=*I* in (3.1) and *t*_{0}=0.

In perturbation theory, the perturbed and the non-perturbed Hamiltonians are always very close to each other in some mathematical sense. It is right that what makes meaningful perturbative computations. It is then natural to determine *T*(*t*,*ϵ*) as a unitary near-identity transformation
Two different ways of guaranteeing these features are the following:

— Introduce a skew-Hermitian operator such that .

— Alternatively, introduce a skew-Hermitian operator

*L*(*t*,*ϵ*) such that*T*(*t*,*ϵ*) is the solution of the operator differential equation 3.4

This second prescription is precisely the approach proposed by Deprit in classical mechanics, and the one we follow here. Somehow, the perturbed Hamiltonian may be regarded as the *evolution*, with respect to the parameter *ϵ*, of the non-perturbed Hamiltonian (which may be thought of as an *initial condition*). Deprit's idea consists in materializing such an abstract *ϵ*-evolution in a manner similar to the usual *t*-evolution, since this way the unitary character of *T* is ensured. To this end, an evolution law similar to the one in (2.2) is introduced. In (3.4), *L*(*t*,*ϵ*) is the (unknown) generator of the transformation *T*, namely, of the *ϵ*-evolution. Equivalently, one has
3.5

The assessment of the situation is the following. Initially, we are given a time-dependent Hamiltonian *H*(*t*,*ϵ*)=*H*_{0}+*H*′(*t*,*ϵ*), and the problem is to solve the Schrödinger equation. However, we are only able to solve the piece *H*_{0}. Then we introduce the idea of an *ϵ*-evolution in such a way that a (unknown) continuous transformation *T*(*t*,*ϵ*) takes *H*(*t*,*ϵ*) to a new Hamiltonian *K*(*t*,*ϵ*) to be determined in such a way that the corresponding Schrödinger equation (3.2) can be solved. The inverse transformation *T*^{†}(*t*,*ϵ*) may be used then to transform the solution of *K*(*t*,*ϵ*) into the one for *H*(*t*,*ϵ*). However, by introducing the law of the *ϵ*-evolution, we have added a further unknown: the generator *L*(*t*,*ϵ*). Thus, the situation may now appear worse than at the beginning: neither *T* nor *L* are known!

This apparently adverse scheme will provide, however, an excellent starting point to deal with power series expansions and preserve the unitary character even if the series are truncated.

We first notice that *L* and *T*^{†} are related by the linear differential equation (3.5). Thus, once *L* has been determined, we can obtain *T* and *T*^{†} by formally applying the Magnus expansion (Magnus 1954) to (3.5). As is well known, his proposal is to solve a generic operator linear equation *dY* (*t*)/*dt*=*A*(*t*)*Y* (*t*) with *Y* (0)=*I* in the form
The exponent *Ω*(*t*) is given by an infinite series of terms that are linear combinations of integrals and nested commutators involving *m*-fold products of the operator *A* at different times (Blanes *et al.* 2009).

Although explicit formulae for *Ω*_{m} exist, it is sometimes more convenient to generate them by a recursive procedure. The one proposed in Klarsfeld & Oteo (1989) is particularly well suited from a computational point of view (Blanes *et al.* 2009):
and
3.6where *B*_{j} stand for Bernoulli numbers and [⋅,⋅] represents the commutator, [*A*,*B*]≡*AB*−*BA*. Notice that this approach can be directly applied to obtain approximate formal solutions of equation (3.5), so that we can write
3.7where *Ω* is a skew-Hermitian operator. Thus, the two above enumerated alternatives to render *T* unitary are, in fact, equivalent.

Next, we derive equation (3.3) with respect to *ϵ*,
3.8Since this expression involves both the operators *T* and *L*, it is convenient to insert (3.7) here. In fact, it is straightforward to show that
3.9Here we have defined the operator ad_{Ω}*B*≡[*Ω*,*B*], with . In consequence, we have
3.10an equation formulated only in terms of *L*,*H* and *K* that is easier to handle than (3.8).

In the present formalism, then, three different problems have to be studied:

Choose the new Hamiltonian

*K*such that equation (3.2) is easy to solve.Compute the skew-Hermitian generator

*L*of the required transformation to the new picture.Construct the unitary transformation

*T*from the generator*L*, or equivalently, the operator*Ω*in .

The first two problems enumerated earlier can be solved perturbatively with equation (3.10), whereas the third one can be tackled independently. We treat these problems in reverse order.

### (b) Determining *T*

We next proceed to compute *Ω*, once *L* has been fixed. To do that, we introduce, in addition to the expression (2.1) for the Hamiltonian *H*, the following series expansions:
3.11and
3.12It turns out that *Ω*(*t*,*ϵ*) in (3.7) can also be determined as a power series in *ϵ*:
3.13with *v*_{n}(*t*) expressed in terms of *L*_{j}(*t*). This can be done algorithmically with the Magnus expansion up to any given order. In fact, the procedure already used to construct numerical integrators based on Magnus expansion (Blanes *et al.* 2009) is of use here also. Specifically, if the series (3.12) is inserted into the recurrence (3.6), we get up to order *ϵ*^{5}
Reordering according with the power of *ϵ*, we have the following expressions for the first terms in the series (3.13):
3.14This procedure can be carried out for higher orders with a symbolic algebra package. On the other hand, if the series (3.13) is inserted into (3.9), then
3.15and, again, this series can be determined algorithmically. The first terms read
3.16It is worth stressing here that the series (3.13) and (3.15) have to be determined only once. In fact, computing these series up to *n*=10 takes only a few minutes on a personal computer. Terms up to *n*=10 and *n*=9, respectively, can be found at the website (http://www.gicas.uji.es/Research/Perturbation.html) as *Mathematica* files.

### (c) Determining *L*

When the series expansions for *K*, *L* and (3.15) are substituted into (3.10), collecting terms of the same power in *ϵ* yields a system of recurrence equations. Specifically, we arrive at *K*_{0}=*H*_{0} and
3.17with
3.18The next step in the procedure is to propose a suitable *K*_{n} and solve (3.17) to determine *L*_{n} and finally the transformation *T*, which is unitary at any order by construction.

Equation (3.17) has the following formal solution:
3.19where *n*>0, and we have taken *t*_{0}=0 for simplicity. Very often we will fix *L*_{n}(0)=0, so that
3.20

## 4. On the choice of *K*

The choice of a particular *K* is a degree of freedom of the method but it may also constitute a nuisance. It might be the case that two particular *H* and *K* are associated in a natural way. Here we will comment about a number of possible choices of *K* in the general case.

### (a) *K*=*H*_{0}

Perhaps, this is the most immediate option. It means that systematically
Equivalently, one tries to construct a unitary transformation in such a way that in the new picture there is no perturbation at all. In that case, and
4.1Notice that with this choice, the whole *ϵ* dependency is contained in *T*.

### (b) *K* diagonal

Suppose that *H*_{0} has a pure non-degenerate point spectrum. In that case, we can choose *K*_{n} to be diagonal. In other words, we perturbatively construct a unitary transformation to a new picture where the Hamiltonian *K* is diagonal and therefore equation (3.2) is easy to solve. Specifically,
and finally

### (c) Time-independent case

The formalism developed earlier works also when *H*′ is time-independent. In that case, it makes sense to construct a time-independent unitary transformation *T*(*ϵ*) so that
A natural choice is then to take the new Hamiltonian so that [*H*_{0},*K*]=0. This is analogous to the standard procedure of constructing the Birkhoff–Gustavson normal form, both in the classical and quantum setting (Ali 1985; Eckhardt 1986). From (3.17), the equations to solve read then
4.2It has been shown in Jauslin *et al.* (2000) that if these equations are defined in the space of bounded skew-Hermitian multiplication operators, then the formal solution can be written in terms of averages as
and
even if *H*_{0} has degenerate eigenvalues. As a matter of fact, these expressions reproduce the standard results (Birtwistle 1928; Wu 1986). If |*ν*,*j*〉 denotes an eigenbasis of *H*_{0}, where *ν* labels the distinct eigenvalues and *j* distinguishes different basis vectors of the degeneracy subspace, then
whereas
The operator evolution in this setting reads
being unitary even upon truncation of the series. This coincides with the formal solution proposed by Kummer (1971) to avoid secular terms in perturbative approximations. In this sense, the formalism developed here could be considered as a generalization to time-dependent problems of Kummer's approach.

## 5. Examples

In this section, we apply the previous formalism to two simple examples with the aim of illustrating its main features.

### Example 5.1

First, we consider the Hamiltonian
5.1where *σ*_{j} denote the Pauli matrices and *ϵ*, *ω*_{0}, *ω*≠*ω*_{0} are real parameters. The exact solution of the Schrödinger equation (2.2) with *U*(0)=*I* is then
Let 1 and 2 denote the spin up and down states, respectively. Then the exact transition probability from state 1 to state 2 is given by
and may serve as a test to analyse the various results in the following.

Our perturbative scheme starts by writing *H*(*t*)=*H*_{0}+*ϵH*_{1}(*t*), with
We have applied the procedure up to *n*=10 twice. First by choosing *K*_{n}=0 for *n*≥1 and second by taking *K*_{n} as a diagonal matrix. In either case, we get
and the corresponding evolution operator *U*_{K}(*t*), whereas
leads to *Ω*(*t*,*ϵ*) up to this order. Finally,

In figure 1, we depict the difference between the exact and approximate transition probabilities as a function of time computed with this procedure up to *n*=4 and *n*=10. Solid lines correspond to choosing *K*=*H*_{0}, whereas dot-dashed curves are generated by taking *K*_{n} as a diagonal matrix. Notice that increasing the number of terms in *L*(*t*,*ϵ*) in both cases provides a more accurate approximation.

### Example 5.2

The second example illustrates the time-independent case. We take the two-level spin Hamiltonian
Then the corresponding equations (4.2) provide for the first orders the following results:
Carrying out the computation up to order *ϵ*^{10}, we get
5.2On the other hand, the exact eigenvalues of the Hamiltonian *H* are in this case and its expansion in power series of *ϵ* agrees with (5.2). Finally, *Ω*(*ϵ*)=i*ω*_{a}*σ*_{2}, with
and thus the approximate solution is given by .

## 6. Discussion

### (a) Connection with the standard time-dependent perturbation theory

The formalism we have developed here contains as a limiting case the usual time-dependent perturbation theory summarized in §2. This indeed corresponds to choosing *K*=*H*_{0} as a new Hamiltonian and expanding the operator *T*(*t*,*ϵ*)=e^{−Ω(t,ϵ)} in power series of *ϵ*. Obviously, the resulting scheme is no longer unitary when the series is truncated.

Let us illustrate this point by computing explicitly the first order in *ϵ*. Since , it is clear from (3.20) that
whereas
But then
with
This, as we know, is precisely the result achieved by standard perturbation theory in the interaction picture defined by *H*_{0}.

The previous treatment can be generalized to any order in *ϵ*. In fact, by following this approach, we have explicitly reproduced up to order *ϵ*^{6} the results achieved by the usual time-dependent perturbation theory for example 5.1.

### (b) Quantum averaging

In a series of papers published in the mid-1990s, Scherer used the averaging method to construct quantum mechanical analogues of the classical perturbation expansions by Poincaré and Kolmogorov, both in the autonomous case and for time-dependent operators (Scherer 1994, 1997, 1998). The idea was further elaborated by Daems *et al.* (2003*a*, 2004) and Sugny *et al.* (2004), and successfully applied in the perturbative treatment of pulse-driven quantum problems. This quantum averaging algorithm can also be reproduced with our formalism as follows.

Scherer's basic proposal also consists of transforming the original Hamiltonian (2.1) with the help of a unitary transformation *T*(*t*,*ϵ*) such that the problem of finding the time evolution of the transformed Hamiltonian *K*(*t*,*ϵ*) is easier than the original one. In particular, he chooses the terms *K*_{n} in the series defining *K* so as to satisfy
6.1Then, it can be shown that the time evolution generated by *K*(*t*,*ϵ*) is given by (*t*_{0}=0)
In other words, the time evolution is reduced to a time-independent problem if one introduces an interaction picture defined by *H*_{0} for the new Hamiltonian *K*. Finally,
6.2whereas equation (6.1) can be solved by averaging, generalizing the autonomous case:
6.3The unitary transformation *T*(*t*,*ϵ*) is then computed order by order from *L*_{n}(*t*). In fact, a recursive procedure is proposed for determining the terms *T*_{n}(*t*) in the series
6.4If not all terms *L*_{n} can be computed, the corresponding truncated series (6.4) provides approximations that are no longer unitary. It is argued, however, that these violations of unitarity do not grow with time, since they do not contain secular terms.

Notice that Scherer's algorithm fits quite naturally in the present scheme. Moreover, with the choice (6.3), we can determine *L*_{n}(*t*) by applying (3.20) and then compute the transformation as with the series (3.13). Now, contrarily to the original formulation, the resulting scheme is unitary upon truncation.

We have applied this modified scheme again to the Hamiltonian (5.1) and compared with the previous choices of *K*_{n}=0 and *K*_{n} diagonal (*n*>0). In figure 2, we represent the corresponding errors in the transition probability obtained with the three procedures when *n*=10 terms are considered in the series. The results illustrate the fairly different character of each approximation. It should be stressed again that the approximation (6.2) with (6.3) considered here is a unitary modification of that originally proposed by Scherer (1997, 1998), i.e. with instead of (6.4).

### (c) Open problems and further comments

There are several issues only touched upon in this work that deserve a much more detailed analysis. First, we have been only concerned here with the formal aspects of the transformations and expansions. It is clear, then, that a more rigorous treatment is necessary, with special attention to the convergence of the expansion. Second, the procedure developed here is fairly flexible since it enables different choices of the new Hamiltonian. We have seen on a simple example that different options lead indeed to different approximations; so the natural question is: for any particular problem, which is the optimal choice? In other words, which particular *K* allows one to get the best approximation once the expansion is truncated at a given order? In this sense, it is worth mentioning that there are several free parameters in the formalism that eventually may be used to tune the algorithm and improve the accuracy (e.g. Daems *et al.* 2004).

In dealing with quantum time-dependent problems, one also finds in the literature the methods that are iterative in character. So, for example, when the Hamiltonian is slowly varying with time, sequences of unitary transformations intended to diagonalize instantaneously the Hamiltonian lead to non-perturbative approaches (Berry 1987, 1990). In turn, the so-called Fer expansion (Blanes *et al.* 2009) is based on the iterative use of the transformation that defines the interaction picture. A perturbation theory can also be used in an iterative way: the Hamiltonian obtained in each step is used as unperturbed Hamiltonian for the next step. This quantum mechanical version of the KAM scheme has been shown to produce very accurate results in treatments, different to the one proposed in this paper, both for time-independent and time-dependent quantum perturbations (Scherer 1997; Daems *et al.* 2003*a*,*b*, 2004). The same idea could be implemented with the procedure developed here. In this work, by following closely the Lie-Deprit method in classical mechanics, we have developed a unitary perturbative algorithm well adapted to quantum problems. The technique provides a, in practice approximate, single unitary transformation designed to simplify the dynamics as much as possible.

## Acknowledgements

This work has been partially supported by MICINN (Spain), under grant no. AYA2010-22111-C03-02 (J.A.O. and J.R.) and MTM2010-18246-C03-02 (F.C.), co-financed by FEDER Funds of the European Union.

- Received June 20, 2011.
- Accepted October 4, 2011.

- This journal is © 2011 The Royal Society