Joackim Bernier

EXACT SPLITTING METHODS FOR SEMIGROUPS GENERATED BY INHOMOGENEOUS QUADRATIC DIFFERENTIAL OPERATORS

We introduce some general tools to design exact splitting methods to compute numerically semigroups generated by inhomogeneous quadratic differential operators. More precisely, we factorize these semigroups as products of semigroups that can be approximated efficiently, using, for example, pseudo-spectral methods. We highlight the efficiency of these new methods on the examples of the magnetic linear Schrödinger equations with quadratic potentials, some transport equations and some Fokker-Planck equations.

1. Introduction 1.1. Motivation. Before presenting the general setting of this paper, let us motivate our splitting methods on a quick example. Assume that we aim at solving numerically the following evolution equation on R n

∂ t u(t, x) = -(|x| 2 -∆) u(t, x), t ≥ 0, x ∈ R n u(0, x) = u 0 (x), x ∈ R n ,
where n ≥ 1 and u 0 is a smooth and well localized function on R n . Since the Hermite functions diagonalize the harmonic oscillator, a first natural approach, to reach spectral accuracy, consists in using a Fast Hermite Transform. However, if the solution has to be known in the space variables (for example on a grid) at each time step, such a method is quite costly. A cheaper classical method consists in splitting the harmonic oscillator. For example, here, a natural splitting would be the following Strang splitting [START_REF] Alphonse | Smoothing properties of fractional Ornstein-Uhlenbeck semigroups andnull-controllability[END_REF] e -δt(|x| 2 -∆) ≃ e -δt|x| 2 /2 e δt∆ e -δt|x| 2 /2 .

At each time step, this method requires only a Fast Fourier Transform to solve the heat equation. Unfortunately, even if it reaches spectral accuracy in space, it is only second order in time. Furthermore, note that, here, the only way to get classical splitting methods of higher order is to use complex time steps (see e.g. [START_REF] Blanes | Optimized high-order splitting methods for some classes of parabolic equations[END_REF][START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF]) and that the higher the order of the method is, the larger the number of FFTs per time-step is (and so the more costly the method is). However, allowing the sub-time-steps of the splitting method to be nonlinear functions of δ t , the following splitting formula can be established (see subsection 7.4 of [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF] for a proof) [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF] e -δt(|x| 2 -∆) = e -1 2 tanh(δt)|x| 2 e 1 2 sinh(2δt)∆ e -1 2 tanh(δt)|x| 2 . Note that, contrary to the Strang splitting [START_REF] Alphonse | Smoothing properties of fractional Ornstein-Uhlenbeck semigroups andnull-controllability[END_REF], this factorization is exact : there is no remainder term. Consequently, it is much more accurate than [START_REF] Alphonse | Smoothing properties of fractional Ornstein-Uhlenbeck semigroups andnull-controllability[END_REF] and the time step can be taken quite large (the only possible restrictions coming from the spatial discretization). Furthermore, it is as cheap as [START_REF] Alphonse | Smoothing properties of fractional Ornstein-Uhlenbeck semigroups andnull-controllability[END_REF] to compute.

Roughly speaking, in this paper, we explain how and why the evolution equations associated with a large class of operators called inhomogeneous quadratic differential operators, can be solved by some exact splitting methods. These exact splittings are similar to (2): they provide numerical methods of infinite order (spectral in space and exact in time) using only a small number of FFTs. Sometimes the coefficients of these splittings are not given by an explicit formula depending on δ t (like in (2)), nevertheless, in any case, we provide an iterative method to compute them efficiently.

We warn the reader that the example of the exact splitting (2) may be somehow misleading about our methods. They are not just a generalization of the classical splitting where the sub-time-steps would become nonlinear functions of δ t (as in [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF]). Most of the time, the exact splitting requires the computation of solutions of evolution equations we would not compute in a classical splitting method. For example, to solve the Fokker-Planck equation

∂ t u(t, x, v) + v∂ x u(t, x, v) = ∂ v (v + ∂ v )u(t, x, v), x, v ∈ R, t ≥ 0, u(0, x, v) = u 0 (x, v),
x, v ∈ R, with an exact splitting, we have to solve some Schrödinger equations (see Proposition 4.7 for details).

1.2. General context and result. We consider the problem of the numerical resolution by splitting methods of linear partial differential equations of the form

(3) ∂ t u(t, x) = -p w u(t, x), t ≥ 0, x ∈ R n u(0, x) = u 0 (x), x ∈ R n
where n ≥ 1, u 0 ∈ L 2 (R n ) and p w is an inhomogeneous quadratic differential operator acting on L 2 (R n ).

First, let us precise what we hear by inhomogeneous quadratic differential operator. This terminology is used (for example by Hörmander in [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF]) to denote the Weyl quantization of a polynomial function of degree 2 or less. More precisely, if p is a polynomial function on C 2n of degree 2 or less (called symbol) whose decomposition in coordinates is given by p(X) = t XQX + t Y X + c where X = t (x 1 , . . . , x n , ξ 1 , . . . , ξ n ), Q is a symmetric matrix of size 2n with complex coefficients, Y ∈ C 2n is a vector and c ∈ C is a constant, the Weyl quantization of p is the operator acting on L 2 (R n ) defined by

p w = t x -i∇ Q x -i∇ + t Y x -i∇ + c.
Note that, for such polynomials, this formula coincides with the usual (and much more general) formula defining the Weyl quantization through oscillatory integrals (see e.g., Section 18.5 in [START_REF] Hörmander | The analysis of linear partial differential operators. III, Pseudo-differential operators[END_REF] or Chapter 1 in [START_REF] Nicola | Global Pseudo-Differential Calculus on Euclidean Spaces[END_REF]).

The Weyl quantization allows in many situations the deduction of some properties of the operator from properties of its symbol. For example, as stated in the following proposition, if the real part of p is bounded below on R 2n then (3) is globally well posed.

Proposition 1.1. If p is a polynomial of degree 2 or less on C 2n whose real part is bounded below on R 2n then -p w generates a strongly continuous semigroup on L 2 (R n ) denoted by (e -tp w ) t∈R + .

This proposition is very classical and relies on the Hille-Yosida Theorem. For example, a proof is given by Hörmander in [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF] (pages 425-426) when p is quadratic, but its proof can clearly be extended to the case where p is inhomogeneous (see e.g. Theorem 4.7 in [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF]).

We recognize that the class of equations given by [START_REF] Ameres | Splitting methods for Fourier spectral discretizations of the strongly magnetized Vlasov-Poisson and the Vlasov-Maxwell system[END_REF], where p is a polynomial of degree 2 or less, may seem too elementary to require the use of specific methods to solve them. However, it is usual to have to solve them as sub-steps of splitting methods for more sophisticated equations. Thus it is crucial to have robust methods to compute them very efficiently. Furthermore, for many of these models, in some relevant regimes, their dynamics are the leading part of the dynamics. So, it is crucial to be able to compute them with as much accuracy as possible.

For example, the linear part of some nonlinear Schrödinger equations describing some rotating Bose-Einstein condensates (see e.g. [START_REF] Bader | Fourier-splitting methods for the dynamics of rotating Bose-Einstein condensates[END_REF][START_REF] Bao | Dynamics of rotating Bose-Einstein condensates and its efficient and accurate numerical computation[END_REF]) is of the form [START_REF] Ameres | Splitting methods for Fourier spectral discretizations of the strongly magnetized Vlasov-Poisson and the Vlasov-Maxwell system[END_REF] where

p w = (-i|ξ| 2 -iV (x) -iBx • ξ) w = i∆ -Bx • ∇ -iV (x).
where V : R n → R is a quadratic external potential and B is a real skew symmetric matrix of size n associated with a constant external magnetic field. The formalism of (3) allows also to consider transport equations associated with affine vectors fields. Indeed, if B is a square real matrix of size n and y ∈ R n then

(iBx • ξ + iy • ξ + TrB) w = (Bx + y) • ∇.
Even if such transport equations are essentially trivial, their resolution is required to compute, using splitting methods, the solutions of many kinetic equations (e.g. the Vlasov-Maxwell equations with a constant external magnetic field, see [START_REF] Bernier | Splitting methods for rotations: application to Vlasov equations[END_REF][START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF]). Finally, equations like (3) can also describe some phenomena of diffusions. For example, the generalized Ornstein-Uhlenbeck operators are of the form

(ξ • Aξ + iBx • ξ + x • Rx + TrB) w = -∇ • A∇ + Bx • ∇ + x • Rx
where A, R are some real nonnegative symmetric matrices of size n and B is a real matrix n. Note that these operators include Fokker-Planck and Kramers-Fokker-Planck operators (see e.g. [START_REF] Dujardin | Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations[END_REF][START_REF] Hérau | Tunnel effect for the Kramers-Fokker-Planck type operators[END_REF]). Some of the equations of the form (3) are easy to solve numerically using pseudo-spectral methods. For example, to solve the heat equation or to compute a shear, it is enough to do some Fast Fourier Transforms. Similarly, in the spirit of the splitting methods, it is not very costly to solve successively some of these equations. More precisely, let us define, in this context, the operators we consider as easily computable using standard pseudo-spectral methods.

Definition 1.2. An operator acting on L 2 (R n ) can be computed by an exact splitting if it can be factorized as a product of operators of the form e α∂x j , e iαx j , e ia(∇) , e ia(x) , e αx k ∂x j , e -b(x) , e b(∇) , e γ where α ∈ R, γ ∈ C, a, b : R n → R are some real quadratic forms, b is nonnegative, j, k ∈ 1, n and k = j. As usual, a(∇) (resp. b(∇)) denotes the Fourier multiplier associated with -a(ξ) (resp. -b(ξ)), i.e. a(∇) = (-a(ξ)) w .

Note that if an operator can be computed by an exact splitting then it is bounded. The following theorem justifies why we focus on splitting methods for semigroups generated by inhomogeneous quadratic differential operators.

Theorem 1.3. If p is a polynomial of degree 2 or less on C 2n whose real part is bounded below on R 2n then e -p w can be computed by an exact splitting.

Remark 1.4. Actually, we prove a slightly stronger result : the quadratic forms associated with a and b in the definition of Definition 1.2 can be chosen diagonal. Nevertheless, it is not necessarily relevant from a numerical point of view because the sub-steps required to diagonalize the quadratic forms can be costly. This theorem is proved in the subsection 5.2 of the Appendix. Unfortunately, its proof does not provide an efficient way to split semigroups (minimizing, for example, the number of sub-steps or the number Fast Fourier Transforms required to approximate the exponentials). Nevertheless, as illustrated in Section 4, on many examples, paying attention to the particular structure of each semigroup, we are able to design optimized exact splittings.

This work can be considered as the theoretical part of a more general study. Indeed, a second work [START_REF] Bernier | Exact splitting methods for kinetic and Schrödinger equations[END_REF], written in collaboration with Nicolas Crouseilles and Yingzhe Li, deals with the implementation of these methods and compares numerically their efficiency, their accuracy and their qualitative properties with respect to the existing methods. We also couple these exact splitting methods with classical methods in order to solve some nonlinear equations and some non quadratic linear equations.

Outline of the work . In Section 2, we develop the notion of exact splitting in a more general framework and specify its links with the classical splittings. It will naturally lead to a general theorem to design efficient exact splittings for many linear ordinary differential equations. In Section 3, we explain how the theory of the Fourier Integral Operators developed by Hörmander in [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF] can be used to transform exact splittings of linear ordinary differential equations into exact splittings of semigroups. And, finally, in Section 4, we apply the results of the previous sections to obtain some efficient exact splittings for the magnetic linear Schrödinger equations with quadratic potentials, some transport equations and some Fokker-Planck equations.

Acknowledgment . The author thanks P. Alphonse, N. Crouseilles and Y. Li for many enthusiastic discussions about this work. Research of the author was supported by ANR project NABUCO, ANR-17-CE40-0025.

Notations and conventions. Let us define some classical notations used in this paper.

• To get convenient notations, most of the time, we denote by t (instead of δ t ) the time-step of our methods. • I n denotes the identity matrix on R n and J 2n denotes the matrix of the canonical symplectic form of R 2n , i.e.

J 2n := I n -I n .
• By convention, the empty spaces in the matrix notations refer to coefficients equal to zero (see for example the definition of J 2n above).

• If A is a matrix t A denotes the transpose of A. • If K ∈ {R, C}, M n (K)
denotes the algebra of the square matrices of size n and S n (K) = {M ∈ M n (K), t M = M } denotes the vector space of the symmetric matrices. • For K ∈ {R, C}, we will use the following classical groups of matrices

GL n (K) = {M ∈ M n (K), M is invertible} SL n (K) = {M ∈ M n (K), det M = 1} Sp 2n (K) = {M ∈ M 2n (K), t M J 2n M = J 2n } O n (R) = {M ∈ M n (R), t M M = I n }.
and their associated Lie algebras

gl n (K) = M n (K) sl n (K) = {M ∈ M n (K), TrM = 0} sp 2n (K) = {J 2n Q, Q ∈ S 2n (K)} so n (R) = {M ∈ M n (R), t M = -M }.
where Tr denotes the trace and the Lie bracket is formally defined through the relation [A, B] := AB -BA.

• A real Lie algebra of matrices is a sub-Lie-algebra of gl n (R) for some n ≥ 0.

• We equip the space of the polynomials of degree 2 or less on C 2n of its structure of Lie algebra induced by the canonical Poisson bracket. This one being defined for two polynomials

p 1 , p 2 on C 2n ≡ C n x × C n ξ by {p 1 , p 2 } := n j=1 ∂ ξ j p 1 ∂ x j p 2 -∂ x j p 1 ∂ ξ j p 2 .
Note that if K ∈ {R, C}, the space of the quadratic forms on K 2n is a Lie algebra naturally isomorphic to sp 2n (K). • If g is a Lie algebra, ad : g → g denotes its adjoint representation, i.e. ∀x, y ∈ g, ad x y := [x, y].

• We consider the natural action of the analytic functions on M n (C) defined by the holomorphic functional calculus (see VII-3 of [START_REF] Dunford | Linear operators. Part I[END_REF] for details). By abuses of notations, if f is analytic on a domain Ω ⊂ C and M ∈ M n (C) has its spectrum included in Ω then (f (z))(M ) is just an other way to denote f (M ). Note that, if f : D(z c , ρ) → C is an analytic function on the complex disk of center z c and radius ρ > 0 and M ∈ M n (C) has its spectrum included in D(z c , ρ) then f (M ) can be defined by the convergent series

f (M ) = k∈N f (k) (z c ) k! (M -z c I n ) k .
• S(R n ) denotes the Schwartz space on R n .

• In a non-commutative setting, the notation to denote a product can be quite ambiguous. In this paper, we adopt the following natural convention. If I is a totally ordered finite set and (g j ) j∈I ∈ M I , where M is a monoid1 , then j∈I g j := g ι 1 . . . g ι ♯I where ι is the increasing bijection from 1, ♯I to I.

Exact splitting methods

The original setting of the splitting methods (see e.g [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]) consists in considering linear equations2 of the form (4)

∂ t u = Lu = L (1) u + • • • + L (k) u u(t = 0) = u 0 .
whose solution is denoted u(t) = exp(tL)u 0 and such that exp(tL (1) ), . . . , exp(tL (k) ) are nicer than exp(tL). In this context, nicer usually means cheaper to compute. Then the approximations of u at times t n = nδ t are got compositing n times an approximation Ψ δt of exp(δ t L) where Ψ δt is a composition of operators of the form exp(α σ j δ t L (σ j ) ) with 1 ≤ σ j ≤ k and α j ∈ R. The most classical methods are the Lie splitting where

Ψ (Lie) δt = exp(δ t L (1) ) . . . exp(δ t L (k) )
and the Strang splitting 1) ).

Ψ (Strang) δt = exp( δ t 2 L (1) ) . . . exp( δ t 2 L (k-1) ) exp(δ t L (k) ) exp( δ t 2 L (k-1) ) . . . exp( δ t 2 L (
These methods are respectively of order 1 and 2. It means that they provide respectively approximations of order δ t and δ 2 t of u(t n ). Similar methods can be derived to obtain splitting methods of arbitrarily high order. However, note that, the higher the order is, the higher the number of step is and so the higher the cost of the method is. Furthermore, the only way to get methods of order higher than 2 where the α σ j are nonnegative is to allow the α σ j to be complex. In this case it is possible to design schemes of high order where the real parts of the α σ j are nonnegative (see [START_REF] Blanes | Optimized high-order splitting methods for some classes of parabolic equations[END_REF][START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF]). Note that this is crucial when irreversible equations are considered.

The setting of the exact splitting is the same as the setting of the classical ones. Nevertheless, we have to assume that exp(tL (1) ), . . . , exp(tL (k) ) are nicer than exp(tL) because L (1) , . . . , L (k) belong to some vector spaces E 1 , . . . , E k . This assumption is natural because usually the exp(tL (j) ) are nice due to a particular structure of each L (j) (e.g. nilpotent, diagonal...). Consequently, the idea of the exact splitting consists in looking for Ψ δt as a product of operators of the form exp(δ t L (σ j ) j,δt ) with 1 ≤ σ j ≤ k, L (σ j ) j,δt ∈ E σ j and to ask that Ψ δt is exactly equal to exp(δ t L). In other words, we are looking for a factorization of exp(δ t L) as a product of operators of the form exp(δ t L

(σ j ) j,δt ) with L (σ j ) j,δt ∈ E σ j .
Remark 2.1. To avoid any possible confusion note that exp(δ t L (σ j ) j,δt ) does not denote the solution of a non-autonomous equation but the solution of ∂ t u = L (σ j ) j,δt u at time t = δ t . Actually, the existence of an exact splitting can be seen as an inverse problem with respect to the classical backward error analysis of splitting methods. In the context of the splitting method, this analysis is realized through the Baker-Campbell-Hausdorff formula (see [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]). It states that [START_REF] Bader | Fourier-splitting methods for the dynamics of rotating Bose-Einstein condensates[END_REF] exp(α

σ 1 δ t L (σ 1 ) ) . . . exp(α σ ℓ δ t L (σ ℓ ) ) = exp(δ t Ω δt,L,α,σ )
where Ω δt,L,α,σ can be expanded in powers of δ t and the operator associated with the power δ n t is obtained as a linear combination of n Lie brackets of the L (j) . For example, we have

Ω δt,L,α,σ = ℓ j=1 α σ j L (σ j ) + δ t 2 i<j α σ j α σ i [L (σ i ) , L (σ j ) ] + O(δ 2 t )
where [L (σ i ) , L (σ j ) ] = L (σ i ) L (σ j ) -L (σ j ) L (σ i ) . Note that, in general, the formula (5) and the expansion of Ω δt,L,α,σ have to be understood in the sense of the formal series in δ t . Nevertheless, if L (1) , . . . L (k) belong to a real Lie algebra of matrices then these expansions converge when δ t is small enough (see e.g. [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]). Now if we consider a product of operators of the form exp(δ t L (σ j ) j,δt ) with 1 ≤ σ j ≤ k, L (σ j ) j,δt ∈ E σ j , where the E j are some vector subspaces of a same real Lie algebra of matrices then the Baker-Campbell-Hausdorff formula states that it is of the form [START_REF] Bader | Fourier methods for the perturbed harmonic oscillator in linear and nonlinear Schrödinger equations[END_REF] exp(α

σ 1 δ t L (σ 1 ) 1,δt ) . . . exp(α σ ℓ δ t L (σ ℓ ) ℓ,δt ) = exp(δ t Ω δt,L,σ ),
where Ω δt,L,σ admits an expansion similar to the expansion of Ω δt,L,α,σ . Consequently, to get an exact splitting method we just have to design L (σ ℓ ) j,δt in order to cancel all the Poisson brackets in the expansion of Ω δt,L,σ , i.e. we have to solve [START_REF] Bambusi | Reducibility of the quantum harmonic oscillator in d-dimensions with polynomial time-dependent perturbation[END_REF] Ω δt,L,σ = k j=1 L (j) .

Since it can be shown that Ω δt,L,σ is smooth with respect to δ t and (L (σ j ) j,δt ) 1≤j≤ℓ , it is natural to try to solve [START_REF] Bambusi | Reducibility of the quantum harmonic oscillator in d-dimensions with polynomial time-dependent perturbation[END_REF] with the Implicit Function Theorem.

In this way, we prove the following theorem for which we have many applications. From now, to get convenient notations, we denote t instead of δ t .

⋆ = b ⋆,1 + • • • + b ⋆,m ∈ b 1 ⊕ • • • ⊕ b m = b is such that (8) (b 1 ⊕ • • • ⊕ b m ) + ad b⋆ (s) is a real Lie algebra
then there exist t 0 > 0 and an analytic function t ∈ (-t

0 , t 0 ) → (b t , s t ) ∈ b × s such that b 0 = b ⋆ and (9)
∀t ∈ (-t 0 , t 0 ), e tb⋆ = e -tst e tb t,1 . . . e tbt,m e tst .

Remark 2.3. Before proving this theorem, let us do some remarks.

• In most of the applications, s can be chosen as a one of the b j . Consequently, assuming that s = b 1 , the notations of this theorem are consistent with the notations previously introduced in this section through the identifications

δ t = t, k = m, E j = b j , b ⋆ = L (1) + • • • + L (m) , ℓ = k + 2, L (1)
1,δt = -s t , L

(j-1) j,δt = b t,j-1 , L (1)
ℓ,δt = s t where j ∈ 2, m + 1 . • Since the proof of this theorem relies on the Implicit Function Theorem, the coefficients of the exact splitting (i.e. s t and b t ) can be efficiently computed by an iterative method.

Unfortunately, this method require some notations introduced in the proof of Theorem 2.2 to be presented. Consequently, it is introduced just after the proof. • In practice, as we will see in Section 4, it may be useful to use Theorem 2.2 just to determine a priori the form of an exact splitting and then to get analytically the associated coefficients (using for example a formal computation software). • The assumption (8) of Theorem 2.2 is a bit too strong. The optimal assumption seems that Ω δt,L,σ (defined implicitly by ( 6)) belongs to the vector space defined in (8) for all L and δ t . Paying attention to the expansion of Ω δt,L,σ given by the Baker-Campbell-Hausdorff formula, it essentially means that we do not need to ask this space to contain Lie brackets of elements belonging to a same space b j . Note that, such a generalization would be necessary to justify a priori the form of the exact splitting of the Example 2.3 of [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF] for the Kramer-Fokker-Planck operator.

Proof of Theorem 2.2. Naturally, following the assumption of the theorem, we introduce the real Lie algebra of matrices, denoted g and defined by

(10) g = (b 1 ⊕ • • • ⊕ b m ) + ad b⋆ (s).
Applying the Baker Campbell Hausdorff formula, we get a neighborhood of the origin in b × s, denoted V , and an analytic function F : V → g such that for all real t and all (b, s) ∈ b × s such that t(b, s) ∈ V , we have e -ts e tb 1 . . . e tbm e ts = exp

  tb + t 2 [b, s] + t 2 2 1≤i<j≤m [b i , b j ] + t 3 F (tb, ts)   .
Consequently, we aim at solving the equation

(11) b ⋆ = b + tad b (s) + t 2 1≤i<j≤m [b i , b j ] + t 2 F (tb, ts).
To solve this equation we have to introduce some notations. First, from the decomposition (10) of g, we deduce naturally that there exists a complementary space to b in g denoted r, a subspace of s denoted s ′ and a vector space isomorphism Ψ : s ′ → r such that ( 12)

Ψ = Π r • ad b⋆
where Π b , Π r are the canonical projections associated with the decomposition

(13) g = b ⊕ r.
Then, let s ⋆ ∈ s ′ be defined by ( 14)

s ⋆ = - 1 2 (Ψ -1 • Π r )   1≤i<j≤m [b ⋆,i , b ⋆,j ]   , let V ⋆ be a neighborhood of (b ⋆ , s ⋆ ) in b × s ′ and let t ⋆ > 0 be such that (-t ⋆ , t ⋆ )V ⋆ ⊂ V.
So, to solve [START_REF] Blanes | Optimized high-order splitting methods for some classes of parabolic equations[END_REF], we are going to apply the Implicit Function Theorem in (0, b ⋆ , s ⋆ ) to the function

G : (-t ⋆ , t ⋆ ) × V ⋆ → b × s ′ (t, b, s) → (G b , G s ′ )(t, b, s) where G b (t, b, s) = b -b ⋆ + t Π b R(t, b, s). with R(t, b, s) = [b, s] + 1 2 1≤i<j≤m [b i , b j ] + tF (tb, ts),
and

ΨG s ′ (t, s (m) , b, s (r) ) = Π r • ad b s + 1 2 Π r 1≤i<j≤m [b i , b j ] + t Π r F (tb, ts).
First, observe that, by construction of G, if G(t, b, s) = 0 then b, s is a solution of [START_REF] Blanes | Optimized high-order splitting methods for some classes of parabolic equations[END_REF]. Indeed, by construction, the equation ( 11) can be rewritten as

(G b + tΨG s ′ )(t, b, s) = 0.
Then observe that, by construction of s ⋆ we have

G(0, b ⋆ , s ⋆ ) = 0.
Consequently, since G is clearly an analytic function, to conclude the proof applying the Implicit Function Theorem, we just have to prove that the partial differential of G with respect to (b, s) in (0, b ⋆ , s ⋆ ) is invertible. Indeed, using a natural matrix representation, this differential is

I b L b,⋆ I s ′ where L b,⋆ : b → s ′ is
an explicit linear map depending on b ⋆ and s ′ . This matrix being triangular, it is clearly invertible and its invert is

I b -L b,⋆ I s ′
Now, let us present an iterative method to determine the coefficients b t and s t given by the Theorem 2.2.

Proposition 2.4. Assume that the assumption of Theorem 2.2 is satisfied, let Ψ be defined by [START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF], s ⋆ be defined by [START_REF] Chen | 3D volume rotation using shear transformation[END_REF] and Π b , Π r be the projections canonically associated with the decomposition [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF]. There exists t 1 ∈ (0, t 0 ), such that for all t ∈ (-t 1 , t 1 ) the sequence (b

(k) t , s (k) t ) k∈N if

well defined by induction as follow :

Initially, we have b

(0) t = b ⋆ and s (0) t = s ⋆ and for k ≥ 0 (15) b (k+1) t = b (k) t + b ⋆ - Π b g (k) s (k+1) t = s (k) t -t -1 Ψ -1 Π r g (k)
where

g (k) = t -1 log(e -ts (k) t e tb (k) t,1 . . . e tb (k)
t,m e ts (k) t ). Furthermore, there exist C > 0 independent of t and k such that for all k ≥ 0 and all t ∈ (-t 1 , t 1 ), we have

|b (k) t -b t | + |s (k) t -s t | ≤ C 2 -k .
Note that, here, as usual, the log function denotes the principal determination of the matrix logarithm. We have chosen to present an iterative method as elementary as possible. However, computing the Lie brackets associated with the operators L s,⋆ and L b,⋆ it would be possible to get a natural iterative method whose convergence rate would be τ k instead of 2 -k with τ a linear function of |t|.

Proof of Proposition 2.4. It follows from the proof of the theorem 2.2 that the sequence is defined by a relation of the kind

(b (k+1) t , s (k+1) t ) = W t (b (k) t , s (k) t )
where for all t ∈ (-t 0 , t 0 ), W t is a smooth function on a same ball B of center (b ⋆ , s ⋆ ) and t → W t is smooth. Consequently, to prove that the sequence is well defined if |t| is small enough and that W t has a fix point we just have to prove that W t is a contraction mapping for a well chosen norm. Indeed, it follows of the proof of the theorem 2.2 that the differential of

W 0 in (b ⋆ , s ⋆ ), denoted dW 0 (b ⋆ , s ⋆ ), is 0 b -L b,⋆ 0 s ′ .
Thus, since it is nilpotent, applying for example the Lemma 5.6.10 of [START_REF] Horn | Matrix Analysis[END_REF], we get a matrix norm

• ⋆ such that dW 0 (b ⋆ , s ⋆ ) ⋆ ≤ 1/4. Now, since W is smooth on (-t 0 , t 0 ) × B and W 0 (b ⋆ , s ⋆ ) = (b ⋆ , s ⋆ ),
t ) converges to this fix point with the rate 2 -k .

Finally, we conclude this proof observing that, by construction and continuity, there exists t 1 ∈ (0, t 2 ) such that for all t ∈ (-t 1 , t 1 ), (b t , s t ) belongs to B ⋆ and is a fix point of W t .

Exact classical-quantum correspondance

As we have seen in Section 2, exact splittings can be designed (using for example Theorem 2.2) to solve linear ordinary differential equations. We aim at designing exact splittings to solve some partial differential equations. So, first, it is natural to focus on linear transport equations. Indeed, the formula of the characteristics ( 16)

e t Bx•∇ u = u • e tB
where u ∈ L 2 (R n ) and B is square matrix of size n, provides a natural way to transform an exact splitting at the level of the linear ordinary differential equations into an exact splitting at the level of the linear transport equations.

For example, we have the following exact splitting for the two-dimensional rotations

(17) exp 0 tan(θ/2) 0 0 exp 0 0 -sin θ 0 exp 0 tan(θ/2) 0 0 = exp(θJ 2 ),
where θ ∈ (-π, π). At the level of the associated transport equation, this formula can be written e t(x 2 ∂x 1 -x 1 ∂x 2 ) = e tan(t/2)x 2 ∂x 1 e -sin(t)x 1 ∂x 2 e tan(t/2)x 2 ∂x 1 . This factorization is very useful to compute rotations since, using semi-Lagrangian methods, it only requires one dimensional interpolations instead of a two dimensional interpolations as we could expect. The formula is well known in image processing and has been used for decades to rotate images (see e.g. [START_REF] Paeth | A fast Algorithm for General Raster Rotation[END_REF]). Two recent papers deal with the applications of this decomposition for the numerical resolution of kinetic equations (see [START_REF] Ameres | Splitting methods for Fourier spectral discretizations of the strongly magnetized Vlasov-Poisson and the Vlasov-Maxwell system[END_REF][START_REF] Bernier | Splitting methods for rotations: application to Vlasov equations[END_REF]). This factorization has also been extended to compute 3 dimensional rotations with only one dimensional interpolations (see e.g. [START_REF] Chen | 3D volume rotation using shear transformation[END_REF][START_REF] Welling | Rotation of 3D volumes by Fourier-interpolated shears[END_REF]). Note that, in Subsection 4.1, we extend this kind of decomposition in any dimension and for more general transforms (including rotations and dilatations).

The inhomogeneous quadratic differential operators and the semigroups they generate enjoy some strong and specific properties providing a more general way to transform an exact splitting at the level of the linear ordinary differential equations into an exact splitting at the semigroup level.

Indeed, usually, the Lie bracket of two pseudo-differential operators is also a pseudodifferential operator. Its symbol is given by the Moyal bracket of the symbols. In many applications, the Moyal bracket admits a natural expansion whose leading part3 is given by the Poisson bracket of the symbols (see e.g. [START_REF] Hörmander | The analysis of linear partial differential operators. III, Pseudo-differential operators[END_REF]). Furthermore, in the particular case of the inhomogeneous quadratic operators all the higher order terms vanish and we have [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] [q w 1 , q w 2 ] = -i{q 1 , q 2 } w where q 1 , q 2 are some polynomials of degree 2 or less on C 2n . The formula [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] is especially useful to realize some changes of unknowns in order to put some operators in normal form (see e.g. [START_REF] Bambusi | Reducibility of the quantum harmonic oscillator in d-dimensions with polynomial time-dependent perturbation[END_REF] for an application to the reducibility that can also be seen as a factorization of exponentials). Note that, applying formally the Baker-Campbell-Hausdorff formula at the level of the operator, [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] suggests a way to transform an exact splitting at the level of the linear ordinary differential equations into an exact splitting at the level of the semigroups generated by inhomogeneous quadratic forms. We refer the reader to the beginning of the Section 3 of [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF] for details about this heuristic. It is used in [START_REF] Chin | Fourth-order algorithms for solving the imaginary-time Gross-Pitaevskii equation in a rotating anisotropic trap[END_REF][START_REF] Bader | Fourier methods for the perturbed harmonic oscillator in linear and nonlinear Schrödinger equations[END_REF][START_REF] Bader | Fourier-splitting methods for the dynamics of rotating Bose-Einstein condensates[END_REF] to design some exact splitting methods to solve linear Schrödinger equations with harmonic potentials and rotating terms. The result suggested by these formal computations is made rigorous in Proposition 3.6 below, and relies on the representation of the semigroups as Fourier Integral Operators (detailed in Theorem 3.5 below).

To present this notion, we need to introduce some basic notations and associated properties.

Definition 3.1. T is a non-negative complex symplectic linear bijection on C 2n , and we denote

T ∈ Sp + 2n (C), if T ∈ Sp 2n (C) and ∀X ∈ C 2n , t X t T (-iJ 2n )T X - t X(-iJ 2n )X ∈ R + .
Note that this set is naturally equipped with a structure of monoid (i.e. it is stable by multiplication and the identity belongs to Sp + 2n (C)). Definition 3.2. If q is a quadratic form on C 2n then its Hamiltonian flow at time t ∈ R, denoted by Φ q t , is the flow at time t of the linear ordinary differential equation

(19) d dt z = -iJ 2n ∇q(z)
where ∇q = t (∂ x 1 q, . . . , ∂ xn q, ∂ ξ 1 q, . . . , ∂ ξn q).

Note that, in particular, the linear ordinary differential equation ( 19) can be solved using an exponential and thus we have [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF] Φ q t = e -2itJ 2n Q where Q is the matrix of q. The following proposition summarizes some elementary properties of these Hamiltonian flows that will be used all along this paper. Proposition 3.3. Let q 1 , q 2 , q be some quadratic forms on C 2n and T ∈ Sp 2n (C) then the following properties holds

i) ∀t ∈ R, Φ q t ∈ Sp 2n (C), ii) {q 1 , q 2 } = 0 ⇐⇒ ∀t ∈ R, Φ q 1 t Φ q 2 t = Φ q 1 +q 2 t , iii) ∀t ∈ R, T -1 Φ q t T = Φ q•T t .
Usually, in this context, the second property is called Noether's theorem because it is also equivalent to have q 1 • Φ q 2 t = q 1 for all t ∈ R. Proofs of these properties can be found, for example, in a nonlinear context, in [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF].

The following classical lemma (whose proof is recalled in the subsection 5.3 of the Appendix) links naturally the two previous definitions. Lemma 3.4. If q is a quadratic form on C 2n such that ℜq is nonnegative on R 2n then ∀t ≥ 0, Φ q t ∈ Sp + 2n (C). The following theorem, proved by Hörmander in [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF] (Thm 5.12 and Prop 5.9), is the main tool we use to realize exact splittings. Theorem 3.5 (Hörmander [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF]). There exists a map

K : Sp + 2n (C) → B/U 2 where U 2 = {+1, -1} and B is the unit ball of L (L 2 (R n ))
, the algebra of the bounded operators on L 2 (R n ), which is monoid morphism, i.e.

∀S, T ∈ Sp

+ 2n (C), K (ST ) = K (S)K (T ) and K (I 2n ) = ±id L 2 .
Furthermore, if q is a quadratic form on C 2n such that ℜq is nonnegative on R 2n then K (Φ q t ) = ±e -tq w . In [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF], K is defined through an explicit but heavy formula that is not relevant for us here. An operator of the form K (T ) with T ∈ Sp + 2n (C) is called a Fourier Integral Operators. It provides a natural way to transform an exact splitting at the level of the linear ordinary differential equations into an exact splitting at the level of the semigroups generated by quadratic differential operators (up to an argument of continuity to remove the uncertainty of sign).

To the best of our knowledge, the idea to use the Fourier integral operators to get an exact splitting has been introduced by Paul Alphonse and the author in [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF]. We aimed at characterizing the regularizing effects of the semigroups generated by non-selfadjoint quadratic differential operator. Our splitting provides a quite explicit representation of the polar decomposition of this semigroup. Note that it is the decomposition we would get applying Theorem 2.2 and Theorem 3.5 with s = {0}, b 1 = i sp 2n (R) and b 2 = sp 2n (R).

In order to give a corollary of Theorem 3.5 well suited to get exact splittings for semigroups generated by inhomogeneous quadratic differential operators, we have to introduce a last elementary notation. If p : C 2n = C n x × C n ξ → C is a polynomial of degree 2 or less that we write in coordinates as

p = t XQX + t Y X + c where X = t (x 1 , . . . , x n , ξ 1 , . . . , ξ n ), Q ∈ S 2n (C), Y ∈ C 2n and c ∈ C, then Pp : C 2n+2 = C n+1 x × C n+1 ξ → C is the complex quadratic form defined by Pp = t XQX + t Y Xx n+1 + cx 2 n+1 .
Proposition 3.6. Let p 1,t , . . . , p m+1,t : C 2n → C be some polynomials of degree 2 or less depending continuously on t ∈ [0, t 0 ] for some t 0 > 0, whose real part is uniformly bounded below on R 2n and satisfying

(21) ∀t ∈ [0, t 0 ], Φ Pp 1,t t . . . Φ Ppm,t t = Φ Pp m+1,t t
then we have [START_REF] Horn | Matrix Analysis[END_REF] ∀t ∈ [0, t 0 ], e -tp w 1,t . . . e -tp w m,t = e -tp w m+1,t .

The proof is given at the end of this section. Let us mention that this proposition is also a corollary of the Theorem 2.3 of [START_REF] Viola | The elliptic evolution of non-self-adjoint degree-2 Hamiltonians[END_REF]. Note that if the polynomials are homogeneous (i.e. if they are some quadratic forms) then the P can be removed. Furthermore, as a corollary of the proof, if we do not assume that the polynomials are uniformly bounded below and that the polynomials depend continuously on t, then [START_REF] Horn | Matrix Analysis[END_REF] holds up to an uncertainty of sign.

It may also be interesting to have exact splittings to compute evolution operators generated by non-autonomous inhomogeneous quadratic differential operators (like it is done for the non-autonomous linear magnetic Schrödinger equations in [START_REF] Bader | Fourier-splitting methods for the dynamics of rotating Bose-Einstein condensates[END_REF]). Here, we chose for conciseness to do not consider the non-autonomous case. Nevertheless, let us mention that it would be possible to generalize Proposition 3.6 to deal with non-autonomous equations using the generalization of the Theorem 3.5 proven in [START_REF] Pravda-Starov | Generalized Mehler formula for time-dependent non-selfadjoint quadratic operators and propagation of singularities[END_REF]. With such a generalization, the exponential of matrices become naturally the solutions of non-autonomous linear ordinary differential equations (which can be studied similarly using Magnus expansions, see [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]).

In order to illustrate Proposition 3.6, let us give an elementary application (Section 4 being devoted to more sophisticated applications). Indeed, the formula (17) used to compute rotations can clearly be extended analytically for θ ∈ -2it where t ∈ R. Since, up to a transposition, this formula can be written as

Φ (tanh t)/(2t) x 2 t Φ (sinh 2t)/(2t) ξ 2 t Φ (tanh t)/(2t) x 2 t = Φ |x| 2 +|ξ| 2 t
we deduce of the Proposition 3.6 that the exact splitting (2) for the harmonic oscillator (presented in the introduction) holds.

Before focusing on the proof of Proposition 3.6, let us present the following Lemma that is useful in practice to establish the factorization [START_REF] Hörmander | The analysis of linear partial differential operators. III, Pseudo-differential operators[END_REF]. It expresses the triangular nature of the equation [START_REF] Hörmander | The analysis of linear partial differential operators. III, Pseudo-differential operators[END_REF]. Its proof is postponed in subsection 5.1 of the appendix. Lemma 3.7. Let p 1 , . . . , p m be some polynomials of degree 2 or less on C 2n whose decompositions are p j = q j + ℓ j + c j where q j quadratic forms, ℓ j are linear forms and c j are complex numbers. Then we have

(23) Φ Pp 1 1 . . . Φ Ppm 1 = Φ Pp m+1 1
if and only if

(24)                Φ q 1 1 . . . Φ qm 1 = Φ q m+1 1 m j=1 ℓ j • Υ q j k>j Φ q k 1 = ℓ m+1 • Υ q m+1 m j=1 c j + κ j + σ j = c m+1 + κ m+1
where, denoting Q (j) the matrix of q j and L (j) ∈ C 2n the matrix of ℓ j

Υ q j = e z -1 z (-2iJ 2n Q (j) ), κ j = k∈N 4 k (2k + 3)! q j ((J 2n Q (j) ) k J 2n t L (j) ), σ j = - i 2 j-1 p=1 L (p) Υ qp   j-1 k=p+1 Φ q k 1   Υ q j J 2n t L (j) .
In this lemma, it is relevant to note that if q 1 , . . . , q m+1 satisfy the first equation of ( 24) then the second equation is just a linear system with respect to ℓ 1 , . . . , ℓ m+1 and that if the two first equations are satisfied then the last equation is linear with respect to c 1 , . . . , c m+1 .

In order to prove Proposition 3.6, we introduce some technical lemmas that will be crucial.

Lemma 3.8. If p : R 2n → R is a real polynomial of degree two or less being bounded below then there exists a nonnegative real quadratic form q : R 2n → R + , Y ∈ R 2n and c ∈ R such that p can be written as

(25) p(X) = q(X -Y ) + c
where X = (x 1 , . . . , x n , ξ 1 , . . . , ξ n ).

Proof of Lemma 3.8. Naturally p can be written in coordinates as

p(X) = t XQX + t ZX + p(0)
where Q ∈ S 2n (R) is the matrix of q and Z ∈ R 2n . Now realizing a canonical factorization, we observe that it is equivalent to prove that p can be written as [START_REF] Pravda-Starov | Generalized Mehler formula for time-dependent non-selfadjoint quadratic operators and propagation of singularities[END_REF] and that Z ∈ Im Q.

In other words, we just have to prove that

∀X 0 ∈ Im Q ⊥ , t X 0 Z = 0. But since if X 0 ∈ Im Q ⊥ we have p(λX 0 ) = λ t X 0 Z + p(0)
and we know, by assumption, that λ → p(λX 0 ) is bounded below then we deduce that t X 0 Z = 0, which conclude this proof.

Corollary 3.9. If p : R 2n → R is a real polynomial of degree two or less being bounded below then we have

P(p -inf R 2n p) ≥ 0
Proof of Corollary 3.9. Applying Lemma 3.8, p can be written as

p(X) = q(X -Y ) + c.
So, by definition, we have

P(p -c) = q(X -x n+1 Y ).
Thus, since q is nonnegative, P(pc) is also nonnegative and c is the infimum of p.

Lemma 3.10. If ψ ∈ S(R n+1 ) and p is a polynomial of degree 2 or less on C 2n whose real part is nonnegative on R 2n then ℜPp is nonnegative on R 2n+2 , e -(Pp) w ψ ∈ S(R n+1 ) and we have (e -(Pp) w ψ)

|R n ×{1} = e -p w ψ |R n ×{1}
Proof of Lemma 3.10. First, observe that Corollary 3.9 proves directly that ℜPp is nonnegative. Consequently, the semigroup generated by -(Pp) w is well defined (see e.g. Proposition 1.1).

Applying the Theorem 4.2 of [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF], we know that

∀t ∈ [0, 1], φ(t) = e -t(Pp) w ψ ∈ S(R n+1 ) and that φ ∈ C ∞ (R + t × R n+1 ) and satisfies ∀t ∈ R + , ∂ t φ = -(Pp) w φ.
Now decomposing p by homogeneity as

p = q + ℓ + c
where q : R 2n → C is a complex valued quadratic form, ℓ : R 2n → C is a complex valued linear form and c ∈ C, we have

∀t ∈ R + , ∂ t φ = -(q w + x n+1 ℓ w + cx 2 n+1
)φ. Since here we deal with smooth functions, this relation can be evaluated in

x n+1 = 1. Thus we get ∀t ∈ R + , ∂ t ρ = -(q w + ℓ w + c)ρ = -p w ρ where ρ = φ |R + t ×R n ×{1} .
Consequently by definition of the semigroup exp(-tp w ), we have

ρ = exp(-tp w )ρ(0, •) = exp(-tp w )ψ |R n ×{1} .
Finally, evaluating this last relation for t = 1, we get the desired relation.

Using these two lemmas, we give the following proof of the Proposition 3.6.

Proof of Proposition 3.6. In order to obtain a factorization of the semigroup using Theorem 3.5, we have to deal with quadratic forms having a nonnegative real part. So, since the polynomials are uniformly bounded below, we apply Corollary 3.9 to get a constant c > 0 such that

(26) ∀j ∈ 1, m + 1 , ∀t ∈ [0, t 0 ], ℜ P(p j,t + c) ≥ 0.
Since the real part of P(p j,t + c) = Pp j,t + cx Now, we can apply the Theorem 3.5 to deduce that for all t ∈ [0, t 0 ], there exists ε t ∈ {+1, -1} such that ∀t ∈ [0, t 0 ], e -t(Pp 1,t ) w -tcx 2 n+1 . . . e -t(Ppm,t) w -tcx 2 n+1 = ε t e -t(Pp m+1,t ) w -tcmx 2 n+1 .

To get the desired factorization, we have to check that ε t = 1 and to prove that this relation can be evaluated at x n+1 = 1.

First, we focus on the sign and we observe that if t = 0 all the exponentials are equal to the identity so we have ε 0 = 1. Thus, since [0, t 0 ] is connected, we just have to prove that t → ε t is continuous to deduce that ε t = 1 for all t ∈ [0, t 0 ]. Let φ ∈ S(R n ) \ {0} be a Schwarz function on R n non identically equals to zero (for example a gaussian). Since e -t(Ppt) w -tcmx 2 n+1 is injective (see Theorem 2.1 and Corollary 7.9 of [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF]), we now that ∀t ∈ [0, t 0 ], e -t(Pp m+1,t ) w -tcmx 2 n+1 φ L 2 > 0.

So we have

ε t = m j=1 e -t(Pp j,t ) w -tcx 2 n+1 φ, e -t(Pp m+1,t ) w -tcmx 2 n+1 φ L 2 e -t(Pp m+1,t ) w -tcmx 2 n+1 φ -2 L 2 .
But it follows for the Theorem 4.2 of [START_REF] Hörmander | Symplectic classification of quadratic forms, and general Mehler formulas[END_REF] that t → e -t(Pp 1,t ) w -tcx 2 n+1 . . . e -t(Ppm,t) w -tcx 2 n+1 φ and t → e -t(Pp m+1,t ) w -tcmx 2 n+1 φ are continuous from [0, t 0 ] to S(R n ). Consequently, t → ε t is continuous as product of two continuous functions and we have ε t ≡ 1, i.e. [START_REF] Rotman | An introduction to the Theory of Groups[END_REF] ∀t ∈ [0, t 0 ], e -t(Pp 1,t ) w -tcx 2 n+1 . . . e -t(Ppm,t) w -tcx 2 n+1 = e -t(Pp m+1,t ) w -tcmx 2 n+1 .

To conclude, we just have to prove that this relation can be evaluated in x n+1 = 1.

Let φ ∈ S(R n ) and choose ψ ∈ S(R n+1 ) such that

ψ |R n ×{1} = φ.
Applying m + 1 times Lemma 3.10 to [START_REF] Rotman | An introduction to the Theory of Groups[END_REF], we deduce that ∀t ∈ [0, t 0 ], e -t(p 1,t ) w -tc . . . e -t(pm,t) w -tc φ = e -tp w m+1,t -tcm φ.

Thus, since it is clear that for all t ∈ [0, t 0 ] and all j ∈ 1, m we have e -t(p j,t ) w -tc = e -t(p j,t ) w e -tc and e -t(p m+1,t ) w -tmc = e -t(p m+1,t ) w e -tmc we deduce that ∀t ∈ [0, t 0 ], e -t(p 1,t ) w . . . e -t(pm,t) w φ = e -tp w m+1,t φ. Finally, since the semigroups are continuous on L 2 (R n ) and that the previous relation holds for all φ ∈ S(R n ) that is dense in L 2 (R n ), we get the desired factorization.

Applications

In this section, we apply the results of the previous sections in order to get some efficient exact splitting methods for the magnetic linear Schrödinger equations with quadratic potentials, some transport equations and some Fokker-Planck equations. Another paper [START_REF] Bernier | Exact splitting methods for kinetic and Schrödinger equations[END_REF], written in collaboration with Nicolas Crouseilles and Yingzhe Li, deals with the implementation of these methods and compares numerically their efficiency, their accuracy and their qualitative properties with respect to the existing methods. 4.1. Transport equations. As we have seen through the formula [START_REF] Dunford | Linear operators. Part I[END_REF] two dimensional rotations can be computed efficiently as products of shear transforms.

In fact, this kind of factorization is much more general since as a classical application of the Gaussian elimination algorithm, we know that each matrix of determinant 1 is a product of shear matrices 5 (see e.g. Lemma 8.7 of [START_REF] Rotman | An introduction to the Theory of Groups[END_REF]). Identifying as usual matrices with linear maps on R n , this factorization is written as follows.

Proposition 4.1. For all G ∈ SL n (R), there exists m ≥ 0, α ∈ R m , k, ℓ ∈ 1, n m such that for all j ∈ 1, m , k j = ℓ j and ∀u ∈ L 2 (R n ), u • G = exp(α 1 x k 1 ∂ x ℓ 1 ) • • • exp(α m x km ∂ x ℓm )u
Note that in the context of the exact splitting, it is more natural to apply this result for G = exp(tB) where B ∈ sl n (R) is a matrix whose trace vanishes.

If we focus on transforms associated with matrices of determinant different than 1, we have to deal with one dimensional dilatations. Indeed, it is clear that each matrix G ∈ GL n (R) can be factorized as a product of a matrix in SL n (R) and the diagonal matrix diag(1, . . . , 1, det G). The following proposition provides a way to deal with positive dilatations with pseudo-spectral methods. Proposition 4.2. For all λ > 0, we have

∀u ∈ L 2 (R), u(λ •) = λ -1/2 e -iε λ α λ ∂ 2 x e -iβ λ x 2 e iε λ β λ ∂ 2 x e iα λ x 2 u,
where

α λ = 1 2 |λ -1 -1|λ -1 , β λ = 1 2 |1 -λ|, ε λ = 1 if λ ≤ 1 and ε λ = -1 else.
Proof of Proposition 4.2. We observe that if t = log λ then u(λ •) = e tx∂x u = e t(ixξ-1 2 ) w u.

Consequently, this proposition is a consequence of the proposition 3.6 and an elementary formal computation which can be checked, for example, with a formal computation software.

We do not provide similar formulas when λ < 0 since we do not know if it may be useful for applications. Nevertheless, since we have (see e.g. Thm 2.2.3 of [START_REF] Nicola | Global Pseudo-Differential Calculus on Euclidean Spaces[END_REF])

(28) ∀u ∈ L 2 (R), u(-x) = -ie i π 2 (x 2 -∂ 2 x ) u(x),
we note that, as a consequence of Theorem 1.3, such formulas exist. The factorization provided by the Gaussian elimination algorithm in Proposition 4.1 is not, in general, the most efficient possible. The following proposition, that is an application of Theorem 2.2, provides some efficient exact splittings generalizing the exact splittings for rotations [START_REF] Dunford | Linear operators. Part I[END_REF]. Proposition 4.3. Let M be a real square matrix of size n ≥ 1 such that [START_REF] Welling | Rotation of 3D volumes by Fourier-interpolated shears[END_REF] ∀i, M i,i = 0 ∃i, ∀j = i, M i,j = 0 5 also called transvection matrices.

then there exist t 0 > 0 and an analytic function (y (ℓ) , (y (k) ) k =j , y (r) ) : (-t 0 , t 0 ) → R n×(n+1) satisfying (30) y

(ℓ) i = y (r) i = 0 ∀k = i, y (k) k = 0
such that for all t ∈ (-t 0 , t 0 ) we have

e tM x•∇ = e t(y (ℓ) (t)•x)∂x i   k =i e t(y (k) (t)•x)∂x k   e t(y (r) (t)•x)∂x i .
For example, this proposition justifies a priori the form of the exact splitting used in [START_REF] Chen | 3D volume rotation using shear transformation[END_REF][START_REF] Welling | Rotation of 3D volumes by Fourier-interpolated shears[END_REF] to compute three dimensional rotations.

Proof of Proposition 4.3. Let B := t M be the transpose of M . We are going to prove that there exists an analytic function (y (ℓ) , (y (k) ) k =i , y (r) ) : (-t 0 , t 0 ) → R n×(n+1) satisfying (30) such that for all t ∈ (-t 0 , t 0 ),

(31) e tB = (I n + t y (ℓ) (t) ⊗ e j )   k =j (I n + t y (k) (t) ⊗ e k )   (I n + t y (r) (t) ⊗ e j ) = e ty (ℓ) (t)⊗e j   k =i e ty (k) (t)⊗e k   e ty (r) (t)⊗e j .
Consequently, Proposition 4.3 will be proven since it is enough to apply the characteristic formula [START_REF] Dujardin | Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations[END_REF] to the transpose of (31). Note that the second equality in (31) comes form the fact that the matrices are nilpotent of order 1.

In order, to prove this factorization applying the Theorem 2.2, we realize the change of unknown y (i) = y (r) + y (ℓ) and thus the first factor of (31) becomes e ty (ℓ) (t)⊗e j = e -ty (r) (t)⊗e j e ty (i) (t)⊗e j .

Up to an irrelevant permutation of indices, without loss of generality we assume that i = 1. In order to apply the Theorem 2.2, we define

b k = {y ⊗ e k | y ∈ R n satisfies y k = 0}, s = b 1 and b = n k=1 b k and b ⋆,k = Be k ⊗ e k (i.e. b ⋆ = B).
Consequently, to prove this proposition applying the Theorem 2.2, we just have to prove that its assumption (8) is satisfied. In our context, we are going to prove that

(32) b + ad B (s) = sl n (R).
So, from now, we just focus on proving (32). But, since the inclusion " ⊂ " is obvious, we just have to prove that

(33) b ∩ ad B (s) = {0}, dim b + dim ad B (s) = dim sl n (R).
Let Ψ : s → R n be the linear map defined by

∀s ∈ s, Ψ(s) = diag • ad B (s)
where diag : gl n (R) → R n is the natural map extracting the diagonal coefficients of a matrix. Since b = Ker diag is the space of the matrices with diagonal coefficients are equal to zero, the first relation of ( 33) is equivalent to have Ker Ψ = {0}. Furthermore, if Ker Ψ = {0}, then ad B is injective on s and we have

dim ad B (s) = dim s = n -1.
Since, for all k, dim b k = n -1 and (n + 1)(n -1) = n 2 -1 = dim sl n (R), we deduce that if Ker Ψ = {0} then the second relation of (33) also holds.

Finally, we just have to verify that Ker Ψ = {0}. Let y ⊗ e 1 ∈ Ker Ψ where y ∈ R n satisfies y 1 = 0. By assumption, if j = 1 then 0 = t e j [B, y ⊗ e 1 ]e j = -t e j y t e 1 Be j = -y j B j,1 .

Since, by assumption, B j,1 = 0 for j = 1, we deduce that y = 0, i.e. Ker Ψ = {0}.

Schrödinger equations.

Let v : R n → R be a real quadratic form and B ∈ so n (R) be a skew symmetric matrix. We aim at solving the following linear Schrödinger equation on R n (34)

i∂ t u(t, x) + 1 2 ∆u(t, x) -v(x)u(t, x) + iBx • ∇u(t, x) = 0.
In this context, a diagonal quadratic form is defined as follow.

Definition 4.4. A quadratic form is diagonal on R n if its matrix in the canonical basis is diagonal.

The following theorem provides an optimized splitting method to solve (34). Its proof is given at the end of this subsection. Theorem 4.5. There exists some quadratic forms v (r) t , a t on R n , a strictly upper triangular matrix U t ∈ M n (R), a strictly lower triangular matrix

L t ∈ M n (R) a diagonal quadratic form v (ℓ)
t on R n , all depending analytically on t ∈ (-t 0 , t 0 ) for some t 0 > 0, such that for all t ∈ (-t 0 , t 0 ) we have

e it(∆/2-v(x))-tBx•∇ = e -itv (ℓ) t (x)   n-1 j=1 e -t(Utx) j ∂x j   e itat(∇)   n j=2 e -t(Ltx) j ∂x j   e -itv (r) t (x)
where a t (∇) denotes the Fourier multiplier of symbol -a t (ξ) and (U t x) j (resp. (L t x) j ) the j st coordinate of U t x (resp. L t x).

The efficiency of this method is optimal, more precisely it is as cheap as a basic Lie splitting. Indeed, considering that the computational cost of this kind of method is proportional to the number of one dimensional Fast Fourier Transform required to implement it, the method of Theorem 4.5 requires 2n 1d-FFTs 6 which is the same as the elementary Lie splitting

e it(∆/2-v(x))-tBx•∇ = t→0 e -itv(x) n j=1 e it 2 ∂ 2 x j -(Bx) j •∂x j + O(t 2 ).
Note that this method is more general and more efficient than the other existing exact splittings for (34). For example, the splitting (2) of Bader in [START_REF] Bader | Fourier-splitting methods for the dynamics of rotating Bose-Einstein condensates[END_REF], designed to solve (34) in dimension n = 2, requires 6 1d-FFTS.

It seems to the author that there does not exist any simple explicit formula to compute the coefficients of the exact splitting given in Theorem 4.5. However, due to some particular properties of nilpotency, a simple and efficient iterative method is available to compute them. Indeed, identifying a quadratic form with its symmetric matrix in order we define, if t is small enough, we define, by induction, the following sequences

         A t,k+1 = A t,k + I n /2 -A t,k L t,k+1 = L t,k + L -L t,k U t,k+1 = U t,k + U -U t,k V (m) t,k+1 = V (m) t,k + V -V (m) t,k + t 2 [D t,k , B] + t 2 2 D 2 t,k
where

(A t,0 , L t,0 + U t,0 , V (m)
t,0 ) = (I n /2, B, V ), L + U = B and 2 V (m) t,k t L t,k + t U t,k + tD t,k L t,k + U t,k + tD t,k 2 A t,k = -t -1 J log(P t,k )
and

P t,k =   n-1 j=1 I n + tU (j) t,k I n -t t U (j) t,k   I n 2tA t,k I n   n j=2 I n + tL (j) t,k I n -t t L (j) t,k   × I n -2tV (m) t,k I n with L (j) t,k = (e j ⊗ e j )L t,k , U (j)
t,k = (e j ⊗ e j )U t,k and (e 1 , . . . , e n ) the canonical basis of R n . Theorem 4.6. There exists τ 0 ∈ (0, t 0 ) such that if 0 < |t| < τ 0 then the preceding sequences are well defined and

|A t -A t,k | + |L t -L t,k | + |U t -U t,k | + |V (ℓ) t + 1 2 D t,k | + |V (r) t -V (m) t,k - 1 2 D t,k | ≤ t τ 0 k .
We could apply directly Theorem 2.2 to prove Theorem 4.5. Indeed, by Proposition 3.6, it is clearly enough to prove that there exist some quadratic forms v (r) t , a t on R n , a strictly upper triangular matrix U t ∈ M n (R), a strictly lower triangular matrix L t ∈ M n (R) and a diagonal quadratic form v (ℓ) t on R n , all depending analytically on t ∈ (-t 0 , t 0 ) for some t 0 > 0, such that for all t ∈ (-t 0 , t 0 ) (35)

Φ i( |ξ| 2 2 +v(x)+Bx•ξ) t = Φ iv (ℓ) t (x) t   n-1 j=1 Φ i(Utx) j ξ j t   Φ iat(ξ) t   n j=2 Φ i(Ltx) j ξ j t   Φ iv (m) t (x) t Φ -iv (ℓ) t (x) t where v (r) t := v (m) t -v (ℓ)
t . Furthermore note that all the factors of (35) are some real symplectic transforms belonging to Sp 2n (R) and that all the right hand side flows are exponentials of nilpotent matrices (it explains why we don't have to compute exponential of matrices in the iterative method). Nevertheless, due to some convenient algebraic cancelations, the iterative method described below is not exactly the one associated with the proof of Theorem 2.2. Consequently, here it is easier to prove Theorem 4.5 and Theorem 4.6 simultaneously. However, the proofs being very similar to the proof of Theorem 2.2, so it is done more informally.

Before giving the proof, let us just mention shortly the algebraic decomposition we would use in order to satisfy the assumption (8) of Theorem 2.2 to prove (35). Recalling that sp 2n (R) is isomorphic to the Lie algebra of the real quadratic forms (equipped with the canonical Poisson bracket), we would simply have to observe that all real quadratic form q : R 2n

= R n x × R n ξ → R written q(x, ξ) = a (q) (ξ) + M (q) x • ξ + b (q) (x)
where a (q) , b (q) are real quadratic forms on R n and M (q) ∈ M n (R), can be decomposed as

q(x, ξ) = n-1 j=1 (U (q) t x) j ξ j + a (q) (ξ) + n j=2 (L (q) t x) j ξ j + b (q) (x) - [D, B]x • x 2 + |ξ| 2 2 + Bx • ξ + v(x), t xD (q) x 2 .
where M (q) = L (q) + D (q) + U (q) is the natural decomposition of M (q) .

Proof of Theorem 4.5 and Theorem 4.6. Similarly to the proof of Theorem 2.2, we prove that the exact splitting can be obtained from the resolution of a nonlinear equation via the Implicit Function Theorem. Consequently, Theorem 4.6 is nothing but the convergence of the natural iterative method associated with the Implicit Function Theorem.

In order to obtain (35), we consider (35) as a nonlinear equation where v

(m) t , a t , U t , L t , v (ℓ) t
are the unknowns. Applying the BCH formula, we get a family of real quadratic forms q t on R 2n (depending analytically on v

(m) t , a t , U t , L t ) such that if t is small enough then   n-1 j=1 Φ i(Utx) j ξ j t   Φ iat(ξ) t   n j=2 Φ i(Ltx) j ξ j t   Φ iv (m) t (x) t = Φ iqt t .
Thus, since Φ -iv

(ℓ) t (x) t
is symplectic, applying Proposition 3.3, we get

Φ iqt•Φ -iv (ℓ) t (x) t t = Φ i( |ξ| 2 2 +v(x)+Bx•ξ) t .
The exponential map begin injective near the origin, we deduce that, (35) is equivalent to (36)

q t • Φ -iv (ℓ) t (x) t = |ξ| 2 2 + v(x) + Bx • ξ. However Φ -iv (ℓ) t (x) t is a shear transform Φ -iv (ℓ) t (x) t (x, ξ) = (x, ξ + 2tV (ℓ) t ).
Consequently, decomposing naturally q t , with notations similar to the ones of the iterative method, as

q t (x, ξ) = v (m) t (x) + ( L t x + U t x + tD t )x • ξ + a t (ξ)
we deduce that (36) is equivalent to the system of equations

a t (ξ) = |ξ| 2 /2, L t = L, U t = U, V (ℓ) t = -D t /2, v (m) t (x) = v(x) -( L t x + U t x + tD t )x • 2tV (ℓ) t x -2t 2 |V (ℓ) t x| 2 .
Noting that by substitution, the last equation writes

v (m) t (x) = v(x) + t 2 [B, D t ]x • x + t 2 2 D t x • D t x,
we conclude that to get the factorization (35), it is enough to choose V

= -D t /2 and to solve the nonlinear equation

F (t, a t , L t , U t , v (m) t ) : = ( a t (ξ), L t , U t , v (m) t (x) - t 2 [B, D t ]x • x - t 2 2 D t x • D t ) = (|ξ| 2 /2, L, U, v(x)).
As in the proof of Theorem 2.2, this nonlinear equation can be solved by the Implicit Function Theorem. Indeed, F (0, |ξ| 2 /2, L, U, v(x)) = (|ξ| 2 /2, L, U, v(x)), F is an analytic function of t, a t , L t , U t , v (m) t , and the partial differential of F in (0, |ξ| 2 /2, L, U, v(x)) with respect to (a t , L t , U t , v (m) t

) is the identity. Note that the natural iterative method associated with the resolution of this nonlinear equation with the implicit function theorem is exactly the one introduced above. 4.3. Fokker-Planck equations. We apply our exact splitting methods to two Fokker-Planck equations. These equations can be used to describe particles system with collisions (in plasma physics or astrophysics). We refer to [START_REF] Hérau | Tunnel effect for the Kramers-Fokker-Planck type operators[END_REF][START_REF] Dujardin | Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations[END_REF][START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF] for more details about these equations.

4.3.1.

Exact splitting for the Fokker-Planck equation. First, we focus on the classical inhomogeneous Fokker-Planck equation

(FP) ∂ t u(t, x, v) + v∂ x u(t, x, v) = ∂ v (v + ∂ v )u(t, x, v).
Note that, since this equation comes from kinetic models, the variable are not denoted x 1 , x 2 but x, v. Implicitly, in this paper, the Fourier variable canonically associated with v is denoted by η.

Proposition 4.7. The following factorization provides an exact splitting for FP ∀t ≥ 0, e -t(v∂x-∂ 2 v -∂vv) = e t/2 e -(e t -1)v∂x e ∇•At∇ e iαt∂ 2 v e -iβtv 2 e -iβt∂ 2 v e iαtv 2 .

where

α t = 1 2 (1 -e -t )e -t , β t = 1 2 √ e t -1, ∇ = t (∂ x , ∂ v
) and A t is the positive matrix defined by

A t = 1 2 e 2t + 2t + 3 -4e t -4 sinh 2 (t/2) -4 sinh 2 (t/2) 1 -e -2t .
A priori, it could seem strange to have to solve some Schrödinger equations in order to compute the solution of FP. However a part of the Fokker-Planck operator is associated with a dilatation, i.e. v∂ v , and as explained in the previous subsection (see Proposition 4.2), the semigroup it generates cannot be computed only by compositions of shear transforms.

Proof. First, let us check that A t is nonnegative. Since the second diagonal coefficient of A t is positive, to check that A t is nonnegative, it is enough to prove that det A t ≥ 0 for t > 0. By a some elementary formal computations, we get

det A t = (3e 4t -12e 3t + (8t + 2)e 2t + 20e t -8t -13)e -2t
16 .

Consequently, we just have to prove that (16e 2t det A t ) ≥ 0. However, realizing the Taylor expansion of 16e 2t det A t in 0 (with a formal computation software), we get

16e 2t det A t = t→0 O(t 4 ).
Thus, it is enough to prove that

∂ 2 t (e -2t ∂ 2 t (16e 2t det A t )) ≥ 0,
which is obvious since

∂ 2 t (e -2t ∂ 2 t (16e 2t det A t )) = 192e 2t -108e t + 20e -t ≥ 0.
Actually, it would be possible to prove a priori that A t is nonnegative. Indeed, in the Fourier variables, the semigroup of FP is associated with a transport equation for which many elementary explicit computations can be done. In particular, it follows of a formula of Kolmogorov (see [START_REF] Kolmogorov | Zufällige Bewegungen (zur Theorie der Brownschen Bewegung)[END_REF][START_REF] Alphonse | Smoothing properties of fractional Ornstein-Uhlenbeck semigroups andnull-controllability[END_REF]) that A t admits an integral representation on which it is obvious to see that A t is nonnegative. Then, we consider the following factorization (whose form could be guessed using Theorem 2.2) that can be checked easily using a formal computation software ∀t ≥ 0, Φ

η 2 +i(vξ-vη) t = Φ i(e t -1)vξ 1 Φ (ξ,η)At t (ξ,η) 1 Φ -ivη t .
Since A t is nonnegative for t ≥ 0, applying Proposition 3.6, we deduce that ∀t ≥ 0, e -t(v∂x-∂ 2 v -∂vv) = e t e -(e t -1)v∂x e ∇•At∇ e tv∂v .

The last factor being a dilatation, we conclude this proof applying Proposition 4.2.

4.3.2.

Exact splitting for the Kramer-Fokker-Planck equation. Now we focus on the Kramer-Fokker-Planck equation

(KFP) ∂ t u(t, x, v) + v 2 u(t, x, v) -∂ 2 v u(t, x, v) + v∂ x u(t, x, v) = 0.
For some discussions about this equation, we refer to [START_REF] Hérau | Tunnel effect for the Kramers-Fokker-Planck type operators[END_REF][START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF].

Proposition 4.8. The following factorization provides an exact splitting for KFP

∀t ≥ 0, e -t(v 2 -∂ 2 v +v∂x) = e -1 2 tanh t v 2 e ∇•At∇ e -tanh t v∂x e -1 2 tanh t v 2
where A t is the nonnegative matrix defined by

A t = 1 2 α t sinh 2 t sinh 2 t sinh 2t , with α t = 1 2 (t -(tanh t)(1 -sinh 2 t)).
Proof. Applying Proposition 3.6, we just have to check that A t is non negative for t ≥ 0 and that

∀t ≥ 0, Φ v 2 +η 2 +ivξ t = Φ 1 2 tanh t v 2 1 Φ (ξ,η)At t (ξ,η) 1 Φ i tanh t vξ 1 Φ 1 2 tanh t v 2 1
.

This factorization, whose form can be deduced a priori using Theorem 2.2, can be verified by a formal computation software. Since the second diagonal coefficient of A t is positive, to check that A t is non negative, it is enough to prove that det A t ≥ 0 for t > 0. Thus, we conclude this proof observing that by Jensen's inequality and some elementary trigonometric computations In this basis the matrix of J 2n+2 is

det A t =
mat B J 2n+2 = J 2n J 2
and the matrix of Pp j is

mat B Pp j =   Q (j) 1 2 t L (j) 1 2 L (j) c j 0   ,
Consequently, in this basis the matrix of the Hamiltonian map is

mat B J 2n+2 Pp j =   J 2n Q (j) 1 2 Υ q j J 2n t L (j) 0 -1 2 L (j) -c j 0   .
Here, it is really relevant to observe the double triangular structure of this matrix. The four block on the top left corner defines an upper triangular matrix by blocks, whereas, considering these four blocks as a single block, the global matrix is lower triangular by blocks. Now observe, through the power expansion series, that if Ψ is an entire function and

M = A B 0
is an upper triangular matrix by blocks then

Ψ(M ) = Ψ(A) Ψ(z)-Ψ(0) z (A)B Ψ(0)
.

Consequently, we have

mat B Φ Pp j 1 =   Φ q j 1 -iΥ q j J 2n t L (j) 1 iL (j) Υ q j 2i κ j + 2ic j 1   where κ j = - i 2 L (j) Θ q j J 2n t L (j) with Θ q j = e z -1 -z z 2 (-2iJ 2n Q (j) ).
At the end of the proof, we will check that κ j = κ j , so for the moment assume that this relation holds.

Thus, realizing a product by block we get by a straightforward induction

mat B m j=1 Φ Pp j 1 =          m j=1 Φ q j 1 -i m j=1   k<j Φ q k 1   Υ q j J 2n t L (j) 1 i m j=1 L (j) Υ q j k>j Φ q k 1 2i m j=1 κ j + c j + σ j 1         
.

Identifying the blocks (1, 1), (3, 1), (3, 2) with those of mat B Φ Pp m+1 1 we get the system [START_REF] Paeth | A fast Algorithm for General Raster Rotation[END_REF]. Conversely, we have to check that if [START_REF] Paeth | A fast Algorithm for General Raster Rotation[END_REF] is satisfied then the blocks (1, 2) are the same. Indeed, consider a complex symplectic matrix M ∈ Sp 2n (C) with a block structure of the form

mat B M =   A B 1 C d 1   .
Note that since M is symplectic, M is invertible and consequently A is also invertible. Since M is symplectic, then is satisfies

(38) M = -J 2n t M -1 J 2n .
However, due to the double triangular nature of M , its invert can be computed easily and we have

mat B M -1 =   A -1 -A -1 B 1 -CA -1 CA -1 B -d 1   .
Consequently, a straightforward block product leads to

-mat B J 2n t M -1 J 2n =   -J 2n t A -1 J 2n -J 2n t A -1 t C 1 -t B t A -1 J 2n d -CA -1 B 1   .
Thus, since M is symplectic, we have

B = -J 2n t A -1 t C.
A fortiori, if two symplectic matrices have this block structure and the same top left corner blocks, if their blocks (3, 1) are equal then their blocks (1, 2) are equal. Consequently, applying this to the symplectic matrices Φ Pp , we deduce that if the system ( 24) is satisfied then we have the factorization [START_REF] Nicola | Global Pseudo-Differential Calculus on Euclidean Spaces[END_REF].

Finally, we just have to check that κ j = κ j . For this computation, we omit the indices j since they are clearly irrelevant. First, we split the even indices from the odd indices in the power expansion of κ :

2i κ = LΘ q J 2n t L = k∈N 1 (k + 2)! L(-2iJ 2n Q) k J 2n t L = k∈N 1 (2k + 2)! L(-2iJ 2n Q) 2k J 2n t L + k∈N 1 (2k + 3)! L(-2iJ 2n Q) 2k+1 J 2n t L := Σ even + Σ odd .
Observing that

(J 2n Q) 2k J 2n = (J 2n Q) k J 2n (QJ 2n ) k , since J 2n is skew-symmetric we have L(J 2n Q) 2k J 2n t L = (-1) k L(J 2n Q) k J 2n t L(J 2n Q) k = 0.
Consequently, Σ even vanishes. Similarly, since we have

(J 2n Q) 2k+1 J 2n = (J 2n Q) k J 2n Q(J 2n Q) k J 2n = (-1) k+1 (J 2n Q) k J 2n Q t ((J 2n Q) k J 2n ),
we get for t ∈ R, j, k ∈ 1, n . Applying Theorem 3.5, if we prove that G = Sp 2n (R), we deduce that L is a product of operators of the form e itx 2 j , e it∂ 2

Σ odd = 2i k∈N 4 k (2k + 1)! L(J 2n Q) k J 2n Q t ((J 2n Q) k J 2n ) t L = 2iκ. 5
x j , e tx j ∂x k (up to the sign). Thus, since Proposition 4.2 states that dilatations (i.e. operators of the form e tx k ∂ k ) can be factorized similarly, we would deduce that L can be computed by an exact splitting.

Consequently, we aim at proving that G = Sp 2n (R). First, let us prove that G contains a neighborhood of the identity in Sp 2n (R).

Indeed, consider the map Ψ : N → sp 2n (R), where N is a neighborhood of the origin in

sp 2n (R), defined for Q ∈ S 2n (R) such that J 2n Q ∈ N by (39) Ψ(J 2n Q) = log   n j=1 Φ iA j,j x 2 j 1/2 1≤j<k≤n Φ iA j,k x j x k 1 n j=1 n k=1 Φ iC j,k x j ξ k 1 1≤j<k≤n Φ iB j,k ξ j ξ k 1 n j=1 Φ iB j,j ξ 2 j 1/2   ,
where the natural block decomposition of Q is

Q = A C t C B with A, B ∈ S n (R) and C ∈ M n (R).
Note that to prove that Ψ takes its values in sp 2n (R), it is enough to apply the Baker-Campbell-Hausdorff formula.

Since the differential of the exponential in the origin and the differential of the logarithm in the identity are equal to the identity, we deduce by composition that the differential of Ψ in the origin is also the identity. Thus, since Ψ vanishes in the origin, we deduce of the Inverse Function Theorem that Ψ defines a local homeomorphism around the origin. Furthermore, we recall 7 that the exponential is an homeomorphism between a neighborhood of the origin in sp 2n (R) and a neighborhood of the identity in Sp 2n (R). Consequently, we deduce that each matrix in Sp 2n (R) close enough to the identity can be written as a product of the form of the product in the logarithm in (39). A fortiori, we have proven that G contains a neighborhood of the identity in Sp 2n (R). Let V denotes this neighborhood.

Since Sp 2n (R) is connected (see e.g. the subsection 4.4 of [START_REF] Arnold | Dynamical Systems IV, Symplectic Geometry and its Applications[END_REF]), to prove that G = Sp 2n (R), we just have to prove that G is closed and open in Sp where q is the quadratic part of p.

2n (R). Indeed, if g ∈ G then gV is a neighborhood of g in Sp 2n (R) and since V is included in G then gV is also included in G. Thus G is open in Sp 2n (R). Conversely, if g / ∈ G then gV is also a neighborhood of g in Sp 2n (R) but since G is a group we have gV ∩ G = ∅. Consequently, the complementary of G in Sp 2n (R) is open, i.e. G is closed in Sp 2n (R), which conclude the proof.
Proof. Let l be the linear part of p. Considering the natural action of the entire functions on M 2n (C), let ℓ = ℓ 1 and c = c 1 where

ℓ t = l • e z -1 z (-2tJ 2n Q) and c t = p(0) -t 2 k∈N (-4) k t 2k (2k + 3)! q((J 2n Q) k J 2n t L)
with Q ∈ S 2n (R), the matrix of q and L the matrix of l. Applying Lemma 3.7, we have

∀t ∈ R, Φ -iPct t Φ -iPq t Φ -iPℓt t = Φ -iPp t .
Consequently, applying Proposition 3.6 at t = 1, we get (40).

Lemma 5.4. If p : R 2n → R is a real valued polynomial of degree 2 or less bounded below, then there exists a real valued linear form ℓ : R 2n → R and c ∈ R such that (41) e -p w = e -c e -iℓ w e -q w e iℓ w , where q is the quadratic part of p.

Proof. Applying Lemma 3.8, we get Y ∈ R 2n such that p = p 1 where p t = q(• -tY ) + p(Y ).

and X = (x 1 , . . . , x n , ξ 1 , . . . , ξ n ).

Let ℓ = -t (J 2n Y )X and c = p(Y ). Let B be the basis introduce in the proof of Lemma 3.7 and defined by (37). Observing that

mat B J 2n+2 ∇Pℓ =   x n+1 Y 0 ℓ   ,
we deduce that

(mat B Φ -iPℓ t )   X x n+1 ξ n+1   =   X -tx n+1 Y x n+1 ξ n+1 -tℓ   .
Consequently, we have P(q + c) Lemma 5.5. If q : R 2n → R is a non-negative real quadratic form then exp(-q w ) can be computed by an exact splitting.

Proof. Applying Theorem 21.5.3 of [START_REF] Hörmander | The analysis of linear partial differential operators. III, Pseudo-differential operators[END_REF], we get a symplectic change of coordinates T ∈ Sp 2n (R) such that there exists m ∈ 1, n and 0 < λ 1 ≤ • • • ≤ λ m some real numbers such that

q • T = m j=1 λ j (x 2 j + ξ 2 j ) + n j=m+1 x 2 j .
Consequently, since T is symplectic, applying Noether theorem, we have Recalling that, as a consequence of the formula (2), the semigroups e λ j (x 2 j -∂ 2

Φ q 1 = T Φ q•T 1 T -1 = T   m j=1 Φ x 2 j +ξ 2
x j ) can be computed by exact splittings, we deduce of Lemma 5.1 that the same applies for e -tq w . Now, we can prove Theorem 1.3.

Proof of Theorem 1.3. Let p be a polynomial of degree 2 or less on C 2n whose real part is bounded below on R 2n . We aim at proving that e -p w can be computed by an exact splitting in the sense of Definition 1.2.

Applying Corollary 3.9, we get a constant c ∈ R such that P(pc) ≥ 0. Then applying Theorem 2.1 and Lemma 3.10 of [START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF], we get t 0 > 0 and two real quadratic forms a t , b t : R 2n+2 → R depending analytically on t ∈ (-t 0 , t 0 ), a t being non-negative, and such that if |t| < t 0 then Φ at t Φ ibt t = Φ Note that since p (r)

t is the restriction of a t on the affine subspace {(x n+1 , ξ n+1 ) = (1, 0)}, it is also non-negative. Now applying Proposition 3.6, we deduce that if 0 ≤ t < t 0 we have (42) e -tp w = e -tc e -t(p (r) t ) w e -it(p (i) t ) w .

Let t ⋆ ∈ (0, t 0 ) be such that there exists n ∈ N satisfying t ⋆ = n -1 . Since (e -tp w ) t≥0 is a semigroup, we have (e -t⋆p w ) n = e -nt⋆p w = e -p w .

Consequently, if e -t⋆p w can be computed by an exact splitting then the same applies for e -p w . Furthermore, from the factorization (42), we deduce that if e -t⋆(p (r) t⋆ ) w and e -it⋆(p (i) t⋆ ) w can be computed by an exact splitting then the same applies for e -t⋆p w . Consequently, we just have to focus on these two semigroups. On the one hand, applying Lemma 5.3, Lemma 5.2, Lemma 5.1 and Theorem 3.5 (to justify that the semigroup generated by a quadratic differential operator is a Fourier Integral Operator), we deduce that exp(-it ⋆ (p (i) t⋆ ) w ) can be computed by an exact splitting. On the other hand, applying Lemma 5.4, Lemma 5.2 and Lemma 5.5, we deduce that exp(-t ⋆ (p (r) t⋆ ) w ) can be computed by an exact splitting. Lemma 3.4. We aim at proving that if ℜq ≥ 0 on R 2n then Φ q t ≡ e -2itJ 2n Q ∈ Sp + 2n (C), where Q is the matrix of q. First, we recall that from Proposition 3.3 we know that Φ q t is a symplectic transformation (i.e. Φ q t ∈ Sp 2n (C)). So we aim at proving that Φ q t is nonnegative, i.e.

Proof of

(43) ∀X ∈ C 2n , t X t Φ q t (-iJ 2n )Φ q t X -t X(-iJ 2n )X ∈ R + .

Since Φ q 0 = I 2n , we have t Φ q t (-iJ 2n )Φ q t -(iJ 2n ) = Φ q s (ℜQ)Φ q s ds.

Since ℜQ is a real symmetric nonnegative matrix, t Φ q s (ℜQ)Φ q s is a Hermitian nonnegative matrix and thus t Φ q t (-iJ 2n )Φ q t -(iJ 2n ) is also a Hermitian nonnegative matrix which proves (43).

Theorem 2 . 2 .

 22 Let m be a positive integer and b 1 , . . . , b m , s be some subspaces of a real Lie algebra of matrices such that b 1 , . . . , b m are complementary. If b

Lemma 5 . 2 .e iL j x j n j=1 e L j ∂x j . Lemma 5 . 3 .

 52j=153 If ℓ : R 2n → R is a real linear form then exp(iℓ w ) be computed by an exact splitting.Proof. Let L ∈ R 2n be the matrix of ℓ and let c = c 1 wherec t = t 2 n j=1 L j L j+n .7 the reader could refer, for example, to Lemma 3.10 in[START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF] for a detailed proof. Proposition 3.6 at t = 1, we get e iℓ w = e ict n j=1 If p : R 2n → R is a real valued polynomial of degree 2 or less then there exists a real linear form ℓ : R 2n → R and c ∈ R, such that (40) e ip w = e ic e iq w e iℓ w ,

1  T - 1 .

 11 Applying Theorem 3.5, we get, at the level of the Fourier Integral Operators, ±e -tq w = K (T )

e 0 e

 0 -2is QJ 2n (-iJ 2n )e -2isJ 2n Q ds = t -2is QJ 2n (-2i QJ 2n )(-iJ 2n ) + (-iJ 2n )(-2iJ 2n Q) e -2isJ 2n Q ds

  W t has an unique fix point in B ⋆ and the sequence (b

	(k) t , s

we deduce of the mean value inequality that there exist t 2 ∈ (0, t 0 ) and a ball B ⋆ , for the norm • ⋆ , centered in (b ⋆ , s ⋆ ) such that for all t ∈ (-t 2 , t 2 ), W t is 1/2 Lipschitzian on B ⋆ , and W t B ⋆ ⊂ B ⋆ .

Consequently,

  Proof of Lemma 3.7. Here, the proof relies essentially on computations by block requiring to introduce a more convenient basis than the canonical basis of C 2n+2 denoted by e 1 , . . . , e 2n+2 . This basis, denoted B, is just a permutation of the canonical basis and is defined by (37) B = (e 1 , . . . , e n , e n+2 , . . . , e 2n , e 2n+1 , e n+1 , e 2n+2 ).

	1 4	tanh(t)(t -tanh(t)) ≥ 0.
		5. Appendix
	5.1.

  .2. Proof of Theorem 1.3. Before proving Theorem 1.3, let us to prove some preparatory lemmas. Let G be the group generated by Φ

	ix 2 j t , Φ	iξ 2 j t , Φ	ix j ξ k t

Lemma 5.1. If L is a bounded operator on L 2 (R n ) such that there exists T ∈ Sp 2n (R) satisfying L = ±K (T ), then L can be computed by an exact splitting.

Proof.

  Now observing that the Hamiltonian Pc commutes (i.e. for the canonical Poisson bracket) with all the other Hamiltonians, we deduce of the Noether theorem (or of Lemma 3.7) that Consequently, applying Proposition 3.6 at t = 1, we get (41).

	However, Φ -iPℓ t	is a symplectic map, so we have
		Φ iPℓ t Φ	P(q+c) t	Φ -iPℓ t	= Φ	P(q+c)•Φ -iPℓ t t	= Φ Ppt t .
		∀t ∈ R, Φ Pc t Φ iPℓ t Φ Pq t Φ -iPℓ t	= Φ Ppt t .

• Φ -iPℓ t = Pp t .

.

  Furthermore, it follows of the Baker-Campbell-Hausdorff formulas and formulas of (3.22) and (3.41) of[START_REF] Alphonse | Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects[END_REF] defining a t , b t that these quadratic forms belong to the complex Lie algebra generated by P(pc) and P(pc). Observing that for all polynomials p 1 , p 2 of degree 2 or less on C 2n we have{Pp 1 , Pp 2 } = P{p 1 , p 2 },the image of P is a Lie algebra. Consequently, a t , b t belong to the image of P, i.e. there exist two real polynomials of degree 2 or less on R 2n depending analytically on t and denoted p

	(r) t and p	(i) t such that
		a t = Pp	(r) t and b t = Pp	(i)
			P(p-c)
			t

t .

i.e. a set equipped with an associative binary operation and an identity element.

Note that as usual this formalism include nonlinear equations, since it is enough to consider the transport equations they generate.

in some specific sense depending on the problem.

An inverse Fourier transform has the same cost as a direct one. For example, to compute the solution of the semigroup generated by the harmonic oscillator with the factorization (2), we need 2 1d-FFTs, one to go on the Fourier side and one other to come back.