Closed-Form Expressions for the Matrix Exponential

We discuss a method to obtain closed-form expressions of f(A), where f is an analytic function and A a square, diagonalizable matrix. The method exploits the Cayley–Hamilton theorem and has been previously reported using tools that are perhaps not sufficiently appealing to physicists. Here, we derive the results on which the method is based by using tools most commonly employed by physicists. We show the advantages of the method in comparison with standard approaches, especially when dealing with the exponential of low-dimensional matrices. In contrast to other approaches that require, e.g., solving differential equations, the present method only requires the construction of the inverse of the Vandermonde matrix. We show the advantages of the method by applying it to different cases, mostly restricting the calculational effort to the handling of two-by-two matrices.


Introduction
Physicists are quite often faced with the task of calculating f (A), where A is an n × n matrix and f an analytic function whose series expansion generally contains infinitely many terms.The most prominent example corresponds to exp A. Usual approaches to calculate f (A) consist in either truncating its series expansion, or else finding a way to "re-summate" terms so as to get a closed-form expression.There is yet another option that can be advantageously applied when dealing with an n × n matrix, and which derives from the Cayley-Hamilton theorem [1].This theorem states that every square matrix satisfies its characteristic equation.As a consequence of this property, any series expansion can be written in terms of the first n powers of A. While this result is surely very well known among mathematicians, it appears to be not so widespread within the physicists' community [2].Indeed, most textbooks on quantum mechanics still resort to the Baker-Hausdorff lemma or to special properties of the involved matrices, in order to obtain closed-form expressions of series expansions [3][4][5].This happens even when dealing with low-dimensional matrices, i.e., in cases in which exploiting the Cayley-Hamilton theorem would straightforwardly lead to the desired result.Such a state of affairs probably reflects a lack of literature on the subject that is more palatable to physicists than to mathematicians.The present paper aims at dealing with the subject matter by using language and tools that are most familiar to physicists.No claim of priority is made; our purpose is to show how well the derived results fit into the repertoire of tools that physicists routinely employ.To this end, we start addressing the simple, yet rich enough case of 2 × 2 matrices.
An archetypical example is the Hamiltonian H = kσ • B that rules the dynamics of a spin-1/2 particle subjected to a magnetic field B. Here, σ = (σ x , σ y , σ z ) denotes the Pauli spin operator and k is a parameter that provides the above expression with appropriate units.The upsurge of research in several areas of physics-most notably in quantum optics-involving two-level systems, has made a Hamiltonian of the above type quite ubiquitous.Indeed, the dynamics of any two-level system is ruled by a Hamiltonian that can be written in such a form.Hence, one often requires an explicit, closed-form expression for quantities such as exp(iαn • σ), where n is a unit vector.This closed-form expression can be obtained as a generalization of Euler's formula exp iα = cos α + i sin α.It reads with I denoting the identity operator.
Let us recall how most textbooks of quantum mechanics proceed to demonstrate Equation (1) (see, e.g., [3][4][5]).The demonstration starts by writing the series expansion exp A = k A k /k! for the case A = iαn • σ.Next, one invokes the following relationship: whose proof rests on σ i σ j = δ ij I + i ijk σ k (summation over repeated indices being understood).Equation (2) implies that (n • σ) 2n = I, and hence (n • σ) 2n+1 = n • σ.This allows one to split the power series of exp(iαn • σ) in two parts, one constituted by even and the other by odd powers of iαn • σ: By similarly splitting Euler's exponential, i.e., one sees that Equation (3) is the same as Equation (1).Although this standard demonstration is a relatively simple one, it seems to be tightly related to the particular properties of the operator n • σ, as well as to our ability to "re-summate" the series expansion so as to obtain a closed-form expression.There are several other cases [6] in which a relation similar to Equation (1) follows as a consequence of generalizing some properties of the group SU (2) and its algebra to the case SU (N ), with N > 2. Central to these generalizations and to their associated techniques are both the Cayley-Hamilton theorem and the closure of the Lie algebra su(N ) under commutation and anti-commutation of its elements [6].As already recalled, the Cayley-Hamilton theorem states that any n × n matrix A satisfies its own characteristic equation p(A) = 0, where is A's characteristic polynomial.From p(A) = 0 it follows that any power A k , with k ≥ n, can be written in terms of the matrices I = A 0 , A, . . ., A n−1 .Thus, any infinite series, such as the one corresponding to exp A, may be rewritten in terms of the n powers A 0 , A, . . ., A n−1 .By exploiting this fact one can recover Equation (1).Reciprocally, given A, one can construct a matrix B that satisfies exp B = A, as shown by Dattoli, Mari and Torre [2].These authors used essentially the same tools as we do here and presented some of the results that we will show below, but leaving them in an implicit form.The aforementioned authors belong to a group that has extensively dealt with our subject matter and beyond it [7], applying the present techniques to cases of current interest [8].A somewhat different approach was followed by Leonard [9], who related the Cayley-Hamilton theorem to the solution of ordinary differential equations, in order to get closed expressions for the matrix exponential.This technique can be applied to all n × n matrices, including those that are not diagonalizable.Untidt and Nielsen [10] used this technique when addressing the groups SU (2), SU (3) and SU (4).Now, especially when addressing SU (2), Leonard's approach seems to be unnecessarily involved.This is because there is a trade-off between the wide applicability of the method and its tailoring to a special case.When dealing with diagonalizable matrices, the present approach may prove more useful.Thus, one exploits not only the Cayley-Hamilton theorem, but the diagonalizability of the involved matrices as well.As a result, we are provided with a straightforward way to obtain closed-form expressions for the matrix exponential.There are certainly many other ways that are either more general [9,11] or else better suited to specific cases [12][13][14][15][16], but the present method is especially useful for physical applications.The rest of the paper is organized as follows.First, we present Leonard's technique in a way that somewhat differs from the approach used in [9].Thereafter, we show how to obtain Equation (1) by using a technique that can be generalized to diagonalizable n × n matrices, thereby introducing the method that is the main subject of the present work.As an illustration of this technique, we address some representative cases that were taken from the repertoire of classical mechanics, quantum electrodynamics, quantum optics and from the realm of Lorentz transformations.While the results obtained are known, their derivations should serve to demonstrate the versatility of the method.Let us stress once again that our aim has been to present this method by following an approach that could be appealing to most physicists, rather than to mathematically oriented readers.

Closed Form of the Matrix Exponential via the Solution of Differential Equations
Consider the coupled system of differential equations, given by with x = (x 1 , . . ., x n ) T and A a constant, n × n matrix.The matrix exponential appears in the solution of Equation ( 6), when we write it as x(t) = e At x(0).By successive derivation of this exponential we obtain D k e At = A k e At .Hence, p(D)e At ≡ (D n + c n−1 D n−1 + . . .+ c 1 D + c 0 )e At = p(A)e At = 0, on account of p(A) = 0, i.e., the Cayley-Hamilton theorem.Now, as already noted, this implies that e At can be expressed in terms of A 0 , A, . . ., A n−1 .Let us consider the matrix M (t) := n−1 k=0 y k (t)A k , with the y k (t) being n independent solutions of the differential equation p(D)y(t) = 0.That is, the y k (t) solve this equation for n different initial conditions that will be conveniently chosen.We have thus that p(D)M (t) = n−1 k=0 p(D)y k (t)A k = 0. Our goal is to choose the y k (t) so that e At = M (t).To this end, we note that It is then clear that we must take the following initial conditions: In such a case, e At and M (t) satisfy both the same differential equation and the same initial conditions.Hence, e At = M (t).
Summarizing, the method consists in solving the n-th order differential equation p(D)y(t) = 0 for n different initial conditions.These conditions read D j y k (0) = δ j k , with j, k ∈ {0, . . ., n − 1}.The matrix exponential is then given by e At = n−1 k=0 y k (t)A k .The standard procedure for solving p(D)y(t) = 0 requires finding the roots of the characteristic equation p(λ) = 0.Each root λ with multiplicity m contributes to the general solution with a term (a 0 + a 1 t + . . .+ a m−1 t m−1 )e λt , the a k being fixed by the initial conditions.As already said, this method applies even when the matrix A is not diagonalizable.However, when the eigenvalue problem for A is a solvable one, another approach can be more convenient.We present such an approach in what follows.

Closed Form of the Matrix Exponential via the Solution of Algebraic Equations
Let us return to Equation (1).We will derive it anew, this time using standard tools of quantum mechanics.Consider a Hermitian operator A, whose eigenvectors satisfy A |a k = a k |a k and span the Hilbert space on which A acts.Thus, the identity operator can be written as for any function F (A) that can be expanded in powers of A.
Let us consider the 2 × 2 case A = n • σ, with n a unit vector.This matrix has the eigenvalues ±1 and the corresponding eigenvectors |n ± .That is, n • σ |n ± = ± |n ± .We need no more than this to get Equation (1).Indeed, from n The operator iαn • σ has eigenvectors |n ± and eigenvalues ±iα.Thus, which is Equation (1).Note that it has not been necessary to know the eigenvectors of A = iαn • σ.
It is a matter of convenience whether one chooses to express exp (iαn • σ) in terms of the projectors |n ± n ± |, or in terms of I and n • σ.
Let us now see how the above method generalizes when dealing with higher-dimensional spaces.To this end, we keep dealing with rotations.The operator exp (iαn • σ) is a rotation operator acting on spinor space.It is also an element of the group SU (2), whose generators can be taken as X i = iσ i /2, i = 1, 2, 3.They satisfy the commutation relations [X i , X j ] = ijk X k that characterize the rotation algebra.The rotation operator can also act on three-dimensional vectors r.In this case, one often uses the following formula, which gives the rotated vector r in terms of the rotation angle θ and the unit vector n that defines the rotation axis: Equation ( 11) is usually derived from vector algebra plus some geometrical considerations [17].We can derive it, alternatively, by the method used above.To this end, we consider the rotation generators X i for three-dimensional space, which can be read off from the next formula, Equation (12).The rotation matrix is then obtained as exp (θn • X), with It is straightforward to find the eigenvalues of the non-Hermitian, antisymmetric matrix M .They are 0 and ±i.Let us denote the corresponding eigenvectors as |n 0 and |n ± , respectively.Similarly to the spin case, we have now We need a third equation, if we want to express the three projectors |n k n k |, k = ±, 0, in terms of I and M .This equation is obtained by squaring M : From Equations ( 13)-( 15 By letting M , as given in Equation ( 12), act on r = (x, y, z) T , we easily see that M r = n × r and Thus, on account of Equation ( 17), r = exp(θM )r reads the same as Equation (11).The general case is now clear.Consider an operator A whose matrix representation is an N × N matrix.Once the eigenvalues a k of A (which we assume nondegenerate) have been determined, we can write the N equations: from which it is possible to obtain the N projectors |a k a k | in terms of I, A, A 2 , . . ., A N −1 .To this end, we must solve the system The matrix in Equation ( 18), with components , is a Vandermonde matrix, whose inverse can be explicitly given [18].Once we have written the |a k a k | in terms of I, A, . . .A N −1 , we can express any analytic function of A in terms of these N powers of A, in particular For the case N = 4, for instance, we have the following result: ) The general solution can be written in terms of the inverse of the Vandermonde matrix V .To this end, consider a system of equations that reads like (18), but with the operators entering the column vectors being replaced by numbers, i.e., |a j a j | → w j , with j = 1, . . ., N , and A k → q k+1 , with k = 0, . . ., N − 1.The solution of this system is given by w j = N −1 k=0 U j,k q k , with U = V −1 , the inverse of the Vandermonde matrix.This matrix inverse can be calculated as follows [18].Let us define a polynomial P j (x) of degree N − 1 as The coefficients U j,k of the last equality follow from expanding the preceding expression and collecting equal powers of x.These U j,k are the components of V −1 .Indeed, setting x = a i and observing that , we see that U is the inverse of the Vandermonde matrix.The projectors |a j a j | in Equation ( 18) can thus be obtained by replacing x → A in Equation (23).We get in this way the explicit solution The above expression can be inserted into Equation (7), if one wants to write F (A) in terms of the first N powers of A.
So far, we have assumed that the eigenvalues of A are all nondegenerate.Let us now consider a matrix M with degenerate eigenvalues.As before, we deal with a special case, from which the general formalism can be easily inferred.Let M be of dimension four and with eigenvalues λ 1 and λ 2 , which are two-fold degenerate.We can group the projectors as follows: It is then easy to solve the above equations for the two projectors associated with the two eigenvalues.We obtain We can then write We will need this result for the calculation of the unitary operator that defines the Foldy-Wouthuysen transformation, our next example.It is now clear that in the general case of degenerate eigenvalues, we can proceed similarly to the nondegenerate case, but solving n < N equations.

Examples
Let us now see how the method works when applied to some well-known cases.Henceforth, we refer to the method as the Cayley-Hamilton (CH)-method, for short.Our aim is to show the simplicity of the required calculations, as compared with standard techniques.

The Foldy-Wouthuysen Transformation
The Foldy-Wouthuysen transformation is introduced [19] with the aim of decoupling the upper (ϕ) and lower (χ) components of a bispinor ψ = (ϕ, χ) T that solves the Dirac equation i ∂ψ/∂t = Hψ, where H = −i cα • ∇ + βmc 2 .Here, β and α = (α x , α y , α z ) are the 4 × 4 Dirac matrices: The Foldy-Wouthuysen transformation is given by ψ = U ψ, with [19] We can calculate U by applying Equation ( 29 The standard way to get this result requires developing the exponential in a power series.Thereafter, one must exploit the commutation properties of α and β in order to group together odd and even powers of θ.This finally leads to the same closed-form expression that we have arrived at after some few steps.

Lorentz-Type Equations of Motion
The dynamics of several classical and quantum systems is ruled by equations that can be cast as differential equations for a three-vector S.These equations often contain terms of the form Ω×. An example of this is the ubiquitous equation Equation ( 35) and its variants have been recently addressed by Babusci, Dattoli and Sabia [20], who applied operational methods to deal with them.Instead of writing Equation (35) in matrix form, these authors chose to exploit the properties of the vector product by defining the operator Ω := Ω×.The solution for the case ∂Ω/∂t = 0, for instance, was obtained by expanding exp(t Ω) as an infinite series and using the cyclical properties of the vector product in order to get S(t) in closed form.This form is nothing but Equation (11) with the replacements r → S(t), r → S(0) and θ → Ωt, where Ω := |Ω|.We obtained Equation (11) without expanding the exponential and without using any cyclic properties.Our solution follows from writing Equation (35) in matrix form, i.e., where M is given by Equation ( 12) with n = Ω/Ω.The solution S(t) = exp(M Ωt)S(0) is then easily written in closed form by applying the CH-method, as in Equation (11).The advantages of this method show up even more sharply when dealing with some extensions of Equation (36).Consider, e.g., the non-homogeneous version of Equation (35): This is the form taken by the Lorentz equation of motion when the electromagnetic field is given by scalar and vector potentials reading Φ = −E • r and A = B × r/2, respectively [20].The solution of Equation ( 37) is easily obtained by acting on both sides with the "integrating (operator-valued) factor" exp(−ΩM t).One then readily obtains, for the initial condition S(0) = S 0 , The matrix exponentials in Equation ( 38) can be expressed in their eigenbasis, as in Equation ( 16).For a time-independent N , the integral in Equation ( 38) is then trivial.An equivalent solution is given in [20], but written in terms of the evolution operator Û (t) = exp(t Ω) and its inverse.Inverse operators repeatedly appear within such a framework [20] and are often calculated with the help of the Laplace transform identity: Â−1 = ´∞ 0 exp(−s Â)ds.Depending on Â, this could be not such a straightforward task as it might appear at first sight.Now, while vector notation gives us additional physical insight, vector calculus can rapidly turn into a messy business.Our strategy is therefore to avoid vector calculus and instead rely on the CH-method as much as possible.Only at the end we write down our results, if we wish, in terms of vector products and the like.That is, we use Equations ( 13)-( 17) systematically, in particular Equation ( 16) when we need to handle exp(θM ), e.g., within integrals.The simplification comes about from our working with the eigenbasis of exp(θM ), i.e., with the eigenbasis of M .Writing down the final results in three-vector notation amounts to expressing these results in the basis in which M was originally defined, cf.Equation (12).Let us denote this basis by {|x , |y , |z }.The eigenvectors |n ± and |n 0 of M are easily obtained from those of X 3 , cf.Equation (12).The eigenvectors of X 3 are, in turn, analogous to those of Pauli's σ y , namely |± = (|x ∓ i|y )/ √ 2, plus a third eigenvector that is orthogonal to the former ones, that is, |0 = |z .In order to obtain the eigenvectors of n • X, with n = (sin θ cos φ, sin θ sin φ, cos θ), we apply the rotation exp(φX 3 ) exp(θX 2 ) to the eigenvectors |± and |0 , thereby getting |n ± and |n 0 , respectively.All these calculations are easily performed using the CH-method.
Once we have |n ± and |n 0 , we also have the transformation matrix T that brings M into diagonal form: T −1 M T = M D = diag(−i, 0, i).Indeed, T 's columns are just |n − , |n 0 and |n + .After we have carried out all calculations in the eigenbasis of M , by applying T we can express the final result in the basis {|x , |y , |z }, thereby obtaining the desired expressions in three-vector notation.Let us illustrate this procedure by addressing the evolution equation In matrix form, such an equation reads The solution is given by S(t) = exp(At)S 0 .The eigenbasis of A is the same as that of M .We have thus The projectors |n k n k | can be written in terms of the powers of A by solving the system Finally, we can write the solution S(t) = exp(At)S 0 in the original basis {|x , |y , |z }, something that in this case amounts to writing M S 0 = n × S 0 and M 2 S 0 = n (n • S 0 ) − S 0 .Equation (39) was also addressed in [20], but making use of the operator method.The solution was given in terms of a series expansion for the evolution operator.In order to write this solution in closed form, it is necessary to introduce sinand cos-like functions [20].These functions are defined as infinite series involving two-variable Hermite polynomials.The final expression reads like Equation ( 11), but with sin and cos replaced by the aforementioned functions containing two-variable Hermite polynomials.Now, one can hardly unravel from such an expression the physical features that characterize the system's dynamics.
On the other hand, a solution given as in Equation ( 45) clearly shows such dynamics, in particular the damping effect stemming from the λ-term in Equation (39), for λ > 0. Indeed, Equation (45) clearly shows that the state vector S(t) = exp(At)S 0 asymptotically aligns with Ω while performing a damped Larmor precession about the latter.The case ∂Ω/∂t = 0 is more involved and generally requires resorting to Dyson-like series expansions, e.g., time-ordered exponential integrations.While this subject lies beyond the scope of the present work, it should be mentioned that the CH-method can be advantageously applied also in this context.For instance, time-ordered exponential integrations involving operators of the form A + B(t) do require the evaluation of exp A. Likewise, disentangling techniques make repeated use of matrix exponentials of single operators [21].In all these cases, the CH-method offers a possible shortcut.

The Jaynes-Cummings Hamiltonian
We address now a system composed by a two-level atom and a quantized (monochromatic) electromagnetic field.Under the dipole and the rotating-wave approximations, the Hamiltonian of this system reads (in standard notation) where δ = ω 0 − ω.A standard way [22] to calculate the evolution operator U = exp(−iHt/ ) goes as follows.One first writes the Hamiltonian in the form H = H In this way, one can simulate Lorentz transformations in the optical laboratory [23].The above formalism is of great help for designing the corresponding experimental setup.

Conclusions
The method presented in this paper-referred to as the Cayley-Hamilton method-proves advantageous for calculating closed-form expressions of analytic functions f (A) of an n×n matrix A, in particular matrix exponentials.The matrix A is assumed to be a diagonalizable one, even though only its eigenvalues are needed, not its eigenvectors.We have recovered some known results from classical and quantum mechanics, including Lorentz transformations, by performing the straightforward calculations that the method prescribes.In most cases, the problem at hand was reshaped so as to solve it by dealing with two-by-two matrices only.
46)Let us denote the upper and lower states of the two-level atom by |a and |b , respectively, and the Fock states of the photon-field by |n .The Hilbert space of the atom-field system is spanned by the basis B = {|a, n , |b, n , n = 0, 1, . ..}.The states |a, n and |b, n are eigenstates of the unperturbed Hamiltonian H 0 = ω 0 σ z /2 + ωa † a.The interaction Hamiltonian V = g a † σ − + aσ + couples the states |a, n and |b, n + 1 alone.Hence, H can be split into a sum: H = n H n , with each H n acting on the subspace Span{|a, n , |b, n + 1 }.Within such a subspace, H n is represented by the 2 × 2 matrix