The flow method for the Baker-Campbell-Hausdorff formula: exact results

Leveraging techniques from the literature on geometric numerical integration, we propose a new general method to compute exact expressions for the BCH formula. In its utmost generality, the method consists in embedding the Lie algebra of interest into a subalgebra of the algebra of vector fields on some manifold by means of an isomorphism, so that the BCH formula for two elements of the original algebra can be recovered from the composition of the flows of the corresponding vector fields. For this reason we call our method the flow method. Clearly, this method has great advantage in cases where the flows can be computed analytically. We illustrate its usefulness on some benchmark examples where it can be applied directly, and discuss some possible extensions for cases where an exact expression cannot be obtained.

Leveraging techniques from the literature on geometric numerical integration, we propose a new general method to compute exact expressions for the BCH formula. In its utmost generality, the method consists in embedding the Lie algebra of interest into a subalgebra of the algebra of vector fields on some manifold by means of an isomorphism, so that the BCH formula for two elements of the original algebra can be recovered from the composition of the flows of the corresponding vector fields. For this reason we call our method the flow method. Clearly, this method has great advantage in cases where the flows can be computed analytically. We illustrate its usefulness on some benchmark examples where it can be applied directly, and discuss some possible extensions for cases where an exact expression cannot be obtained.

Introduction
The Baker-Campbell-Hausdorff (BCH) formula is a well known result in the theory of Lie groups and Lie algebras [1,19,26] that found wide application in both mathematics and physics: from Lie theory [26] to numerical integration [14,37], quantum mechanics [5,36] and statistical mechanics [39]. Roughly speaking, it links the composition between two elements of a Lie group in a connected neighbourhood of the identity to an element of the Lie algebra whose exponential map is the group element which results from the composition. In other words, it provides a means to compute the map Z : g × g → g defined by Z(A, B) := log(e A e B ) , where e A , e B ∈ G are the corresponding elements in the associated Lie group given by the exponential map. When the Lie algebra is commutative, the BCH formula is simply However, Lie algebras are in general not commutative and finding closed-form expressions or even accurate estimates for Z can become a daunting task. The general BCH formula can be expressed (formally) in terms of the integral formulation where ad A is the adjoint representation of A and I is the identity matrix. To approximate (2), the most common method is to employ the series expansion For a detailed treatment of the topic, we refer the reader to [19,26,42]. The problem of computing a closed form for (2) for a given particular algebra has already been addressed in various settings, employing a number of different algebraic and analytic techniques [6,12,19,31,32,40,45]. In this paper, we present a dynamic approach to such problem: we compute closed forms of the BCH formula by exploiting an isomorphism between the given Lie algebra and an appropriate subalgebra of the infinitesimal contactomorphisms on a contact manifold. The latter in turn is also isomorphic to the Lie algebra of the Hamiltonian functions endowed with the Jacobi bracket that naturally arises from the contact structure [29]. This helps us to simplify the calculations. The mapping between functions and vector fields is always an isomorphism for contact Hamiltonians [29], differently from the symplectic case [18]. Therefore, even in the case of some algebras which arise naturally in the symplectic setting (e.g. the Heisenberg group), the use of the contact structure is more natural, as it guarantees a one-to-one correspondence between contact Hamiltonian functions and contact Hamiltonian vector fields.
We prove the usefulness of our method using some interesting examples arising from the recent literature in contact geometry and Hamiltonian systems, namely, the Heisenberg algebra [26], the contact Heisenberg algebra [12], the quadratic symplectic and the quadratic contact algebras, and the complexification of su (2). Furthermore, we also include an analysis of contact splitting numerical integrators [14,48] for the damped harmonic oscillator [11]. More precisely, in this example our method will be used in order to obtain a priori estimates of the numerical errors for different splitting integrators.
The paper is organized as follows: in Section 1 we briefly review some concepts in contact geometry and contact Hamiltonian mechanics. In Section 2 we describe the flow method. Section 3 is devoted to the aforementioned examples, presented in detail, while Section 4 presents the conclusions and some additional directions for further research. Additional results, including a matrix representation of the algebras treated in the paper, are presented at the end of the paper in Appendix A.

Contact geometry
We summarize here the essential concepts of contact geometry that are needed for this work and refer the reader to [4,7,10,11,13,23,24,29,33] for additional details, including the relevant proofs. The main defining object for us is a special differential 1-form, called the contact form, which induces both a volume form on the manifold and a maximally non-integrable distribution on its tangent bundle, called the contact structure.
In this case η is called a contact form and its kernel at each point defines a maximally non-integrable distribution of hyperplanes D, called the contact structure.
important to stress that the definition of the horizontal distribution is invariant under the transformation η → f η, since ker(f η) = ker(η), while the vertical distribution is not [33]. The vertical distribution is 1-dimensional and is described in terms of the so-called Reeb vector field, which is unique once a contact form is fixed. Proposition 1.1. Given a contact manifold (M, η) there exists a unique vector field R, the Reeb vector field, such that η(R) = 1 and dη(R, ·) = 0 .
Diffeomorphisms that preserve the contact structure are of particular interest. We have seen that two contact forms define the same contact structure if they are equal up to multiplication by a non-vanishing function: this allows to classify contactomorphisms into two classes.
Similarly, we can also classify the corresponding infinitesimal transformations.
Finally, we want to give a coordinate description of a contact manifold. In Definition 1.1 we have seen that the distribution is described by a 1-form. The following theorem provides local canonical coordinates. Theorem 1.1 (Darboux's theorem). Consider a contact manifold (M, η) and a point x ∈ M . Then there exist local coordinates (q i , p i , s) in a neighbourhood of x such that the contact form is written as where here and in the following Einstein's summation convention over repeated indices will be always assumed.
The coordinates in Darboux's theorem are called canonical coordinates. Note that in these coordinates we have dη = dq i ∧ dp i and R = ∂ ∂s .

Contact hamiltonian flows and Jacobi structures
The contact form allows to define an isomorphism between functions and vector fields. We call these vector fields contact gradients or contact Hamiltonian vector fields [33].
We denote the map that assigns to a function H its associated contact Hamiltonian vector field as a η , so that X H = a η (H) and the function H is called contact Hamiltonian function.
Any contact Hamiltonian vector field is an infinitesimal contactomorphism, as shown by a direct application of Cartan's magic formula: In particular, it is strict if and only if R(H) = 0. Remarkably, one can prove also the opposite, i.e. any infinitesimal contactomorphism is a contact Hamiltonian vector field, with the corresponding (contact) Hamiltonian function being recovered by using the first condition in (3).
In canonical coordinates (q i , p i , s), the contact Hamiltonian vector field is given by [10,11] and its integral curves satisfy the contact Hamiltonian equations of motion In what follows we will denote the contact Hamiltonian flow of X H , or equivalently the solution of the above system, by where x 0 ∈ M is the initial condition. Differently from its symplectic counterpart, a contact Hamiltonian function is not conserved along the flow of X H . Instead, its dynamics depends on R(H) as can be noted in the following This means, in particular, that the surface defined by H = 0 is invariant under the contact Hamiltonian evolution and that this is generally not true for the other level sets [10,11,30].
Classical Hamiltonian systems can be studied in terms of the algebra of Hamiltonian vector fields on a symplectic manifold and in terms of the algebra of Hamiltonian functions with the Poisson bracket. This duality persists also in contact Hamiltonian mechanics. The set of contact Hamiltonian functions on a contact manifold (M, η) has a structure of a (local) Lie algebra in the sense of Kirillov [27]: the bracket can be defined by [14] {f, g} η : which, in local canonical coordinates, has the form This bracket is bilinear, antisymmetric and satisfies the Jacobi identity, but not the Leibniz rule. Instead, one has {f, gh} η = g{f, h} η + h{f, g} η + gh{f, 1} η .
For our purposes, a fundamental difference between the symplectic and the contact settings is the fact that in the contact case the map a η , which assigns to each function the corresponding Hamiltonian vector field (see Definition 1.4), is an isomorphism of Lie algebras, while in the symplectic case it is known to be just a homomorphism, as constant functions all belong to the kernel [29]. This is the content of the next result [2,29].
is an isomorphism between the algebra of functions on M with the Jacobi bracket and the algebra X C (M ) of infinitesimal contactomorphisms with the standard Lie bracket. Its inverse is provided by the contraction with the contact form: η : Remark 1.1. We can compare this with the usual Hamiltonian structures on symplectic manifolds. A symplectic manifold is a pair (M, ω), where M is an even-dimensional manifold and ω is a closed, non-degenerate, 2-form on M called the symplectic form [18,33]. The symplectic form allows to associate to each function H : M → R a vector field X H , called the Hamiltonian vector field or the symplectic gradient, by ω(X H , ·) = −dH.
In canonical coordinates (q i , p i ), the symplectic form is ω = dp i ∧ dq i and thus the integral curves of X H give Hamilton's equations in their standard form Finally, the Poisson bracket [33] of two functions is defined as The side-to-side comparison of the symplectic and contact structures is summarized in Table 1.

The flow method
In this section we present the main idea of this work, proceeding in order of decreasing generality and increasing applicability. Given a Lie algebra g and two elements A, B ∈ g, in the following definition we outline a method to compute (1) without recurring to the general BCH formula (2). Definition 2.1. Given a Lie algebra g, the general flow method consists of the following:

Symplectic Hamiltonian Systems Contact Hamiltonian Systems
1. Find a Lie algebra isomorphism φ between g and some subalgebra g M of X(M ) for some manifold M ; 2. Given A, B ∈ g, replace the product e A e B by the composition of the flows A graphical representation of the idea behind steps 2. and 3. of the flow method is depicted in Fig. 1. Remark 2.1. Clearly the flow method, whenever applicable, is a way to compute (1) without recurring to the BCH formula. An advantage of the flow method is that in some cases computing the composition of the flows e tX A • e tX B can be very easy, whereas computing the corresponding product e A e B can be rather cumbersome. More importantly, once we have the composed flow e tX A •e tX B , then it is a standard calculation to recover the infinitesimal generator Z(X A , X B ). Therefore, step 3. yields a substantial simplification. However, the main difficulty in the process is to find M and g M such that (i) there exists an isomorphism between g and g M and (ii) the flows e tX A and e tX B can be explicitly found. In general this is not even possible. However, as we show below, in some special cases one can work out the flow method and thus obtain closed-form expressions for the corresponding BCH formula.
Remark 2.2. The flow method was already used in a similar fashion in [16] in the context of finding the modified Hamiltonian for different symplectic numerical methods for the harmonic oscillator, as we will see in Example 1. However, the usefulness of this method has not been established before, perhaps because, as we will show shortly, the symplectic context is not the most useful one for the general purpose at hand.
To implement the flow method, it turns out to be particularly convenient to choose M to be either a symplectic or a contact manifold and g M to be a subalgebra of the Hamiltonian vector fields on M , which we will denote by h M . To avoid restating this all the time, we will denote the flow method with this particular choice of manifold and subalgebra the Hamiltonian flow method, overloading the word Hamiltonian to denote both the contact and symplectic choices.
3. An obvious advantage of the Hamiltonian flow method is the fact that in this case we are provided with an additional Lie algebra homomorphism ψ between h M and the corresponding Hamiltonian functions. It is much easier to specify the desired Lie algebras in terms of Hamiltonian functions rather than in terms of vector fields directly. For instance, the Heisenberg Lie algebra is often written as the algebra of linear functions on a symplectic manifold together with their Poisson bracket, but rarely expressed as the corresponding Lie algebra of vector fields (and there is a reason for it, since the two in this case are not isomorphic, see below). Now, let {H 1 , . . . , H n } denote a basis for ψ(h M ), where ψ is the Lie algebra homomorphism mentioned in Remark 2.3 above. Then a general element in ψ(h M ) takes the form This means that, given 2 elements we need to find the coefficients z i as functions of a i and b i . At this point, we need to add the following assumption.
ii. the flow e tX H of the general H (with unspecified coefficients h i ) can be obtained by quadratures.
If (and only if) Assumption 2.1 is satisfied, then points 3. and 4. in the general flow method (Definition 2.1) can be replaced bỹ 3a. equate e tX Z(A,B) = e tX A • e tX B to obtain a system of algebraic equations relating the known coefficients a i and b i to the unknown coefficients z ĩ 3b. solve the above system to find the coefficients z i as functions of a i and b i , say z i (a, b), where ψ is the isomorphism in Assumption 2.1 and φ the one from step 1. of the flow method.
Definition 2.2. We call the method given by the steps 1., 2.,3a.,3b.,4. the algebraic Hamiltonian flow method, or AHFM. The symplectic approach to the AHFM will be specified by SAHFM, while the contact approach by CAHFM.
Let us see how the AHFM works in the example, taken from [16], of finding the modified Hamiltonian for a 1st order splitting integrator for a classical harmonic oscillator. The purpose of this example is to develop the intuition on the flow method. At the same time, it also provides a first practical application (in discriminating among different numerical integrators).
Example 1 (Splitting integrators for the harmonic oscillator). Let us consider the 1dimensional harmonic oscillator, with M = R 2 , ω = dp ∧ dq, and Hamiltonian function The symplectic Hamiltonian vector fields corresponding to these functions are From the theory of geometric numerical integration (see Section 3.6 or [14,16,46]), we know that is a (symplectic) numerical integrator of 1st order for the flow of H, meaning that, for small values of τ , the timestep, one has S 1 (τ )(x 0 ) = exp{τ X H }(x 0 ) + O(τ ) for any initial condition x 0 . We also know that, being symplectic, S 1 (τ ) is Hamiltonian with respect to some other Hamiltonian functionH(q, p; τ ), called the modified Hamiltonian (see also e.g. [25,28] for more details). Finally, to findH(q, p; τ ) we need to use the BCH formula (2) with A = X T and B = X V . Instead, in the following we will apply the flow method.
We begin by observing that the commutation relations of our functions of interest with respect to the Poisson bracket { · , · } P B , see Remark 1.1, are , qp} provides a set of generators of a Lie subalgebra of ( C ∞ (M, R), { · , · } P B ). Thus we can argue that the modified HamiltonianH(q, p; τ ) must be (for any τ ) a linear combination of the elements of the basis of this Lie subalgebra, that is, Now, to compute the unknown coefficients a(τ ), b(τ ), c(τ ) we observe that, for a fixed timestep τ , we can solve the differential equation The solution, evaluated at t = τ , is The map , the numerical integrator map. Here we have used ω = −c 2 +4ab. We know from (6) that the same function is given by the composition of the flows of X T and X V at t = τ , i.e. the maps We can now compute a(τ ), b(τ ), c(τ ) by comparing which should hold for any arbitrary initial condition (q 0 , p 0 ). Solving the algebraic equations (8), we obtain which in turn give the modified Hamiltonian (7). We notice that knowledge of the modified Hamiltonians for different integrators can help identify the most suitable one for the particular problem, as we will explain in more detail in Section 3.6.
The purpose of the next example is to show that the flow method really relies on the contact algebra and cannot, in general, be used just by employing symplectic Hamiltonians.
Example 2 (The Heisenberg algebra: SAHFM vs CAHFM). Consider the Heisenberg algebra H. It is well-known [38] that it can be represented as the Lie algebra of linear functions on the 2-dimensional symplectic manifold (R 2 , ω = dq ∧ dp), that is, where { · , · } P B denotes the standard Poisson bracket. The only non-vanishing commutator is { p , q } P B = 1, as it should be. However, in this case there is no isomorphism between H and a subalgebra of h M of symplectic Hamiltonian vector fields, as required by Definition 2.1. This is because the symplectic gradient has a non-trivial kernel, consisting of all the constant functions [18]. A direct way to see this is by looking at Hamilton's equations in the symplectic case (equations (5) in Remark 1.1). Clearly, two Hamiltonians differing by an additive constant generate the same system of equations, and thus the same vector field. Therefore, in this example (and in the general case), step4. of the SAHFM is problematic. Notice that this problem does not appear if we use the CAHFM. This is because in general the contact gradient defines an isomorphism, as we have seen in Theorem 1.2. Also in this case, looking back at the contact Hamiltonian equations can help clarify this statement. Indeed, one can verify from equations (4) that in this case two Hamiltonians differing by an additive constant generate different vector fields.

Closed-form expressions for the BCH formula for various algebras
In this section we show how the flow method can provide explicit closed forms of the BCH formula in a number of examples, namely, the classical Heisenberg algebra, the contact Heisenberg algebra introduced in [12], the quadratic symplectic algebra, and the quadratic contact algebra. Finally, we use the latter to compute the modified Hamiltonians of several possible 1st order contact splitting integrators for the damped harmonic oscillator [14].

The Heisenberg algebra
so that its associated contact Hamiltonian vector field is The flow of this vector field for an initial condition (q 0 , p 0 , s 0 ) is explicitly computed as Our aim is to compose the flows of two different contact Hamiltonians h = aq +bp+z1 andh =āq +bp +z1 and exploit the fact that such flow is a contact flow of some contact Hamiltonian h belonging to the same algebra. Ultimately, this means that we are left to compute the explicit composition of the flows of h andh and equate the coefficients with those of the flow of a general h = αq + βp + ζ1 in order to find the appropriate h. In practice, this amounts to solving a system of algebraic equations obtained as follows: the composition of the flows of h andh at some time t, namely e tXh • e tX h , is computed explicitly as 2 t 2 (ab +ā(2b +b)) − t(q 0 (a +ā) + z +z) + s 0 .
. By equating this with the flow of h, namely at t = 1, we obtain the following defining relations for the coefficients which indeed is the standard closed-form expression for the BCH formula for the Heisenberg algebra. In Figure 2 we give a graphical representation of the method.

The contact Heisenberg algebra
The contact Heisenberg algebra (CHA) is a natural extension of the classical Heisenberg algebra [12] in the contact phase space (R 3 , η = ds − pdq), being defined as Therefore, any contact Hamiltonian function in this algebra has the following form: and the corresponding contact Hamiltonian vector field is Again, we can compute the flow explicitly as and thus consider the composition of the flows of two contact Hamiltonians f = aq + bp + cs + z1 andf =āq +bp +cs +z1; and compare it with the flow generated by f = αq + βp + γs + ζ1 at t = 1. This leads to the following solution which again is the correct closed-form expression for the BCH formula for this algebra (cf. [12]). The method is represented graphically in Figure 3.

The quadratic symplectic algebra
The quadratic symplectic algebra is the algebra in the two-dimensional symplectic manifold (R 2 , ω = dp ∧ dq) generated by the three quadratic functions, that is, Notice that this case includes Example 1 and that here we do not have constant functions as generators of the algebra, differently from Example 2. Therefore, there is an isomorphism between this Lie algebra and the algebra of the corresponding symplectic Hamiltonian vector fields. Indeed, to each f (q, p) = aq 2 + bp 2 + cqp, a, b, c ∈ R, we can associate uniquely Hamilton's equations for X f are given by the linear ordinary differential equation whose general solution at t = 1 for an initial condition x(0) = (q 0 , p 0 ) can be expressed as q(1) p(1) = cosh(det(A))I + sinh(det(A)) det (A) A q 0 p 0 .
Now we want to compose two of these flows, generated respectively by H A = a 1 q 2 +b 1 p 2 + c 1 qp and H B = a 2 q 2 + b 2 p 2 + c 2 qp, which correspond to the flows given by (10) with the matrices To keep the following expressions compact, it is convenient to define the notations Moreover, since A and B are Hamiltonian matrices, their product can be decomposed as follows AB = JS + λI , .
Substituting, and exploiting J −1 = −J, we can rewrite the right-hand side of (11) as Form the second equation we can obtain immediately the determinant of C. Moreover, using the inverse function of the hyperbolic cosine, we can also obtain the corresponding value of the hyperbolic sine by Now, using the second equation in (12), we obtain the final result (14) which is the closed-form expression for the BCH formula for the quadratic symplectic algebra (in matrix form). To our knowledge this result has not been reported elsewhere previously and therefore it is interesting in itself.
The graphical result of this procedure is plotted in Figure 4. Since we have been working on the two-dimensional manifold (R 2 , ω = dp ∧ dq), one may be surprised by the appearance of the (q, s) and (p, s) planes in Figure 4. Indeed, we are describing there the flow method applied to the extension of the same algebra to the contact manifold (R 2 × R, η = ds − pdq). In practice, the (p, q)-projection (leftmost plot) is exactly what we have computed here, and for the rest of the pictures we only need to extend (10) with the additional differential equationṡ = bp 2 − aq 2 . The relevance of this example will be clear in the next section on the quadratic contact algebra. There we will need to use a combination of the ideas presented here and in the previous examples.
This algebra is interesting since it contains the quadratic symplectic algebra, the conformal symplectic one (generated by {q 2 , p 2 , qp, s}), and the Hamiltonian function of the harmonic oscillator (with and without dissipation). We will come back to the latter in the next subsection, where we will rely on the results of this section to compute the modified Hamiltonian for different choices of splitting numerical integrators for such system. Hamiltonian functions in the algebra (15) have the form and their corresponding contact gradients are This system of differential equations is integrable by quadratures, since the first two equations are linear and the third one is linear in s once the first two are solved. In the same way as in the previous sections, we compose the flows of two contact Hamiltonians in the algebra, f = aq 2 + bp 2 + ds + cqp + z1 andf =āq 2 +bp 2 +ds +cqp +z1, and compare the result with the flow of a third Hamiltonian The flow method immediately provides an expression for the unknowns δ and ζ from the solution s(t), reading To obtain the remaining parameters α, β and γ we repeat the same computation as in the previous section: observe that the solution in the plane (q, p) is independent of the dynamics in s, so we only need to solve the dynamical system The solution of equation (18) at t = 1 can be expressed as whereÃ is the matrixÃ which is a Hamiltonian matrix, as in the previous section. When we compose the flows generated by the matrix A in (18)  This has to be compared with the flow generated bỹ which, at t = 1, has the same form as equation (19). Since we already know the equation for the parameter δ from (17), we can reduce this condition to (12) from the symplectic case: the exponential term exp(d/2 +d/2) simplifies and the remaining terms are the flows generated by Hamiltonian matrices of the form (20). This leads to the solution for C being given bỹ where x is defined as in (13) usingÃ andB. Once again this equation gives the closedform expression for the BCH formula for this algebra. Note that for z =z = 0 we recover a closed-form expression for the conformal symplectic algebra. To our knowledge this result has not been reported elsewhere previously and therefore it is interesting in itself. A graphical test of the results is shown in Figure 5.

The complexified su(2)
In this section we will show how to map the complexification of su(2), su(2) C ≃ sl(2, C) ≃ su(2)⊗C, to a complexified version of the quadratic symplectic Lie algebra of Section 3.3.
To keep the presentation short we will omit the details of the complexification procedure, for which we refer the readers to [26, Chapters 3.6 and 4.6]. The Lie algebra su(2) can be represented as the vector space spanned by the Pauli matrices σ 1 , σ 2 , σ 3 . The structure coefficients of the algebra are well known, and with our notation they correspond to [Σ k , Σ l ] = 2ϵ k l m Σ m , k, l, m = 1, . . . , 3.
The complexification of su(2) leads to the Lie algebra su(2) C ≃ su(2) ⊗ C spanned by These new generators satisfy the commutation relations [26,Chapter 4.6]: To find a mapping between the su(2) and the symplectic algebra, we start from the commutation relation It is a matter of explicit computation to find the mapping and its inverse With this mapping, an element of the su(2) algebra is mapped to an element of the complexified symplectic quadratic algebra by: To obtain the corresponding BCH formula, is then enough to consider the matrix and the respective matrix B, in equation (14) and then invert the mapping (22). This also provides an alternative derivation of the well-known parameterization of the composition law of the group SU (2) in terms of Pauli matrices. Now, using the second equation in (12), we obtain the final result

Splitting numerical integrators
The flow method can immediately be put to use in the computation of the modified Hamiltonian of certain splitting integrators. Contact splitting integrators are a family of numerical algorithms to integrate flows of Hamiltonian functions that can be written as sums of separate terms, i.e., It is convenient to assume that each h i is a Hamiltonian function whose flow can be explicitly integrated. Then, for a small timestep 0 < τ ≪ 1, the flow ϕ H x (τ ) of H can be approximated to first order in τ by the map where ⃝ n i=1 ψ i := ψ 1 • ψ 2 • · · · • ψ n , and to second order by The modified Hamiltonian of the numerical integrator is the Hamiltonian function that generates the integrator, that is, the Hamiltonian of the approximate flow (see e.g. [14] for more details). As previously mentioned, the 1-dimensional damped harmonic oscillator belongs to the quadratic contact algebra Span R {p 2 , q 2 , q p, s, 1} presented in the previous section. Indeed, its Hamiltonian is In the rest of this section we will focus on the six possible 1st order integrators for the Hamiltonian (23), that is, given by all the possible permutations of the three terms that form the Hamiltonian: the kinetic energy, the quadratic potential and the dissipation term. One can immediately check that their modified Hamiltonians must have the form H mod (q, p, s; τ ) = a(τ )q 2 + b(τ )p 2 + c(τ )qp + d(τ )s, and equation (17) implies that d(τ ) = γ. The remaining three parameters a(τ ), b(τ ) and c(τ ) can be found from (21). They are plotted in Figure 6 as functions of the timestep τ .
In order to provide a quantitative comparison of the different modified Hamiltonians, we exploit the natural bilinear form induced on a Lie algebra by the Killing form [19]: where Tr is the trace operator and ad X is the adjoint representation of the contact quadratic algebra given in Appendix A. Since b(X, Y ) in general fails to be positivedefinite, we use b(X, Y ) 2 as a (possibly degenerate) pseudo-distance on the space of contact Hamiltonian functions of the form (16). One can think of it as an analogue of the trace distance in quantum information theory and therefore we also refer to it as the trace distance with a slight abuse of terminology. In Figure 7 we plot the trace distance  (24). The rows are a(τ ), b(τ ) and c(τ ), respectively, while along the columns we vary the parameter γ ∈ 1 2 , 2, 4 .
between the original Hamiltonian and the modified Hamiltonian as a function of the time step τ . Plotting the parameters and the distance from the original Hamiltonian already provides a visual way to select better-performing integrators: in almost all cases the integrator S τ T CV (x) has the least distance from the original Hamiltonian (23), and as such it is expected to perform better. At the same time, we can observe that for small values of γ, S τ CV T (x) could become a slightly better choice. Of course, at this stage, this is still a heuristic, but already hints at the potential usefulness of the method in the analysis of the performance of certain numerical integrators.

Conclusions
In this paper we presented a remarkably simple method, the flow method, to compute closed-form expressions for the Baker-Campbell-Hausdorff formula for Lie algebras that are finite-dimensional subalgebras of the Lie algebra of vector fields on an appropriate manifold. In particular, we have focused on Lie algebras that can be expressed as algebras of contact Hamiltonian functions endowed with the Jacobi bracket. Our original motivation for this investigation was to improve the error analysis of contact splitting numerical integrators [14,48], which we presented in the last section for the special example of the damped harmonic oscillator. Having explicit modified Hamiltonians allows fine estimates on the numerical errors and a precise a priori comparison of the performance of different numerical algorithms.  (24). Along the columns we vary the parameter γ ∈ 1 2 , 2, 4 .
Even though the flow method in its current form has only limited applicability (due to the fact that one needs to work with algebras of vector fields whose flows are integrable by quadratures), there are several possibilities for future research. We believe that the examples presented in the second part of the paper showcase its promises also in other contexts, for example it allowed us to recover in a remarkably simple way the result from [12] and extend it to a larger class of algebras. An adaptation of the techniques presented here to other Lie algebras is currently under investigation and we expect to present the results in a following publication.
The Baker-Campbell-Hausdorff formula is so pervasive in mathematics and physics [5,26,36,39] that it would not be surprising if more applications will emerge from different fields. For instance, one could use the flow method in the context of the Kepler problem, after rewriting it as a harmonic oscillator [41], or for the numerical integration of timedependent contact Hamiltonians. The latter should appear in a future work in the context of the numerical error analysis in the integration of contact Hamiltonians with timedependent forcing and driving. Moreover, squeezed states in quantum optics are quantum states constructed using the unitary representation of the symplectic group elements acting on the vacuum state of the quantum harmonic oscillator [8,22]. Our results could thus be used to enhance the analysis and the implementation of quantum gates, which can be considered as the product of two unitary symplectic operators [8,22]. They could also be extended to explore quantum gates in the context of polymer quantum mechanics once adapted to the polymer representation of the symplectic group [15,20,21]. Moreover, the analysis and use of symplectic and contact integrators is a very active and prolific field e.g. in mechanics [3,14,17,25,47,48], general relativity [43,44] and plasma physics [34,35], and thus we consider that our results can be helpful in these contexts as well.

A Matrix representations
An alternative (and equivalent) approach to compute the Baker-Campbell-Haussdorf formula by the flow method is to consider the matrix representation of the Lie algebra. In this Appendix we treat the examples in Section 3 from this point of view, giving particular attention to the adjoint representation.
If we start from a Hamiltonian (either symplectic or contact) representation of a finite dimensional algebra g, with generators the Hamiltonian functions {f 1 , · · · , f n }, we can represent a general element of g as the vector Exploiting the corresponding (Poisson or Jacobi) bracket {·, ·}, we can associate to each element a ∈ g the operator ad a := {a, ·} , which is called the adjoint operator corresponding to a. Notice that ad a is a linear operator from g to itself and therefore (for finite-dimensional algebras) it provides a way to associate to every element a ∈ g a matrix (once a basis of g has been fixed). The adjoint operator is always a morphism of Lie algebras, and thus provides a matrix representation of g. This representation is particularly useful when g is centerless, as shown by the following proposition [19].
Proposition A.1. The adjoint representation is faithful, meaning that it is an injective morphism, if and only if the algebra is centerless.
If the algebra has a non-trivial center, then the adjoint representation is not faithful, since the center is mapped to the zero operator. This is the case of the Heisenberg algebra in the symplectic representation, for which the center is the one-dimensional subspace generated by the Hamiltonian function H = 1. On the contrary, the adjoint representation of the same algebra as a subalgebra of contact Hamiltonian functions is faithful.

A.1 Matrix representation of the CHA
There is a representation of the contact Heisenberg algebra [12] realised in the space of the real 4 × 4 matrices, provided by the following identifications: where [a, b] := ab − ba and all the remaining commutators vanish. Another representation is the adjoint representation, which can be found as described in the previous section by exploiting the Jacobi bracket. It leads to the new form of the commutation relations (25) (z, a, b, c), (z,ā,b,c) η = (z 1 , a 1 , b 1 , c 1 ), where [12] z 1 = +cz − cz +āb −ba, a 1 = −āc + ac, b 1 = 0, c 1 = 0.
We can use this to extract the matrix representation of the adjoint since M f · (z,ā,b,c) = (z 1 , a 1 , 0, 0), that is, For the sake of simplicity we present only the computations using the first representation, but the results are equivalent if one uses the representation (27). We do not need to explicitly compute the matrix logarithm of the product of two such exponentials and instead we can work directly with the group elements: since the aforementioned matrix logarithm is again an element of the algebra, the coefficients α, β, γ and ζ in f = αq + βp + γs + ζ1 can be identified by equating term by term, immediately leading to (9).
In this case, the matrix representation of the adjoint is In this case, however, it is much easier to apply the flow method: the computation of the matrix exponential of a 5 × 5 matrix is much more cumbersome. One can already foresee how, with larger algebras, this difference will start getting larger and the flow method could become increasingly advantageous.