Analytical evaluation of the numerical coefficients in the Zassenhaus product formula and its applications to quantum and statistical mechanics

This paper studies the exponential of the sum of two non-commuting operators as an infinite product of exponential operators involving repeated commutators of increasing order. It will be shown how to determine two coefficients in front of the nested commutators in the Zassenhaus formula. The knowledge of one coefficient is enough to generate a closed formula that has several applications in solving problems ranging from linear differential equations, quantum mechanics to non-linear differential equations.


Introduction
A crucial problem of mathematical physics is the development of the exponential of the sum of two non-commuting operators. A very deep and important application of the Baker-Campbell-Hausdorff formula [1,2,3] and its dual, the Zassenhaus formula [4,5,6], is related to quantum mechanics. Quantum mechanics is the kingdom of the operators and the Baker-Campbell-Hausdorff formula, jointly with the Zassenhaus formula, are an important mathematical tool. Even the simplest case of a timeindependent one-dimensional Hamiltonian shows the difficulty of the problem. If considered by a general point of view, both the time evolution of a state |ψ(t) and the eigenvalues problem associated with it are formidable problems. Indeed the following well-known Schrödinger's equations where V (x) is a generic potential, are an unsolved problem. Exact solutions are known for a particular choice of the potential function only. Schrödinger's equations is a particular case of a linear differential equation of order n and with variable coefficients y (n) (z) + an−1(z)y (n−1) (z) + · · · + a0(z)y(z) = 0, where y (j) means the jth derivative of y(z). It can be rewritten as d dz Y (z) = Y (z)M (z), Y (z) ≡ y(z) · · · , y (n−1) (z) .
Reducing Eq. (3) to a linear differential equation for a linear operator, Eq. (4), then we may write a formal solution via the Magnus expansion [4]. This last case is directly related to the solution of the problem of a timedependent Hamiltonian that is an intensive research field [7,8,9,10,11,12,13,14,15,16]. In particular in Refs. [13,15] the investigation is performed using the Zassenhaus formula. In the present paper, we will find an approximated formula obtained from the Zassenhaus expansion based on the knowledge of the coefficients in front of the nested commutators. Another important active research field is Statistical Mechanics. We will focus on the liouvillian approach to statistical mechanics, i.e. ∂ ∂t ρ(q1 · · · qn, t) = − n i=0 ∂ ∂qi [ai(q1 · · · qn, t)ρ(x, t)].
The associated Langevin equations (Ref. [17] for an historical point of view and [18] for a generalization) are dqi dt = ai(q1 · · · qn, t), i = 1, · · · n, where ai(q1 · · · qn, t) are given functions that in general can be stochastic. Its connection to statistical mechanics can be found via the Van Kampen's lemma [19] that relates the liouvillian density to the probability density through the relation P (q, t) = ρ(q, t) . It is virtually impossible to list a complete bibliography on the topic. For this reason, the reader is referred to few exemplary textbooks or paper collections [20,21,22,23,24]. As we may infer by direct inspection, Eqs. (6) are a set of non-linear equations whose solution is generally unknown. Using the results of the paper we will find an approximate analytical solution of Eqs. (6) and consequently of Eq. (5). To realize that we will use the Liouville equation or more in general, a linear partial differential equation. Solving nonlinear problems through the solution of linear problems is a technique used in literature (see Ref. [25] on this topic) and a recent method based on spectral decomposition is presented in Ref. [26]. Applying the results of the paper we will be able to evaluate the argument of the liouvillian density and through the method of characteristics, we will give a formula for the solution of Eqs. (6).
The paper is organized as follows: In Sec. 2 we will obtain an analytical expression for the coefficients of the nested commutators [A, [A, · · · , [A, B]]] and [[[A, B], · · · , B, B]. In Sec. 3, using the result of the previous section, we will find an approximated closed formula for linear differential equations. In Sec. 4 we will apply the previous results to recover the fundamental energy level of a quantum system in a generic shallow potential well and we will shortly sketch a solution for the time evolution of a quantum state. In Sec. 5, using the result of Sec. 2, we will find an approximated closed formula for the first order non-linear differential equations applied to statistical mechanics.

The Zassenhaus product formula
In this section, we will study an analytical approach to determine the coefficients in the Zassenhaus exponents in front of the nested commutators. The knowledge of these coefficients inspired several recent works (see for example [28,29,30]). As we will see, the result will have several important applications. Our starting point is the analytical solution of the operatorial differential equation obtained in Ref. [27]. Given the equation its solution is where M (t) is a time dependent operator, such as for example a n × n matrix, and O is defined as The symbol I is the identity matrix and it is understood that O acts on the left, i.e. on M (t0). The symbol exp [uO] M (t0) |t 0 =0 means that it must be evaluated at t0 = 0. Equivalently for the equation we have where it is understood that O acts on the right. The aim of this section is to rigorously show that the coefficient in the Zassenhaus product formula associated to the term [A, As starting point, we will prove that when the matrix M (t) and its derivatives commute at any value of the parameter t than we will recover the solution of the differential equation (7). To fix the ideas let us assume that M (t) is a diagonal matrix. Using the Trotter product formula [31] we may rewrite solution (11) as Note that the operator exp u N ∂t 0 translates of a quantity u N the argument of M (t0). Without loss of generality we set t0 = 0. Performing the translation and exploiting the commutating properties of the matrix M we may write Since M j u N commutes for any value of j u N , we rewrite Eq. (13) as We now use the fact that for N → ∞ we can write the sum in the continuous limit, i.e.
Plugging the result into Eq. (14) we finally obtain We will use this result to determine the coefficient in front of the term and setting A = ∂t and B = M (t) we have that Formally solution of Eq. (11) can be rewritten as Since M (t) and its derivatives are commutating matrices ∀t, all the commutators present in the formula vanish except for those that contain only once the matrix M . Using the properties of the translational operator exp [u∂t] we may write the following Exploiting again the commutating properties of M (t), the exponential matrix product is commutative and we may write Eq. (20) as This result has to be the same result obtained in Eq. (16). This is so if and only if Integrating by part repeatedly the left side of Eq.
Equating the coefficients of the power laws of the two series we infer that We may use the same development in terms of nested commutators identifying A with M and B with ∂t. In this case Eq. (19) may be written as As before the solution has to be the solution given in Eq. (16). But in this case the identification of the correspondent coefficient is not so straightforward. We assume the same conditions on the matrix M as before but we consider that its elements are polynomial. Without loss of generality, we set its matrix element as Since the matrices commute, we may perform the sum of the exponentials exp [cm [[M (t), ∂t] · · · ∂t]] and, according to Eq. (22), we must have For each power of the polynomial we determine the coefficientcm. We may write the following recursive relationship It is straightforward to show that Eq. (27) correspond tō Formulas (24) and (28) give an analytical expression for the coefficients in respectively. The exactness of those formulas can be checked in Ref. [28,29,30].

Ordinary differential equations with small variable coefficients
In this section, we will use the result of Sec. 2 to provide an approximate solution of a differential equation in the form where ε is a small parameter and an−1(t, ε) → 0 for ε → 0. The above equation can be written as a first order differential equation in the form of Eq. (10). We will consider the problem set in Eq. (29) when the matrix M may be written as M (t) = εN (t), where N (t) is a matrix with finite values of its elements. We assume that it holds the following equality where k is a given integer. Using the results of the previous section we obtain at ε order where the coefficients cn are given by Eq. (24). In particular if M (t) commutates starting from the second derivative on, we have the following approximate formula As example we consider the following second order differential equation Note that the small coefficient ε can be eliminated by rescaling the variable t and leading to an alternative formulation of the problem The above equations can be generated by a two dimensional system and, applying Eq. (32), we obtain where It is no hard to see that a sufficient condition so that Eqs. (36) and (37) are a valid approximation in where K is a finite positive number. For an infinite interval, t ∈ [0, ∞), a sufficient condition is given by f (t) = k + g(t) where k is a finite value parameter and g(t) is a finite value function such that g(t) → 0 faster than 1/t for t → ∞. Less strict conditions can be found according to the case under study. As example, we consider the differential equation The coefficient f (t) is on purpose taken with an exponentially decaying term. This will show that the contribution of such a term to the correction to the unperturbed solution (x0(t), x ′ 0 (t)) = (cos εt, −ε sin εt) is not negligible. The approximate solution given by Eq. (36) is visually indistinguishable and it is more convenient to plot the percent error. Since the solution is oscillating, to plot the percent error we can not use the traditional definition | xex(t) − xapp(t) | / | xex(t) |. This because when the exact solution xex(t) vanishes, the error would be divergent, regardless of how the approximate value is near to the exact one. To overcome this difficulty we first define the percent error as the "distance" between two functions with respect one of the two. In mathematical symbols where y2(t) is the reference function, y1(t) is the approximated function and · L 1 is defined as In Fig. 1 it has been plotted the percent error between the exact solution and the "unperturbed" solution i.e. x0(t) = cos εt. The percent error, given by xex(t) − x0(t) L 1 / xex(t) L 1 100, is of the order of 10%. In Fig. 2 it has been plotted the percent error between the exact solution and x(t) given by Eq. (36). The percent error is of the order of 2%. We end the section noting that if we consider the derivative of the solution then the percent error between the exact solution and the unperturbed  solution x ′ 0 (t) = −ε sin εt ranges from 10% up to 50%. Taking the derivative of Eq. (36), the percent error between the exact solution and the approximate solution is of the order of 1%. Finally, we note that one can build more accurate solutions rearranging x(t) the derivative of the function x(t) and y(t) according to the initial values (x0, y0). A detailed study of this topic is out of the scope of this paper.

Applications to quantum mechanics: shallow potential well and time evolution of a quantum state
Many works in quantum mechanics are devoted to finding approximate solutions of Schrödinger's equation. Among them, several methods use Zassenhaus expansion formula [13,15]. As already noticed in the previous sections, this paper is more focused on the wide range of possible applications of the proposed method rather than its accuracy. In spite of the fact that we are working with a first order approximation, the accuracy is enough to find both the energy level of a quantum particle in a shallow potential well and the time evolution of a quantum state. First we will study the stationary Schrödinger's equation to evaluate the energy of a particle in a shallow potential well. The one dimension Schrödinger's equation reads as where for shortness we dropped off the argument of the wave function. Without loss of generality, we assume that the potential is of the form V (x) = −V0g(x/a) where a is a length scale parameter and g(z) a positive function. Performing the variable change z = x/a and defining e = E/V0 we may rewrite Eq. (43) as We will focus on a shallow potential well. The literature on shallow potential well is vast and we limit ourselves to Ref. [32] as textbook. In particular, we consider a symmetric potential with a finite value for the minimum as shown in Fig. 3. The condition for shallow potential well is  (36) and (37), we obtain the necessary conditions to have a decaying solution at infinity. At the zero and first order on g(z), we have Vanishing the determinant of the above system we have the condition or, in terms of the energy E, We recovered the result given at pag. 163 in Ref. [32]. We end this section shortly studying the time evolution of a quantum state. To apply directly the formalism of Sec. 2 we consider the time evolution of ψ(t)|. Setting for brevityh = 1, we have where H0 is a time-independent Hamiltonian and εH1(t) is a small timedependent correction to H0. Passing to the interaction representation [33] we may always transform Eq. (50) into with Using the results of Sec. 2 at ε order we may write the approximate solution where the coefficient cn are given by (24). Without to specify the commutator [H

Statistical mechanics and non-linear equations
Generally speaking, a non-linear first order differential equation is an unsolved problem. There is not a general method to solve such a differential equation. One can use a recursive method such as the method of successive approximations (see for example Ref. [34]). These kind of methods are quite hard to handle from an analytical point of view. In this section, we will provide an approximate closed formula for first order non-linear differential equations. Let us consider the non-linear first order equation where x is a vector and a(x, t) is a given vectorial function. Its connection to statistical mechanics can be found simply interpreting Eq. (54) as Langevin equation [17,18]. The corresponding Liouville equation is If a(x, t) is a stochastic function then, via the Van Kampen's lemma [19], the liouvillian density is related to the probability density by the relation P (x, t) = ρ(x, t) . The difficulty to solve Eq. (54) is inherited by Eq. (55) since also for Eq. (55) there is not a general method for solving it. But, using physical arguments, if we interpret Eq. (54) as the motion equation of a particle, then the solution of Eq. (55) is where δ(z) is the Dirac delta function and xP (t) is the solution of Eq. (54). An important advantage of considering Liouville equation, or, more in general, a linear partial differential equation, is that for the partial differential equation we may use a linear theory. For our purposes, we can exploit the theory developed in Sec. 3. For sake of simplicity, we shall focus on a onedimensional case with a(x, t) a real function. The one-dimensional case can be related to anomalous diffusion. In the past thirty years, anomalous diffusion has been intensively investigated [35,36,37,38,39]. We will study the following non-linear equation [40,41] where, in a over-damped regime, V (x) represent a potential and F (t) can be interpreted as a force and more in general as a stochastic force. The above equation can be always rewritten as where now t, x, v(x) and f (t) are dimensionless quantities and ε is a parameter that we will take small. Alternatively, making the variable change t → τ /ε, we may consider also equations like The corresponding Liouville equation is Due to the formalism developed in Sec. 3 it is more convenient to study the equivalent equation with the understanding that the partial derivative with respect to x acts on the left (stressed by the arrow on top of the partial derivative symbol), while the partial derivative with respect to t acts on the right. Applying Eq. (32) we obtain the approximate solution To perform the action of the exponential operators we make a change of variable as following and we have where we used the definition of ∆(u) given in Eq. (38). We may rewrite Eq. (63) as With enough accuracy we may substitute the partial derivative ∂w with ∂u.
Integrating by part, taking in account thatΦ [w(x − ε∆(t)) + t] |t=0 = Φ(w(x)) = Φ(x, 0) we finally obtain The solution of the linear problem is giving us the solution of the nonlinear problem (58) setting as constant the argument of functionΦ(z), i.e.
The constant can be determined using the initial condition x = x0 at t = 0. It is worthy to stress the following points: a) When f (t) is a constant Eq. (68) gives the known exact solution.
b) The difficulties to invert Eq. (68) with respect to x are the same of the exact case, i.e. when f (t) = constant.
c) The derivative of the potential function v ′ (x) is subjected to the less restrictive condition |v ′ (x)| < K with K a positive constant, while f (t) must satisfy the more strict condition (39).
We now will consider an example that will illustrate the points (a), (b), (c). We will study Eq. (68) using as potential a periodic potential [that does not satisfy (39)] in presence of an arbitrary force f (t). This problem is known as the diffusion in the egg-carton potential [42,43]. The inversion of Eq. (70) is straightforward and x(t) can be found explicitly. This step is left to the reader. In Fig. 4 we compare the numerical and the analytical solution in the case of a slowly damped harmonic force f (t) f (t) = f0 + b cos(Ωt) The agreement is excellent and the percent error evaluated with formula (41) is smaller than 1%. A detailed study of the solution is out of the scope of this paper.

Conclusion
In this paper, it has been shown how to determine two coefficients in the Zassenhaus expansion. The analytical expression of these coefficients allowed a series of applications ranging from nth order linear equations, evaluation of eigenvalues, to first order non-linear equations. We have been able to build an approximate closed formula at first order of the expansion parameter ε in the case of linear differential equations and non-linear first order differential equation. We showed the applications regarding the evaluation of eigenvalues and the statistical mechanics. Finally, we gave a formula implying an infinite matrix product of exponential operators for the time evolution of a quantum state. According to the expression of the Hamiltonian, the product of the of exponential operators can be summed and expressed as a closed analytical formula. The presented approach could be extended to evaluate other coefficients in the Zassenhaus expansion and to solve higher order non-linear differential equations. This is left for a future work.