Evelyn Buckwar
email: evelyn.buckwar@jku.at

Adeline Samson
email: adeline.leclercq-samson@univ-grenoble-alpes.fr

Massimiliano Tamborrino
email: massimiliano.tamborrino@warwick.ac.uk

Irene Tubikanec
email: irene.tubikanec@jku.at

A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model

Keywords: Stochastic differential equations, Locally Lipschitz drift, Hypoellipticity, Ergodicity, FitzHugh-Nagumo model, Splitting methods, Mean-square convergence AMS subject classifications 60H10, 60H35, 65C20, 65C30

In this article, we construct and analyse an explicit numerical splitting method for a class of semi-linear stochastic differential equations (SDEs) with additive noise, where the drift is allowed to grow polynomially and satisfies a global one-sided Lipschitz condition. The method is proved to be mean-square convergent of order 1 and to preserve important structural properties of the SDE. First, it is hypoelliptic in every iteration step. Second, it is geometrically ergodic and has an asymptotically bounded second moment. Third, it preserves oscillatory dynamics, such as amplitudes, frequencies and phases of oscillations, even for large time steps. Our results are illustrated on the stochastic FitzHugh-Nagumo model and compared with known mean-square convergent tamed/truncated variants of the Euler-Maruyama method. The capability of the proposed splitting method to preserve the aforementioned properties may make it applicable within different statistical inference procedures. In contrast, known Euler-Maruyama type methods commonly fail in preserving such properties, yielding ill-conditioned likelihood-based estimation tools or computationally infeasible simulation-based inference algorithms.

Introduction

The aim of this article is to construct and analyse a splitting method for semi-linear stochastic differential equations (SDEs) of additive noise type dX(t) = F (X(t))dt + ΣdW (t) := AX(t) + N (X(t)) dt + ΣdW (t), X(0) = X 0 ,

where the diffusion matrix Σ may be degenerate and the drift F satisfies a global one-sided Lipschitz condition and is allowed to grow polynomially. Coefficients with these properties appear in many applications [START_REF] Hutzenthaler | Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF], ranging from physics [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF] over population growth problems [START_REF] Hutzenthaler | Ergodic behavior of locally regulated branching populations[END_REF][START_REF] Khasminskii | Long term behavior of solutions of the Lotka-Volterra system under small random perturbations[END_REF] to neuroscience [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF] and others. As an illustrative equation from this class of SDEs, we discuss the stochastic FitzHugh-Nagumo (FHN) model [START_REF] Berglund | Mixed-mode oscillations and interspike interval statistics in the stochastic FitzHugh-Nagumo model[END_REF][START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF][START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF], a well-known neuronal model describing the generation of spikes of single neurons at the intracellular level. This model is given by the 2-dimensional SDE

d   V (t) U (t)   :=X(t) =   1 V (t) -V 3 (t) -U (t) γV (t) -U (t) + β   := F (X(t)) dt +   σ 1 0 0 σ 2   :=Σ dW (t), (2)
where σ 1 may be zero. The V -component of the system describes the evolution of the membrane voltage of the neuron and the U -component is a recovery variable. Our aim is to construct a numerical method for [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF], and (2) in particular, which is easy to implement and also applicable across different disciplines in the broad field of statistical inference. This implies that the method needs to meet several requirements:

• Statistical applications require strong approximations of SDEs. Thus, we focus on the concept of mean-square convergence [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF][START_REF] Milstein | Stochastic Numerics for Mathematical Physics. Scientific computation[END_REF][START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF]. Since it was shown in [START_REF] Hutzenthaler | Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF] that the standard Euler-Maruyama method does not converge in the mean-square sense under the above assumptions on the drift F , the development of mean-square convergent variants of this method has received much attention. In particular, tamed [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients[END_REF][START_REF] Sabanis | Euler approximations with varying coefficients: the case of superlinearly growing diffusion coefficients[END_REF][START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF][START_REF] Zhang | Order-preserving strong schemes for SDEs with locally Lipschitz coefficients[END_REF] and truncated [START_REF] Guo | The partially truncated Euler-Maruyama method and its stability and boundedness[END_REF][START_REF] Hutzenthaler | Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF][START_REF] Mao | The truncated Euler-Maruyama method for stochastic differential equations[END_REF][START_REF] Mao | Convergence rates of the truncated Euler-Maruyama method for stochastic differential equations[END_REF] Euler-Maruyama methods have been proposed. They all aim to control the unbounded growth arising from the non-globally Lipschitz drift by enforcing a rescaling modification to the drift and/or diffusion coefficients.

• Simulation-based statistical methods require to generate paths of SDEs as computationally efficient as possible, see, e.g., [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF]. Using explicit numerical methods is a first step to achieve sufficiently low computational cost. While the aforementioned mean-square convergent variants of the Euler-Maruyama method are explicit, they commonly fail in preserving important structural properties of the SDE. The major key to computational efficiency, however, is to construct explicit methods which are capable to preserve the underlying properties for time steps as large as possible. This leads to the next point.

• An important issue in the field of (stochastic) numerical analysis is the preservation of structural properties of the considered SDE by the numerical method used to approximate it. Geometric Numerical Integration is a well-established framework in this context [START_REF] Hairer | Geometric Numerical Integration[END_REF]. Here, we discuss the preservation of hypoellipticity, geometric ergodicity and oscillatory dynamics such as amplitudes, frequencies and phases of oscillations:

-The diffusion matrix Σ of SDE (1) may be of full rank or degenerate, where in the latter case the SDE may be hypoelliptic, depending on the drift F . The case of degenerate noise naturally occurs in many applications [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF][START_REF] Leimkuhler | Molecular Dynamics: With Deterministic and Stochastic Numerical Methods[END_REF][START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF][START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF], with the hypoelliptic property ensuring that the solution of the SDE admits a smooth transition density [START_REF] Nualart | The Malliavin Calculus and Related Topics[END_REF]. This means that the noise is propagated through the whole system via the drift of the SDE, even though it does not directly act on all components. In many inference approaches using discrete approximations of SDEs, it is necessary that a discrete analogue of the hypoelliptic property holds at each iteration step. In particular, considering a discretised time interval with equidistant time steps ∆ = t i -t i-1 , the distribution of the numerical solution X(t i ) of (1) at time t i given the previous value X(t i-1 ) must admit a smooth density, a property that we term 1-step hypoellipticity. It is known that Euler-Maruyama type methods do not satisfy this if the SDE is not elliptic but only hypoelliptic. Thus, they yield ill conditioned likelihood-based inference methods [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF][START_REF] Melnykova | Parametric inference for hypoelliptic ergodic diffusions with full observations[END_REF][START_REF] Pokern | Parameter estimation for partially observed hypoelliptic diffusions[END_REF]. Higher-order Taylor approximation methods [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] may be 1-step hypoelliptic [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF]. However, since such methods may neither be mean-square convergent in the case of superlinearly growing coefficients [START_REF] Hutzenthaler | Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF], nor preserve other structural properties, they may lead again to ill-posed statistical problems. -The analysis of the asymptotic behaviour of the process is of further crucial interest. In particular, if SDE [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF] possesses an underlying Lyapunov structure, it may be geometrically ergodic [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. This property ensures that the distribution of the process converges exponentially fast to a unique limit for any starting value X 0 , and has two important statistical implications. First, the choice of the initial value X 0 is negligible since its impact on the distribution of the process decreases exponentially fast. This is relevant, especially when the process is only partially observed, since X 0 is usually not known. Second, there is a correspondence of "time averages along trajectories" and "space averages across trajectories" of geometrically ergodic systems, see, e.g., [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Da Prato | Ergodicity for infinite dimensional systems[END_REF]. This means that quantities related to the distribution of the process can be estimated from a single path simulated over a sufficiently large time horizon instead of relying on repeated simulations of trajectories. For the importance of this feature in statistical inference algorithms, we refer, e.g., to [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF]. Euler-Maruyama type methods may not provide these features as they tend to lose the Lyapunov structure of the SDE [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. In particular, here we illustrate that they react sensitively to the initial condition X 0 , and that they may yield poor approximations of the underlying invariant density of the process. -The last structural properties we are focusing on are features linked to oscillatory dynamics such as amplitudes, frequencies and phases of oscillations. Already in the deterministic scenario it has been observed that Euler type methods may not preserve amplitudes and frequencies of oscillations, see, e.g., [START_REF] Chen | Structure-preserving numerical integrators for Hodgkin-Huxley-type systems[END_REF][START_REF] Hairer | Geometric Numerical Integration[END_REF]. Similar findings have been made for Euler-Maruyama type methods in the stochastic case. For example, it has been proved that the Euler-Maruyama method does not preserve the growth rate of the second moment of linear stochastic harmonic oscillators, overshooting the amplitudes of the underlying oscillations, even for arbitrarily small time steps ∆ [START_REF] Strømmen Melbø | Numerical simulation of a linear stochastic oscillator with additive noise[END_REF]. Similar nonpreserving results of oscillation amplitudes have been observed for non-linear, ergodic and higher-dimensional stochastic oscillators [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Chevallier | Diffusion approximation of multi-class Hawkes processes: Theoretical and numerical analysis[END_REF][START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF]. Taming/truncating perturbations do not improve this behaviour. Even worse, taming perturbations may also lead to a non-preservation of frequencies of oscillations [START_REF] Kelly | Adaptive time-stepping strategies for nonlinear stochastic systems[END_REF][START_REF] Kelly | Adaptive Euler methods for stochastic systems with non-globally Lipschitz coefficients[END_REF]. This lack of amplitude and frequency preservation is also confirmed by our numerical experiments on the FHN model [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF]. Moreover, we find that Euler-Maruyama type methods may also not preserve phases of oscillations. This poor behaviour is linked to the non-preservation of geometric ergodicity.

Here, we propose to apply the splitting technique, an approach that addresses all previously listed issues. The general idea of this method is to split the SDE of interest into exactly solvable subequations, to derive their solutions, and to compose them in a suitable way. We refer to [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF][START_REF] Mclachlan | Splitting methods[END_REF] for a thorough discussion of splitting methods for ordinary differential equations (ODEs) and to [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF][START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Alamo | A technique for studying strong and weak local errors of splitting stochastic integrators[END_REF][START_REF] Bou-Rabee | Cayley splitting for second-order Langevin stochastic partial differential equations[END_REF][START_REF] Bréhier | Analysis of some splitting schemes for the stochastic Allen-Cahn equation[END_REF][START_REF] Chevallier | Diffusion approximation of multi-class Hawkes processes: Theoretical and numerical analysis[END_REF][START_REF] Leimkuhler | Molecular Dynamics: With Deterministic and Stochastic Numerical Methods[END_REF][START_REF] Milstein | Quasi-symplectic methods for Langevin-type equations[END_REF][START_REF] Misawa | A Lie algebraic approach to numerical integration of stochastic differential equations[END_REF][START_REF] Petersen | A general implicit splitting for stabilizing numerical simulations of Itô stochastic differential equations[END_REF][START_REF] Shardlow | Splitting for dissipative particle dynamics[END_REF] for articles considering extensions to SDEs. Note that it is often possible to split the differential equation under consideration into different sets of subequations, the choice of the most useful set depending on the problem to be solved. For the class of SDEs with additive noise, where the drift F consists of a linear and a non-linear term, i.e., F (X(t)) = AX(t)+N (X(t)), as in Equation ( 1), the idea is to exclude the nonlinear term N (X(t)) into a deterministic differential subequation and to treat the remaining linear term AX(t) via a stochastic differential subequation. In this article, we exploit properties of the exact solutions of both subequations, as illustrated on the FHN model [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF]. When N is globally Lipschitz continuous and uniformly bounded, this idea has been applied to the Jansen and Rit neural mass model in [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF], and for locally Lipschitz N it has been applied to the Allen-Cahn equation in [START_REF] Bréhier | Analysis of some splitting schemes for the stochastic Allen-Cahn equation[END_REF].

We illustrate that splitting methods may be able to deal with mean-square convergence issues arising from superlinearly growing coefficients. In particular, in Lemma 2, we prove the boundedness of the moments of the proposed splitting method. This result is the key to establish its mean-square convergence. In Theorem 2, we use Tretyakov's and Zhang's "Fundamental meansquare convergence theorem for SDEs with locally Lipschitz coefficients" [START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF] to prove that the splitting method converges with mean-square order 1. This is in agreement with the convergence rate of comparable known splitting methods in the globally Lipschitz scenario [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Alamo | A technique for studying strong and weak local errors of splitting stochastic integrators[END_REF][START_REF] Milstein | Quasi-symplectic methods for Langevin-type equations[END_REF], and with standard methods such as the Euler-Maruyama method in the case of additive noise [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. Moreover, we address the fact that splitting methods may also be able to tackle the problems arising from degenerate noise structures. In particular, in Theorem 3, we show that the proposed splitting method is 1-step hypoelliptic and yields non-degenerate multivariate normal transition distributions for any time step ∆, provided that the stochastic subequation of the splitting framework is chosen to be hypoelliptic. This may be beneficial for likelihood-based inference. Furthermore, in Theorem 4, we prove that the constructed method satisfies a discrete Lyapunov condition, and is thus geometrically ergodic, for any time step ∆. This result requires an assumption on the solution of the deterministic subequation defined via N (X(t)) and that e A∆ < 1, where the matrix norm is induced by the Euclidean norm. Moreover, in Corollary 1, we show that the second moment of the splitting method is asymptotically bounded by a constant which is independent of the time step size ∆ and the number of time steps i. This result holds if, in addition, the logarithmic norm [START_REF] Söderlind | The logarithmic norm. History and modern theory[END_REF][START_REF] Ström | On logarithmic norms[END_REF] of the matrix A is strictly negative. In the one-dimensional case, some of the involved expressions simplify such that, in Corollary 2, we obtain a precise closed-form (asymptotic) bound of the second moment of the proposed splitting method. This bound is illustrated on a cubic one-dimensional model problem with drift given by F (X(t)) = -X 3 (t) [START_REF] Hutzenthaler | Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF][START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. In addition, we illustrate the proposed splitting method on the stochastic FHN model ( 2) and show through a variety of numerical experiments that it preserves the qualitative dynamics of neuronal spiking, in particular, amplitudes, frequencies and phases of the underlying oscillations even for large time steps ∆.

The article is organised as follows. In Section 2, we introduce necessary mathematical preliminaries and notations, and we discuss equations of interest and relevant properties. In Section 3, we present the proposed splitting method. In Section 4, we establish its mean-square convergence. In Section 5, we prove its 1-step hypoellipticity, establish its geometric ergodicity, derive an (asymptotic) second moment bound and illustrate these results on a one-dimensional cubic model problem. In Section 6, we apply the proposed splitting approach to the stochastic FHN model [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF]. In Section 7, we provide a variety of numerical experiments, illustrating the theoretical results and reporting comparisons with different tamed/truncated variants of the Euler-Maruyama method. Conclusions are given in Section 8.

Model and properties

Throughout, the following notations are used. Notation 1. Let x, y ∈ R d be two generic vectors. Then x l denotes the l-th entry of x, x the transpose of x, x = (x 2 1 +. . .+x 2 d ) 1/2 the Euclidean norm of x and (x, y) = x 1 y 1 + . . . + x d y d the scalar product of x and y. Further, let A, B ∈ R d×d be two generic matrices. Then a lj denotes the component in the l-th row and j-th column of A, A the transpose of A, 0 d the ddimensional zero vector and I d the d × d-dimensional identity matrix. Moreover, we denote by A = λ max (A A) the matrix norm which is induced by the Euclidean norm, where λ max (A) is the largest eigenvalue of A, and with µ(A) = λ max ((A + A )/2) the real-valued logarithmic norm which results from the Euclidean norm and its induced matrix norm.

Let (Ω, F, P) be a complete probability space with filtration (F(t)) t∈[0,T ] . Further, let (W (t)) t∈[0,T ] be a m-dimensional Wiener process defined on that space and adapted to (F(t)) t∈[0,T ] . We consider the d-dimensional autonomous SDE of additive noise type ( 1)

dX(t) = F (X(t))dt + ΣdW (t) := AX(t) + N (X(t)) dt + ΣdW (t), X(0) = X 0 , where t ∈ [0, T ], T > 0, A ∈ R d×d , Σ ∈ R d×m , F : R d → R d and N : R d → R d are locally
Lipschitz continuous. The initial value X 0 is an F(0)-measurable R d -valued random variable which is independent of (W (t)) t∈[0,T ] and such that E X 0 2p < ∞ for all p ≥ 1.

Conditions required to ensure the existence of a unique strong solution of SDE [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF], which is regular in the sense of [START_REF] Khasminskii | Stochastic Stability of Differential Equations[END_REF], i.e., it is defined on the entire interval [0, T ] such that sample paths do not blow up to infinity in finite time, are, e.g., discussed in [START_REF] Alyushina | Euler polygonal lines for Itô equations with monotone coefficients[END_REF][START_REF] Khasminskii | Stochastic Stability of Differential Equations[END_REF][START_REF] Krylov | A simple proof of the existence of a solution of Itô's equation with monotone coefficients[END_REF][START_REF] Mao | Stochastic Differential Equations and Applications[END_REF]. Here, we follow the setting in [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients[END_REF][START_REF] Kelly | Adaptive time-stepping strategies for nonlinear stochastic systems[END_REF][START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF] and suppose that the drift satisfies a global one-sided Lipschitz condition and is allowed to grow polynomially at infinity. It suffices to place these conditions on N : Assumption 1. (A1) The function N is globally one-sided Lipschitz continuous, i.e., there exists a constant c 1 > 0 such that

(x -y, N (x) -N (y)) ≤ c 1 x -y 2 , ∀ x, y ∈ R d .
(A2) The function N grows at most polynomially, i.e., there exist constants c 2 > 0 and χ ≥ 1 such that

N (x) -N (y) 2 ≤ c 2 (1 + x 2χ-2 + y 2χ-2 )||x -y|| 2 , ∀ x, y ∈ R d .
Assumption 1 also ensures the finiteness of the moments of the solution of (1) [START_REF] Higham | Strong convergence of Euler-type methods for nonlinear stochastic differential equations[END_REF][START_REF] Khasminskii | Stochastic Stability of Differential Equations[END_REF][START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF][START_REF] Zhang | Order-preserving strong schemes for SDEs with locally Lipschitz coefficients[END_REF]. In particular, there exists a constant K(T, p) > 0 such that

E sup 0≤t≤T X(t) 2p ≤ K(T, p) 1 + E X 0 2p . (3)
Moreover, the process (X(t)) t∈[0,T ] is a Markov process. Denoting B(R d ) the Borel sigmaalgebra on R d , its transition probability is defined as

P t (A, x) := P (X(t) ∈ A|X(0) = x) , (4)
where A ∈ B(R d ). This corresponds to the probability that the process reaches a Borel set A ⊂ R d at time t, provided that it started in x ∈ R d at time 0 < t.

Noise structure: ellipticity and hypoellipticity

Depending on the noise structure, two classes of models are obtained. The first class is called elliptic and corresponds to SDEs with a non-degenerate diffusion matrix, i.e., ΣΣ is of full rank.

In particular, we consider the case d = m and a diagonal matrix Σ = diag[σ 1 , . . . , σ d ] with entries σ j > 0 for j = 1, . . . , d.

The second class corresponds to SDEs with degenerate diffusion matrix, as it naturally occurs in many application models. Following the notion in [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF], we consider m = d -1 and Σ given by

Σ :=   0 d-1 Γ   , (5)
where

Γ = diag[σ 1 , . . . , σ d-1 ] ∈ R (d-1)×(d-1
) is a diagonal matrix with entries σ j > 0 for j = 1, . . . , d -1. The first component of the solution (X(t)) t∈[0,T ] is called smooth, since it is not directly affected by the noise. The remaining d -1 components are called rough, because the noise acts directly on them. In this scenario, SDE (1) is often hypoelliptic. It means that the transition probability (4) admits a smooth density, even though ΣΣ is not of full rank. This is the case when the SDE satisfies the weak Hörmander condition, based on the concept of Lie-brackets [START_REF] Nualart | The Malliavin Calculus and Related Topics[END_REF]. In [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF], it was shown that a necessary and sufficient condition for the SDE to meet the weak Hörmander condition is that at least one of the rough coordinates of the process (X(t)) t∈[0,T ] appears in the first component F 1 (X(t)) of the drift, that is

∀x ∈ R d , ∂ r F 1 (x), σ j = 0 for at least one j = 1, . . . , d -1, (6)
where σ j denotes the j-th column vector of Γ and ∂ r F 1 (x) := (∂ x2 F 1 (x), . . . , ∂ x d F 1 (x)) is the vector of partial derivatives of the first entry of the drift with respect to the rough components. This setting can be extended to multiple smooth coordinates by requiring that at least one of the rough coordinates enters those components of the drift, where noise does not directly act upon.

Lyapunov structure: geometric ergodicity

Here, a particular interest lies in SDEs of type [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF] where the drift F (X(t)) satisfies the following dissipativity condition

(F (x), x) ≤ α -δ x 2 , ∀ x ∈ R d , (7)
where α, δ > 0. Condition [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF] ensures that the function L :

R d → [1, ∞) defined by L(x) := 1 + x 2
is a Lyapunov function for [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF], see [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. That is L(x) → ∞ as x → ∞, and there exist constants ρ, η > 0 such that

L{L(x)} ≤ -ρL(x) + η, ( 8
)
where L is the generator of the SDE given by

L{g(x)} = d l=1 F l (x) ∂g ∂x l (x) + 1 2 d l,j=1 ΣΣ lj ∂ 2 g ∂x l ∂x j (x),
for sufficiently smooth functions g : R d → R. The existence of a Lyapunov function satisfying [START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF] is the key to establish the geometric ergodicity of the solution of (1). This property means that the distribution of the Markov process (X(t)) t∈[0,T ] converges exponentially fast to a unique invariant distribution π, satisfying

π(A) = R d P t (A, x)π(dx), ∀ A ∈ B(R d ), t ∈ [0, T ].
In particular, if SDE (1) is elliptic, the existence of a Lyapunov function meeting Condition [START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF] suffices to establish the geometric ergodicity of (X(t)) t∈[0,T ] . If SDE (1) is not elliptic, the process is geometrically ergodic, if, in addition to fulfilling Condition [START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF], it is hypoelliptic and satisfies the irreducibility condition P t (A, x) > 0 for all open sets A ∈ B(R d ) and x ∈ R d . The reader is referred to [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF] and the references therein for further details.

Splitting method

Consider a discretised time interval [0, T ] with equidistant time steps ∆ = t i -t i-1 ∈ (0, ∆ 0 ], ∆ 0 ∈ (0, 1), i = 1, . . . , n, where t 0 = 0 and t n = T . Throughout, we denote by ( X(t i )) i=0,...,n a numerical solution of SDE (1), approximating the process (X(t)) t∈[0,T ] at t i , where X(0) := X 0 . A numerical splitting solution is obtained based on the following three steps [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF][START_REF] Mclachlan | Splitting methods[END_REF]:

(i) Split the equation of interest into exactly solvable subequations, which may consist of deterministic and/or stochastic dynamical systems;

(ii) Derive the exact solutions of these subequations;

(iii) Compose the derived solutions in a proper way.

In this section, the following splitting strategy is proposed for SDEs of type [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF].

Step (i): Choice of the subequations To make use of the treatable underlying stochastic linear dynamics, we propose to split Equation (1) into the following two subequations dX [1] (t) = AX [1] (t)dt + ΣdW (t), X [1] (0) = X

[1] 0 , t ∈ [0, T ], (9)
dX [2] (t) = N (X [2] (t))dt, X [2] (0

) = X [2] 0 , t ∈ [0, T ]. ( 10
)
This splitting strategy is an extension of the method presented in [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF], where the authors consider a globally Lipschitz Hamiltonian type equation with uniformly bounded non-linear terms.

Our method considers a more general class of coefficients N (X(t)), including functions which are allowed to grow polynomially at infinity according to Assumption 1.

Step (ii): Exact solution of the subequations In the following, we discuss the subequations ( 9) and [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF], and denote by ϕ t (X 0 ), k = 1, 2, their exact solutions (flows) at time t and starting from X 0 . The first subequation ( 9) is a linear SDE. It can be solved exactly, even when the dimension d is large and independent of whether the equation has an elliptic or hypoelliptic noise structure [START_REF] Arnold | Stochastic Differential Equations: Theory and Applications[END_REF][START_REF] Mao | Stochastic Differential Equations and Applications[END_REF]. In particular, the exact solution of ( 9) is given by

X [1] (t) = e At X [1] 0 + t 0 e A(t-s) Σ dW (s). ( 11
)
The Itô integral in ( 11) is normally distributed with mean 0 d . Moreover, using Itô's isometry and the fact that the components of the Wiener process are independent, its d × d-dimensional covariance matrix is given by

C(t) = t 0
e A(t-s) ΣΣ (e A(t-s) ) ds.

Hence, paths of (9) can be simulated exactly at the discrete time points t i . In particular, ϕ

∆ (X [1] (t i-1 )) := X [1] (t i ) = e A∆ X [1] (t i-1 ) + ξ i-1 , i = 1, . . . , n,

where the ξ i-1 are independent and identically distributed d-dimensional Gaussian vectors with mean 0 d and covariance matrix C(∆) given by [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF]. Due to (A1) of Assumption 1, the global solution of the second subequation (10) exists, i.e., it is defined on the entire interval [0, T ] such that it does not blow up to infinity in finite time [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF]. At the discrete time points t i , we have

ϕ [2]
∆ (X [2] (t i-1 )) := X [2] (t i ) = f (X [2]

(t i-1 ); ∆), i = 1, . . . , n, (14)
where f : R d → R d denotes the exact solution of Equation [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF].

Remark 1. To establish the boundedness of the moments (Lemma 2), the Lyapunov condition (Theorem 4) and the asymptotic second moment bound (Corollary 1 and 2), we exploit properties of the exact solution f of Equation [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF]. These will be illustrated on a cubic model problem and the FHN model (2) in Section 5 and Section 6, respectively. However, note that some of the results presented in the following are formulated with conditions not involving f directly (see Lemma 1 on the mean-square consistency and Theorem 3 on the 1-step hypoellipticity). Therefore, these results would also hold when a numerical method to approximate the solution of Equation (10) is used. We refer to [START_REF] Hairer | Solving Ordinary Differential Equations: Nonstiff Problems[END_REF][START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF] for an exhaustive discussion of numerical methods for locally Lipschitz ODEs.

Step (iii): Composition of the exact solutions To finally obtain a numerical solution of SDE (1), the exact solutions ( 13) and ( 14) of the subequations ( 9) and [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] are composed in every iteration step. In particular, we investigate the following explicit method

X LT (t i ) = ϕ [1] ∆ • ϕ [2] ∆ X LT (t i-1 ) = e A∆ f ( X LT (t i-1 ); ∆) + ξ i-1 , (15)
which is based on the Lie-Trotter (LT) composition approach [START_REF] Trotter | On the product of semi-groups of operators[END_REF]. Note that the matrix exponential e A∆ and the covariance matrix C(∆) required in [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] have to be precomputed only once for a given time step ∆, and the normal random variables ξ i-1 , i = 1, . . . , n, can be obtained using a Cholesky decomposition of the covariance matrix C(∆).

Mean-square convergence

In this section, mean-square convergence of order 1 is proved for the constructed splitting method. It has been observed in the globally Lipschitz case that splitting methods have the same convergence order as the standard Euler-Maruyama method, i.e., order 1 in the case of additive noise, see, e.g., [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Alamo | A technique for studying strong and weak local errors of splitting stochastic integrators[END_REF][START_REF] Milstein | Quasi-symplectic methods for Langevin-type equations[END_REF]. We extend this result to the one-sided Lipschitz case.

Throughout this section, K denotes a generic constant, which may depend on T , p and ∆ 0 , but is independent of ∆ and i.

Required background

To establish mean-square convergence, we rely on Theorem 2.1 of [START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF], which provides an extension of Milstein's fundamental theorem on the mean-square order of convergence for globally Lipschitz coefficients [START_REF] Milstein | A theorem on the order of convergence of mean-square approximations of solutions of systems of stochastic differential equations[END_REF] (see also Theorem 1.1 in [START_REF] Milstein | Stochastic Numerics for Mathematical Physics. Scientific computation[END_REF]) to the considered setting specified in Assumption 1.

To facilitate the illustration of our results, we recall this statement in Theorem 1 below, after defining the required ingredients of mean-square consistency and boundedness of moments.

Let X ti-1,x (t i ) denote the true solution at time t i starting from x at time t i-1 , i.e., X(t i-1 ) = x, and X ti-1,x (t i ) the one-step approximation used to construct a numerical solution X(t i ). In particular, the one-step approximation of the numerical method discussed in the previous section is defined by [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF], where X LT (t i-1 ) is replaced by x.

Definition 1. The one-step approximation X ti-1,x (t i ) of a numerical solution X(t i ) of SDE (1) is mean-square consistent of order q 2 -1/2, if for some p ≥ 1, there exist α ≥ 1, ∆ 0 > 0 and K > 0 such that for arbitrary t i , i = 1, . . . , n, x ∈ R d , and all ∆ ∈ (0, ∆ 0 ], it holds that

E X ti-1,x (t i ) -X ti-1,x (t i ) ≤ K 1 + x 2α 1/2 ∆ q1 , E X ti-1,x (t i ) -X ti-1,x (t i ) 2p 1/(2p) ≤ K 1 + x 2αp 1/(2p) ∆ q2 , with q 2 ≥ 1/2 and q 1 ≥ q 2 + 1/2.
Besides mean-square consistency, the boundedness of the moments of the numerical solution has to be proved. In the globally Lipschitz case, this is guaranteed by the linear growth bounds of the coefficients.

Definition 2. A numerical solution X(t i ) of SDE (1)
has bounded moments, if for any p ≥ 1, there exist ∆ 0 > 0 and K > 0 such that for all ∆ ∈ (0, ∆ 0 ] and i = 0, . . . , n, it holds that

E X(t i ) 2p ≤ K 1 + E X 0 2p .
Based on the above defined ingredients, the following theorem guarantees mean-square convergence.

Theorem 1 (Theorem 2.1 in Tretyakov and Zhang (2013) [START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF]). Let X(t i ) denote a numerical solution of SDE (1) at time t i starting at X 0 , constructed using the one-step approximation X ti-1,x (t i ).

Further, let Assumption 1 hold. If

(i) The one-step approximation X ti-1,x (t i ) is mean-square consistent of order q 2 -1/2 in the sense of Definition 1.

(ii) The numerical method X(t i ) has bounded moments in the sense of Definition 2.

Then the numerical method X(t i ) is mean-square convergent of order q 2 -1/2, i.e., for any n and i = 0, . . . , n, the following inequality holds:

E X(t i ) -X(t i ) 2p 1/(2p) ≤ K 1 + E X 0 2pc 1/(2p) ∆ q2-1/2 ,
where K > 0 and c ≥ 1.

Mean-square convergence of the splitting method

In the following, we prove the required Conditions (i) and (ii) of Theorem 1 for the constructed splitting method. Condition (i) can be proved in a similar fashion as Lemma 2.1 in [START_REF] Milstein | Quasi-symplectic methods for Langevin-type equations[END_REF] (globally Lipschitz case) and Lemma 3.2 in [START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF][START_REF] Zhang | Order-preserving strong schemes for SDEs with locally Lipschitz coefficients[END_REF] (locally Lipschitz case). These proofs rely on the mean-square consistency of the Euler-Maruyama method, which is given by

X EM (t i ) = X EM (t i-1 ) + F ( X EM (t i-1 ))∆ + Σ √ ∆ψ i-1 , (16)
where the ψ i-1 ∼ N (0 m , I m ), i = 1, . . . , n, are independent and identically distributed mdimensional standard Gaussian vectors [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF][START_REF] Milstein | Stochastic Numerics for Mathematical Physics. Scientific computation[END_REF].

Lemma 1 (Mean-square consistency). Let X LT ti-1,x (t i ) be the one-step approximation of the splitting method defined through [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] and let Assumption 1 hold. Further, assume that the drift F (x) has continuous first and second order derivatives in x which satisfy a polynomial growth condition of the form (A2). Then X LT ti-1,x (t i ) is mean-square consistent of order 1 in the sense of Definition 1. In particular, q 1 = 2 and q 2 = 3/2.

Proof. Consider first the one-step approximation

X EM ti-1,x (t i ) = x + F (x)∆ + Σ √ ∆ψ i-1 (17)
of the Euler-Maruyama method ( 16) applied to SDE (1). Since we consider the case of additive noise and (A2) holds for the drift F (x) = Ax + N (x) by assumption of Lemma 1, it follows from the proof of Lemma 3.2 in [START_REF] Zhang | Order-preserving strong schemes for SDEs with locally Lipschitz coefficients[END_REF] that (17) satisfies Definition 1 with q 1 = 2, q 2 = 3/2. Thus, it suffices to compare the splitting and Euler-Maruyama methods.

Since the drift of the stochastic subequation (9) of the splitting framework grows linearly, using again Lemma 3.2 in [START_REF] Zhang | Order-preserving strong schemes for SDEs with locally Lipschitz coefficients[END_REF] (see also Section 1.1.5 in [START_REF] Milstein | Stochastic Numerics for Mathematical Physics. Scientific computation[END_REF]), its solution can be expressed as

ϕ [1] ∆ (x) = x + Ax∆ + Σ √ ∆ψ i-1 + r s (x, ∆), (18)
where r s (x, ∆) satisfies the inequalities of Definition 1 for α = 1, q 1 = 2 and q 2 = 3/2. In particular, we have that

E [r s (x, ∆)] ≤ K 1 + x 2 1/2 ∆ 2 , E r s (x, ∆) 2p 1/(2p) ≤ K 1 + x 2p 1/(2p) ∆ 3/2 . (19)
Similarly, since N (x) satisfies (A2), the solution of the deterministic subequation [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] of the splitting framework can be expressed as

ϕ [2] ∆ (x) = x + N (x)∆ + r d (x, ∆), (20)
where

r d (x, ∆) ≤ K 1 + x 2α 1/2 ∆ 2 , (21)
for some α ≥ 1. The one-step approximation of the Lie-Trotter splitting method is then obtained by composing the above expressions ( 18) and ( 20), yielding

X LT ti-1,x (t i ) = (ϕ [1] ∆ • ϕ [2] ∆ )(x) = x + N (x)∆ + r d (x, ∆) + Ax∆ + AN (x)∆ 2 + Ar d (x, ∆)∆ +Σ √ ∆ψ i-1 + r s x + N (x) + r d (x, ∆), ∆ .
Thus, the difference between the splitting and Euler-Maruyama methods becomes

r LT (x, ∆) := X LT ti-1,x (t i ) -X EM ti-1,x (t i ) = r d (x, ∆) + AN (x)∆ 2 + Ar d (x, ∆)∆ + r s x + N (x) + r d (x, ∆), ∆ .
Using (A2) and the inequalities ( 19) and ( 21), it follows that r LT (x, ∆) also satisfies the inequalities in Definition 1 for q 1 = 2, q 2 = 3/2 and some α ≥ 1. This concludes the proof. Now, we establish the boundedness of the moments of the splitting method. Intuitively, this is guaranteed by the use of the global exact solution of the locally Lipschitz ODE [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF], which is defined on the entire interval [0, T ] without any explosion occuring in finite time. Thus, the iterative composition of this function with the solution of the linear SDE via the Lie-Trotter method (15) does not cause an explosion of the moments in finite time either. The formal proof of this result, provided in Lemma 2, is done in the spirit of the proof of Proposition 3 in [START_REF] Bréhier | Analysis of some splitting schemes for the stochastic Allen-Cahn equation[END_REF].

Lemma 2 (Boundedness of moments). Let X LT (t i ) be the splitting method defined through [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] and let Assumption 1 hold. Then X LT (t i ) is mean-square bounded in the sense of Definition 2.

Proof. Consider the linear SDE

dZ(t) = AZ(t)dt + ΣdW (t), Z(0) = Z 0 = 0 d .

Its exact solution is given by

Z(t) = t 0 e A(t-s) ΣdW (s),
where Z(t) is normally distributed with mean vector 0 d and covariance matrix C(t) as defined in [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF]. Consequently, the moments of Z(t) are bounded, i.e., for any p ≥ 1 there exists K Z (T, p) > 0 such that

E sup 0≤t≤T Z(t) 2p ≤ K Z (T, p). (22)
Now, define the process R(t i ) := X LT (t i ) -Z(t i ). It suffices to prove the boundedness of the moments of R(t i ). Note that in a discretised regime we have that

Z(t i ) = e A∆ Z(t i-1 ) + ξ i-1 . Thus, R(t i ) = e A∆ f ( X LT (t i-1 ); ∆) -Z(t i-1 ) = e A∆ f (R(t i-1 ) + Z(t i-1 ); ∆) -f (Z(t i-1 ); ∆) + f (Z(t i-1 ); ∆) -Z(t i-1 ) .
Using that e A∆ x ≤ e A∆ x ≤ e µ(A)∆ x for all x ∈ R d , we obtain

R(t i ) ≤ e µ(A)∆ f (R(t i-1 ) + Z(t i-1 ); ∆) -f (Z(t i-1 ); ∆) +e µ(A)∆ f (Z(t i-1 ); ∆) -Z(t i-1 ) .
Since the function N : R d → R d satisfies (A1), using the continuous Gronwall Lemma, the function f : R d → R d fulfils the following global Lipschitz condition

f (x; ∆) -f (y; ∆) ≤ e c1∆ x -y , ∀ x, y ∈ R d ,
where the constant c 1 is the same as in Assumption (A1), see, e.g., Theorem 1.2.17 in [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF]. Moreover, using the Taylor series expansion (20) f (x; ∆) = x + N (x)∆ + r d (x, ∆), and applying (A2) and ( 21), we obtain

f (x; ∆) -x = ∆N (x) + r d (x, ∆) ≤ ∆ N (x) + r d (x, ∆) ≤ c(1 + x ĉ)∆,
where c and ĉ are positive constants. Thus, defining c := max{|µ(A)|, c 1 } > 0, we get that

R(t i ) ≤ e c∆ R(t i-1 ) + e c∆ ∆c 1 + Z(t i-1 )
ĉ .

Now, we can perform back iteration, obtaining

R(t i ) ≤ e cti R 0 + c∆ i k=1 e ck∆ 1 + Z(t i-k ) ĉ ≤ e cT X 0 + c 1 + sup 0≤l≤i-1 Z(t l ) ĉ ∆ i k=1 e ck∆ ,
where we used that R 0 = X 0 , since

Z 0 = 0 d . Using that ∆ i k=1 e ck∆ = (e cti -1) ∆e c∆ e c∆ -1 ≤ (e cT -1) ∆ 0 e c∆0 e c∆0 -1 , ∀ ∆ ∈ (0, ∆ 0 ],
we get that

R(t i ) ≤ e cT X 0 + c(e cT -1) ∆ 0 e c∆0 e c∆0 -1 1 + sup 0≤l≤i-1 Z(t l ) ĉ .
Thus, there exists a constant K(T, ∆ 0 ) > 0 such that

R(t i ) ≤ K(T, ∆ 0 ) 1 + X 0 + sup 0≤l≤i-1 Z(t l ) ĉ .
Considering the 2p-th moments and using [START_REF] Hairer | Geometric Numerical Integration[END_REF] concludes the proof.

Based on the above results, we establish the mean-square convergence of the splitting method in the following theorem.

Theorem 2 (Mean-square convergence). Let X LT (t i ) be the splitting method defined through [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] and let the assumptions of Theorem 1, Lemma 1 and Lemma 2 hold. Then X LT (t i ) is mean-square convergent of order 1.

Proof. The result is a direct consequence of Theorem 1, Lemma 1 and Lemma 2.

Note that, in contrast to ODEs [START_REF] Hairer | Geometric Numerical Integration[END_REF], the mean-square convergence order of splitting methods for SDEs cannot be increased by using compositions based on fractional steps. Indeed, to achieve this in the stochastic scenario, higher-order stochastic integrals would be required [START_REF] Milstein | Quasi-symplectic methods for Langevin-type equations[END_REF]. Thus, the splitting method

X S (t i ) = ϕ [2] ∆/2 • ϕ [1] ∆ • ϕ [2] ∆/2 X S (t i-1 ) = f e A∆ f X S (t i-1 ); ∆/2 + ξ i-1 ; ∆/2 , ( 23
)
which is based on the Strang (S) composition approach [START_REF] Strang | On the construction and comparison of difference schemes[END_REF], is expected to also have mean-square order 1. Nevertheless, it has been observed that Strang methods may perform better than Lie-Trotter methods in numerical experiments, possibly due to the symmetry of this composition method, see, e.g., [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Chen | Structure-preserving numerical integrators for Hodgkin-Huxley-type systems[END_REF][START_REF] Chevallier | Diffusion approximation of multi-class Hawkes processes: Theoretical and numerical analysis[END_REF][START_REF] Tubikanec | Qualitative properties of different numerical methods for the inhomogeneous geometric Brownian motion[END_REF]. Thus, the Strang method ( 23) is also considered in the numerical experiments reported in Section 7.

Structure preservation

The mean-square convergence discussed in the previous section is a limit result for the time discretisation step ∆ going to zero over a finite interval. This result does not carry any information about the quality of the numerical method under the use of strictly positive time steps ∆, as always required when implementing any numerical method. In the following, we discuss the preservation of important structural properties, focusing on hypoellipticity and ergodicity.

Preservation of noise structure and 1-step hypoellipticity

To obtain a discrete analogue of the transition probability (4) introduced in Section 2, we define the k-step transition probability of a numerical solution X(t i ) of SDE (1) as follows

P t k (A, x) := P( X(t k ) ∈ A| X(0) = x), (24)
where A ∈ B(R d ) and x ∈ R d . Now, assume that SDE (1) is hypoelliptic, i.e., its transition probability (4) has a smooth density even though ΣΣ is not of full rank, see Section 2.1. We introduce a discrete version of this property in the subsequent definition.

Definition 3 (k-step hypoellipticity). Let X(t i ) be a numerical solution of (1) and k ∈ N be the smallest k such that its transition probability (24) has a smooth density. Then,

X(t i ) is called k-step hypoelliptic.
This means that the numerical method propagates the noise into the smooth component after k iteration steps. The preservation of this property is not an issue when using the numerical method to simulate paths of the SDE over a large enough time horizon, as standard methods usually satisfy it for some k. For example, the Euler-Maruyama method has been observed to be 2-step hypoelliptic, see, e.g., Corollary 7.4 in [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. However, the case k = 1, where we also use the notation

P ∆ (A, x) := P( X(t i ) ∈ A| X(t i-1 ) = x), (25)
is of crucial relevance when using the numerical method within statistical applications. For example, in the field of likelihood-based parameter estimation, explicit numerical methods are used to approximate transition densities [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF][START_REF] Melnykova | Parametric inference for hypoelliptic ergodic diffusions with full observations[END_REF][START_REF] Pokern | Parameter estimation for partially observed hypoelliptic diffusions[END_REF]. In this regard, a particular interest lies in situations where [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] corresponds to a non-degenerate multivariate normal distribution, i.e., X(t i ) given X(t i-1 ) is normally distributed with a covariance matrix that reflects the propagation of the noise to the smooth components. This is not the case for the Euler-Maruyama method [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF], as it yields a degenerate multivariate normal transition distribution with conditional covariance matrix given by Cov(

X EM (t i )| X EM (t i-1 )) = ∆ΣΣ .
Note that the same degenerate covariance matrix is obtained by tamed/truncated variants of the Euler-Maruyama method (see Section 7.1). In contrast, the conditional covariance matrix of the Lie-Trotter splitting (15) coincides with C(∆), as defined in [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF]. Thus, if the stochastic linear subequation (9) of the splitting framework is hypoelliptic, the proposed splitting method yields a non-degenerate multivariate normal transition distribution.

Assumption 2. The matrix A is such that SDE (9) is hypoelliptic.

Theorem 3 (1-step hypoellipticity). Let X LT (t i ) be the splitting method defined through [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] and let Assumption 2 hold. Then, X LT (t i ) is 1-step hypoelliptic according to Definition 3. Moreover, X LT (t i ) given X LT (t i-1 ) admits a non-degenerate normal distribution with mean vector and covariance matrix given by

E X LT (t i )| X LT (t i-1 ) = e A∆ f ( X LT (t i-1 ); ∆), Cov( X LT (t i )| X LT (t i-1 )) = C(∆),
respectively, where C(∆) is defined in [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF].

Proof. The fact that X LT (t i ) given X LT (t i-1 ) is normally distributed with the corresponding mean vector and covariance matrix is an immediate consequence of formula [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF], recalling that the ξ i are Gaussian random vectors with null mean and covariance matrix C(∆). Moreover, the linear SDE ( 9) is hypoelliptic by assumption. Thus, its solution (X [1] (t)) t∈[0,T ] has conditional covariance matrix C(t) = Cov(X [1] (t)|X [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF] 0 ) [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF] which is of full rank. Since the covariance matrix of the Lie-Trotter splitting equals C(∆), this method is 1-step hypoelliptic according to Definition 3, and thus the normal distribution is non-degenerate.

Remark 2. Note that the 1-step hypoellipticity of the Lie-Trotter splitting [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] and the fact that this method yields a non-degenerate normal transition distribution according to Theorem 3 hold without requiring any conditions on the nonlinearity of SDE (1). Moreover, even though the transition distribution of the Strang splitting (23) is not explicitly available in general, this numerical method is expected to be 1-step hypoelliptic too, since it also benefits from the covariance matrix C(∆) (12).

Preservation of Lyapunov structure and geometric ergodicity

We now assume that SDE ( 1) is geometrically ergodic. The main task to establish the geometric ergodicity of a numerical solution of ( 1) is to prove a discrete analogue of the Lyapunov condition [START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF] introduced in Section 2.2.

Definition 4 (Discrete Lyapunov condition). Let L be a Lyapunov function for SDE [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF]. A numerical solution X(t i ) of (1) satisfies the discrete Lyapunov condition if there exist ρ ∈ (0, 1) and η ≥ 0 such that

E L( X(t i ))| X(t i-1 ) ≤ ρL( X(t i-1 )) + η, ∀ i ∈ N.
Analogously to the continuous case, this condition implies geometric ergodicity of the numerical method if SDE (1) is elliptic. If the equation is not elliptic, in addition to the discrete Lyapunov condition, k-step hypoellipticity and a discrete irreducibility condition are required. For further details, the reader is referred to [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF].

Euler-Maruyama type methods do not preserve this property, especially when the drift of SDE (1) is only locally Lipschitz continuous. In particular, the problem does not lie in the preservation of hypoellipticity and irreducibility, but in preserving the Lyapunov structure [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]. Consider, for example, the cubic one-dimensional SDE

dX(t) = -X 3 (t)dt + σdW (t), X(0) = X 0 . ( 26
) Since F (x)x = -x 4 ≤ 1 -x 4 ≤ 2 -2x 2
, this SDE satisfies the dissipativity condition [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF]. Thus, L(x) = 1 + x 2 is a Lyapunov function satisfying [START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF] and the process (X(t)) t∈[0,T ] is geometrically ergodic. However, it is shown in Lemma 6.3 of [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF] that, if E[X 2 0 ] ≥ 2/∆, the second moment of the Euler-Maruyama method goes to infinity as the time t i grows, since

E X EM (t i ) 2 ≥ E[X 2 0 ] + t i .
Thus, for any fixed time step ∆ > 0 (even when it is chosen to be arbitrarily small), one can find a starting value X 0 such that the Euler-Maruyama method does not converge to a unique invariant distribution. This also means that for any ∆ > 0 and X 0 , there is a positive probability of blow-up, as discussed in [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF].

In contrast, splitting methods may preserve the Lyapunov structure [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF][START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF]. This is also proved for the proposed splitting method and the Lyapunov function L(x) = 1 + x 2 , under an additional Assumption on the function f and the matrix A, respectively. Assumption 3. There exists a constant c 3 ≥ 0 such that for any x ∈ R d , it holds that

f (x; ∆) 2 ≤ x 2 + c 3 ∆, ∀∆ ∈ (0, ∆ 0 ].
Assumption 4. The matrix A is such that e A∆ < 1 for all ∆ ∈ (0, ∆ 0 ].

Theorem 4 (Discrete Lyapunov condition). Let X LT (t i ) be the splitting method defined through [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF], and let Assumptions 3 and 4 hold. Then X LT (t i ) satisfies the discrete Lyapunov condition of Definition 4 with Lyapunov function L(x) = 1 + x 2 .

Proof. We have that

X LT (t i ) 2 = e A∆ f ( X LT (t i-1 ); ∆) + ξ i-1 2 = f ( X LT (t i-1 ); ∆) (e A∆ ) (e A∆ )f ( X LT (t i-1 ); ∆) +f ( X LT (t i-1 ); ∆) (e A∆ ) ξ i-1 + ξ i-1 e A∆ f ( X LT (t i-1 ); ∆) + ξ i-1 ξ i-1 .
Denoting the diagonal entries of the covariance matrix C(∆) (12) by c jj (∆), taking the expectation, using the fact that X LT (t i-1 ) and ξ i-

1 are independent, that E[ξ i-1 ] = 0 d , that E[ξ i-1 ] = 0 d and that C(∆) := d j=1 c jj (∆) = E[ξ i-1 ξ i-1 ],
we get that

E X LT (t i ) 2 = E e A∆ f ( X LT (t i-1 ); ∆) 2 + C(∆). (27)
Considering

L( X LT (t i )) = 1 + X LT (t i ) 2 ,
and using [START_REF] Hutzenthaler | Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF] and Assumption 3, we get that

E L( X LT (t i ))| X LT (t i-1 ) = 1 + e A∆ f ( X LT (t i-1 ); ∆) 2 + C(∆) ≤ 1 + e A∆ 2 X LT (t i-1 ) 2 + c 3 ∆ + C(∆) + e A∆ 2 = e A∆ 2 L( X LT (t i-1 )) + 1 + e A∆ 2 c 3 ∆ + C(∆),
where we added e A∆ 2 in the inequality. Thus, applying Assumption 4, the discrete Lyapunov condition of Definition 4 is satisfied for ρ = e A∆ 2 < 1 and η = 1 + e A∆ 2 c 3 ∆ + C(∆) > 0, which proves the result.

In the following corollary of Theorem 4, we show that the second moment of the splitting method is asymptotically bounded by a constant which is independent of T , ∆ and i. In particular, this bound is reached exponentially fast, independently of the choice of X 0 , in agreement with the geometric ergodicity of the splitting method. This result also requires Assumption 3 and an assumption related to the matrix A.

Assumption 5. The matrix A is such that the logarithmic norm µ(A) < 0.

Note that Assumption 5 implies Assumption 4, since e A∆ ≤ e µ(A)∆ [START_REF] Ström | On logarithmic norms[END_REF]. However, the converse is not true in general. Assumption 5 is, e.g., satisfied for normal matrices, where all eigenvalues have strictly negative real part [START_REF] Söderlind | The logarithmic norm. History and modern theory[END_REF]. Matrices contained in this class are, e.g., diagonal ones with strictly negative diagonal entries.

Corollary 1 (Asymptotic second moment bound). Let X LT (t i ) be the splitting method defined through [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF], and let Assumptions 3 and 5 hold. Then, there exists a constant K LT ∞ > 0, which is independent of T , ∆ and i, such that

lim ti→∞ E X LT (t i ) 2 ≤ K LT ∞ .
Proof. Recalling ( 27) from the proof of Theorem 4, and using Assumption 3 and the logarithmic norm, we further obtain

E X LT (t i ) 2 ≤ e 2µ(A)∆ E X LT (t i-1 ) 2 + c 3 ∆ + C(∆).
Now, we can perform back iteration, yielding

E X LT (t i ) 2 ≤ e 2µ(A)ti E X 0 2 + c 3 ∆ i k=1 e 2µ(A)k∆ + C(∆) i-1 k=0 e 2µ(A)k∆ .
Using that

i k=1 e 2µ(A)k∆ = 1 -e 2µ(A)ti e 2µ(A)∆ 1 -e 2µ(A)∆ , i-1 k=0 e 2µ(A)k∆ = 1 -e 2µ(A)ti 1 1 -e 2µ(A)∆ ,
we obtain

E X LT (t i ) 2 ≤ e 2µ(A)ti E X 0 2 + 1 -e 2µ(A)ti c 3 ∆e 2µ(A)∆ 1 -e 2µ(A)∆ + C(∆) 1 -e 2µ(A)∆ . ( 28
)
Applying Assumption 5, yields

lim ti→∞ E X LT (t i ) 2 ≤ c 3 ∆e 2µ(A)∆ 1 -e 2µ(A)∆ + C(∆)
1 -e 2µ(A)∆ . Now, we have that

∆e 2µ(A)∆ 1 -e 2µ(A)∆ ≤ - 1 2µ(A) , ∀ ∆ > 0 and ∆ 1 -e 2µ(A)∆ ≤ ∆ 0 1 -e 2µ(A)∆0 , ∀ ∆ ∈ (0, ∆ 0 ]. (29)
Moreover, recalling that e A∆ = I d + ∆A + O(∆ 2 ), it follows from (12) that C(∆) = O(∆). This implies the result. Remark 3. Theorem 4 and Corollary 1 can be proved similarly for the Strang splitting method [START_REF] Hairer | Solving Ordinary Differential Equations: Nonstiff Problems[END_REF].

The one-dimensional case Consider the case d = 1, Σ = σ > 0 and A = -a < 0. In this case, the solution of the linear SDE (9) corresponds to the Ornstein-Uhlenbeck process X [1] (t) = e -at X [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF] 0 + σ t 0 e -a(t-s) dW (s), [START_REF] Hutzenthaler | Ergodic behavior of locally regulated branching populations[END_REF] with variance [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF] given by

C(t) = σ 2 2a (1 -e -2at ). (31)
Thus, due to the specific form of (31), the previously derived bound can be expressed in closedform for any time t i . In particular, the following (asymptotic) bound for the second moment of the splitting method ( 15) is obtained.

Corollary 2 (Closed-form (asymptotic) second moment bound). Let d = 1, Σ = σ > 0 and A = -a < 0. Further, let X LT (t i ) be the splitting method defined through [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF], and let Assumption 3 be satisfied. Then, it holds that

E ( X LT (t i )) 2 ≤ K LT (t i , X 0 ) := e -2ati E X 2 0 + (1 -e -2ati ) c 3 2a + σ 2 2a , lim ti→∞ E ( X LT (t i )) 2 ≤ K LT ∞ := c 3 2a + σ 2 2a .
Proof. Using [START_REF] Kamke | Differentialgleichungen Lösungsmethoden und Lösungen: Gewöhnliche Differentialgleichungen[END_REF] and noting that C(∆) = C(∆) and that µ(A) = -a < 0, the result is a direct consequence of Corollary 1.

Note that, for t i = 0, the bound K LT (0, X 0 ) in Corollary 2 coincides with E[X 2 0 ]. Moreover, since a > 0, the distribution of (30) converges to a unique limit

X [1] (t) D ---→ t→∞ N 0, σ 2 2a . (32)
Intuitively, this fact, combined with Assumption 3, guarantees the geometric ergodicity of the splitting method obtained via Theorem 4, and thus the existence of the asymptotic bound for the second moment reported in Corollary 2.

Cubic model problem For an illustration of the derived bound, consider again SDE [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF]. We propose to rewrite this equation as

dX(t) = -X(t) + X(t) -X 3 (t) dt + σdW (t),
and to choose

A = -1 < 0, N (X(t)) = X(t) -X 3 (t). ( 33
)
The exact solution of the resulting linear SDE ( 9) is then given by (30) for a = 1, and that of ODE ( 10) is given by

X [2] (t) = f (X [2] 0 ; t) = X [2] 0 e -2t + (X [2] 0 ) 2 (1 -e -2t ) . ( 34
)
This choice guarantees that all required assumptions are satisfied.

Proposition 1. Let A, N and f be as in [START_REF] Kelly | Adaptive Euler methods for stochastic systems with non-globally Lipschitz coefficients[END_REF] and (34), respectively. Then, Assumptions 1-5 are satisfied.

Proof. The proof is given in Appendix A.

Therefore, the proposed splitting method (15) applied to SDE [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF] is not only mean-square convergent, but also geometrically ergodic. In particular, while even for arbitrarily small ∆ one can find X 0 such that the second moment of the Euler-Maruyama method explodes (see the beginning of Section 5.2), the second moment of the splitting method is bounded by K LT (t i , X 0 ), which converges to the constant K LT ∞ = 1/4 + σ 2 /2 exponentially fast and for any choice of the initial value X 0 , see Corollary 2.

In Figure 1, we illustrate the derived second moment bound K LT (t i , X 0 ) (grey solid line) of Corollary 2 for SDE (26) as a function of the time t i , in comparison with E[X 2 (t i )] (red dashed line). The latter is estimated based on 10 4 paths generated under the Lie-Trotter splitting [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF]. The asymptotic bound K LT ∞ of Corollary 2 is indicated by the blue dotted line.

0 1 2 3 4 5 6 0 1 2 3 4 K LT (t,X 0 ) K ∞ LT E[X 2 (t)]
t Remark 4. For SDE [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF], an immediate choice of the subequations of the splitting framework would also be N (X(t)) = -X 3 (t) and A = 0. For this choice, Assumption 1 related to the locally Lipschitz function N is satisfied, and thus the resulting splitting method (15) is mean-square convergent. Also, Assumptions 2 and 3 are satisfied. However, since e A∆ = 1 and µ(A) = 0, Assumptions 4 and 5 do not hold, asymptotic bounds cannot be derived, and the preservation of ergodicity remains an open question. In particular, in contrast to the proposed approach (see formulas [START_REF] Hutzenthaler | Ergodic behavior of locally regulated branching populations[END_REF] and (32)), the distribution of the solution X [1] (t) = X [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF] 0 + σW (t) of the resulting linear SDE (9) does not converge to a unique limit as t tends to infinity.

Stochastic FitzHugh-Nagumo model

In this section, the proposed splitting strategy is illustrated on the stochastic FHN model, a widely used neuronal model. It is given by the 2-dimensional SDE (2) with solution X(t) := (V (t), U (t)) for t ∈ [0, T ]. This equation has been used to model the firing activity of single neurons [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF]. If the membrane voltage of the neuron is sufficiently high, it releases an action potential, also called spike. The first component (V (t)) t∈[0,T ] describes the membrane voltage of the neuron at time t, while the second component (U (t)) t∈[0,T ] corresponds to a recovery variable modelling the channel kinetics. The parameter > 0 corresponds to the time scale separation of the two components and β ≥ 0 and γ > 0 are position and duration parameters of an excitation, respectively.

Properties of the FHN model

If both noise intensities σ 1 and σ 2 are strictly positive, the model is elliptic. If σ 1 = 0, the diffusion term becomes ΣdW (t) = (0, σ 2 ) dW 2 (t), corresponding to the notation in [START_REF] Arnold | Stochastic Differential Equations: Theory and Applications[END_REF]. In this case, due to the U -component entering the first entry of the drift F (X(t)), the model is hypoelliptic. This is confirmed by the fact that

∂ u F 1 (x)σ 2 = - σ 2 = 0, (35)
guaranteeing Condition [START_REF] Berglund | Mixed-mode oscillations and interspike interval statistics in the stochastic FitzHugh-Nagumo model[END_REF]. We refer to [START_REF] Berglund | Mixed-mode oscillations and interspike interval statistics in the stochastic FitzHugh-Nagumo model[END_REF][START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF][START_REF] Kelly | Adaptive Euler methods for stochastic systems with non-globally Lipschitz coefficients[END_REF][START_REF] Muratov | Noise-induced mixed-mode oscillations in a relaxation oscillator near the onset of a limit cycle[END_REF] and to [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF][START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF] for the consideration of the elliptic and hypoelliptic FHN model, respectively, and to [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] for an investigation of both cases.

Moreover, it has been proved that the FHN model is ergodic, see, e.g., [START_REF] Bonaccorsi | Analysis of the stochastic FitzHugh-Nagumo system[END_REF][START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. Here, we study this property under a restricted parameter space, for which SDE (2) satisfies the dissipativity condition [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF] such that the function L(x) = 1+ x 2 is a Lyapunov function meeting Condition (8).

Proposition 2 (Dissipativity of the FHN model). Let

γ - 1 < 2 min 1 , 1 -τ ,
for some arbitrarily small τ ∈ (0, 1). Then, the drift F of the FHN model (2) satisfies the dissipativity condition [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF].

Proof. We have that

(F (x), x) =   1 (v -v 3 -u) γv -u + β   ,   v u   = 1 (v 2 -v 4 ) + vu(γ - 1 ) -u 2 + βu. Defining c := |γ -1/ |, using 2vu ≤ v 2 + u 2 and v 2 -v 4 ≤ 1 -v 2 , applying Young's inequality βu ≤ u 2 τ /2 + β 2 /2τ
, for some arbitrarily small τ > 0, and setting τ = τ /2, we obtain

(F (x), x) ≤ 1 (1 -v 2 ) + c 2 (v 2 + u 2 ) -u 2 + τ u 2 + β 2 2τ = -v 2 ( 1 - c 2 ) -u 2 (1 -τ - c 2 ) + 1 + β 2 2τ . Since 1 - c 2 > 0 and 1 -τ - c 2 > 0,
by assumption, it follows that (F (x), x) ≤ α -δ x 2 ,

where α = 1/ + β 2 /2τ > 0 and δ = min{1/ -c/2, 1 -τ -c/2} > 0.

Note that the condition on the model parameters in Proposition 2 is satisfied for parameter settings which may be relevant in applications, see Section 7. For example, it is met when γ = 1/ .

Splitting method for the FHN model

The FHN model ( 2) is a semi-linear SDE of type [START_REF] Ableidinger | Splitting integrators for the stochastic Landau-Lifshitz equation[END_REF]. The choice of the matrix A and the function N (X(t)) is not unique. While the locally Lipschitz term -V 3 (t)/ and the constant β have to enter into N (X(t)), the goal is to allocate the remaining terms such that as many of the introduced assumptions as possible are satisfied. For the splitting method to satisfy Assumption 2, and thus to be 1-step hypoelliptic, the term -U (t)/ of the first component of the drift must enter into AX(t). Moreover, shifting the term γV (t) to AX(t) leads to a decoupling of the resulting ODE [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] such that its global solution can be derived exactly and proved to satisfy Assumption 3. Thus, there are four strategies left, depending on whether the remaining terms V (t)/ and -U (t) enter into AX(t) or N (X(t)). The only case where the matrix A meets Assumption 4 (under a restricted parameter space) is when -U (t) appears in AX(t) and V (t)/ in N (X(t)). Similar to the proposed splitting of SDE [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF], the resulting linear SDE ( 9) is then geometrically ergodic. In particular, it corresponds to a version of the well-studied damped stochastic harmonic oscillator whose matrix exponential e At and covariance matrix C(t) have manageable expressions. Therefore, we propose to choose the matrix A and the function N as follows

A =   0 -1 γ -1   , N (X(t)) =   1 V (t) -V 3 (t) β   . ( 36
)
The resulting linear damped stochastic harmonic oscillator [START_REF] Bou-Rabee | Cayley splitting for second-order Langevin stochastic partial differential equations[END_REF] with A as in ( 36) is weakly-, critically-or over-damped, depending on whether

κ := 4γ -1 (37)
is positive, zero or negative, respectively. This terminology, along with the choice of κ, is linked to the roots of the characteristic function of the underlying deterministic equation, see, e.g., Chapter 5 in [START_REF] Weiglhofer | Ordinary Differential Equations and Applications[END_REF]. In particular, the sign of κ determines the shape of the exponential of the matrix A. If κ = 0,

e At = e -t 2   1 + t 2 -t t 4 1 -t 2   . If κ > 0, e At = e -t 2   cos( 1 2 √ κt) + 1 √ κ sin( 1 2 √ κt) -2 √ κ sin( 1 2 √ κt) 2γ √ κ sin( 1 2 √ κt) cos( 1 2 √ κt) -1 √ κ sin( 1 2 √ κt)   .
If κ < 0, the sine and cosine terms of the above expressions can be rearranged using the relations cos 1 2

√ κt = cosh 1 2 √ -κt and 1 √ κ sin 1 2 √ κt = 1 √ -κ sinh 1 2 √ -κt . (38)
Moreover, the covariance matrix C(t) (12) also depends on the sign of κ and is given as follows.

If κ = 0, If κ > 0,

c 11 (t) = e -t 4
c 11 (t) = e -t 2γκ - 4γ 2 (σ 2 1 γ + σ 2 2 1 ) + κe t (σ 2 1 (1 + γ ) + σ 2 2 1 2 ) + σ 2 1 (1 - 3γ ) + σ 2 2 1 2 cos( √ κt) - √ κ(σ 2 1 (1 - γ ) + σ 2 2 1 2 ) sin( √ κt) , c 12 (t) = c 21 (t) = e -t 2κ σ 2 1 κe t - 2 (σ 2 1 γ + σ 2 2 1 ) + σ 2 1 (1 - 2γ ) + 2σ 2 2 1 2 cos( √ κt) -σ 2 1 √ κ sin( √ κt) , c 22 (t) = e -t 2κ (σ 2 2 1 + σ 2 1 γ) cos( √ κt) - 4γ + κe t + (σ 2 2 1 -σ 2 1 γ) √ κ sin( √ κt) .
If κ < 0, the relations [START_REF] Leimkuhler | Molecular Dynamics: With Deterministic and Stochastic Numerical Methods[END_REF] can again be used to rewrite the above expressions accordingly. Note that parameter configurations typically considered in the literature fulfill κ > 0, see, e.g., [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF][START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF][START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF]. This is in agreement with the fact that, under κ > 0, SDE (9) models a weakly damped system which describes oscillatory dynamics. The exact solution of the resulting ODE [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] with N (X(t)) as in (36) reads as

X [2] (t) = f (X [2] 0 ; t) =    V [2] 0 e -2t +(V [2] 0 ) 2 1-e -2t βt + U [2] 0    . ( 39
)
The corresponding Lie-Trotter splitting method for the FHN model ( 2) is then given by [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF], where the matrix exponential e A∆ , the covariance matrix C(∆) and the function f are as reported above.

Properties of the splitting method for the FHN model

In the following proposition, we verify Assumptions 1-4.

Proposition 3. Let A, N and f be as in [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] and (39), respectively. Then the following statements hold.

(i) N satisfies Assumption 1.

(ii) A satisfies Assumption 2.

(iii) If β = 0, then f satisfies Assumption 3.

(iv) If γ = 1/ , then A satisfies Assumption 4.

Proof. The proof is given in Appendix B.

Therefore, the proposed splitting method [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] applied to the FHN model ( 2) is mean-square convergent of order 1, according to Theorem 2. Applying Theorem 3, the method is also 1-step hypoelliptic and yields a non-degenerate Gaussian distribution with covariance matrix C(∆) reported above. This matrix is thus of full rank, even if σ 1 = 0 and independently of the value of κ.

Moreover, for β = 0 and γ = 1/ , L(x) = 1 + x 2 is a Lyapunov function for the FHN model (2) according to Propositon 2 and the method satisfies a discrete Lyapunov condition via Theorem 4. Combined with the 1-step hypoellipticity and a discrete irreducibility condition, which can be proved in the same way as done, e.g, in [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF][START_REF] Chevallier | Diffusion approximation of multi-class Hawkes processes: Theoretical and numerical analysis[END_REF][START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF], the splitting method is geometrically ergodic. Intuitively, the Lyapunov structure of the FHN model is kept by the numerical solution, since the linear SDE (9) determined by the matrix A in ( 36) is geometrically ergodic, implying that the process (X [1] (t)) t∈[0,T ] converges to a unique invariant distribution given by X [1] (t)

D ---→ t→∞ N   0 0   ,   5 2 σ 2 1 + 2 2 σ 2 2 2 σ 2 1 2 σ 2 1 2 8 σ 2 1 + 1 2 σ 2 2   ,
for κ = 0, and

X [1] (t) D ---→ t→∞ N   0 0   ,   2γ (σ 2 1 + γ σ 2 1 + 1 2 σ 2 2 ) 2 σ 2 1 2 σ 2 1 1 2 ( γσ 2 1 + σ 2 2 )   ,
for κ = 0. Since this fact holds without any restrictions of the parameters, it is expected that the splitting method preserves this property for any values of γ > 0 and > 0. This is confirmed by our numerical experiments (see Section 7). Note also that, under γ = 1/ , the logarithmic norm µ(A) = 0. Thus, Assumption 5 is not fulfilled and the asymptotic bound of Corollary 1 cannot be derived.

Remark 5. Another plausible choice of the subequations is

A =   0 -1 γ 0   , N (X(t)) =   1 V (t) -V 3 (t) -U (t) + β   .
For this choice, Assumption 1 related to the locally Lipschitz function N is satisfied, and the splitting method is mean-square convergent. In addition, since the term -U (t)/ still enters into AX(t), Assumption 2 holds, and the method is thus also 1-step hypoelliptic. Moreover, Assumption 3 would also hold under β = 0. However, in this case, the linear SDE (9) corresponds to a version of the simple (undamped) harmonic oscillator which is not ergodic. In particular, its matrix exponential is given by

e At =   cos( √ γt √ ) -1 √ γ sin( √ γt √ ) √ γ sin( √ γt √ ) cos( √ γt √ )   ,
with e A∆ ≥ 1 and e A∆ = 1 for γ = 1/ in particular.

Numerical experiments for the FHN model

We now illustrate the performance of the Lie-Trotter [START_REF] Clairon | Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations[END_REF] and Strang [START_REF] Hairer | Solving Ordinary Differential Equations: Nonstiff Problems[END_REF] splitting methods in comparison with Euler-Maruyama type methods through a variety of numerical experiments carried out on the FHN model [START_REF] Ableidinger | A stochastic version of the Jansen and Rit neural mass model: analysis and numerics[END_REF]. First, the proved mean-square convergence order 1 is illustrated numerically. Second, the ability of the different numerical methods to preserve the qualitative dynamics of neuronal spiking is analysed, in particular their ability to reproduce the correct amplitudes and frequencies of the underlying oscillations when the time step ∆ is increased. Third, the robustness of the numerical methods to changes in the initial condition X 0 , and how the choice of X 0 may influence the preservation of the phases of the underlying oscillations are analysed. All simulations are carried out in the computing environment R [START_REF]R: A Language and Environment for Statistical Computing[END_REF]. Before we present the simulation results, different Euler-Maruyama type comparison methods, proposed for superlinearly growing coefficients, are recalled.

Revision of Euler-Maruyama type methods

In [START_REF] Hutzenthaler | Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF], it has been shown that the Euler-Maruyama method ( 16) is not mean-square convergent if at least one of the coefficients of the SDE grows superlinearly, as this results in unbounded moments of the iterates. Since then, several explicit variants of this method have been proposed, which aim to control this unbounded growth. The first variant, designed for polynomially growing and one-sided Lipschitz drift and globally Lipschitz diffusion coefficients, has been introduced in [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients[END_REF]. It is based on a taming perturbation which avoids large values caused by the superlinearly growing drift. The method is defined through the iteration

X TEM (t i ) = X TEM (t i-1 ) + F ( X TEM (t i-1 ))∆ 1 + F ( X TEM (t i-1 )) ∆ + Σ √ ∆ψ i-1 , (40)
and proved to be mean-square convergent of order 1/2 (order 1) for SDEs with multiplicative noise (additive noise), i.e., it yields the same convergence rate as achieved by the Euler-Maruyama method in the globally Lipschitz case [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF].

Another variant, aiming to tame both the drift and the diffusion term, has been suggested in [START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF]. The method is defined via

X DTEM (t i ) = X DTEM (t i-1 ) + F ( X DTEM (t i-1 ))∆ + Σ √ ∆ψ i-1 1 + F ( X DTEM (t i-1 )) ∆ + Σ √ ∆ψ i-1 , (41)
and is designed for the broader class of equations where also the diffusion coefficient is allowed to grow polynomially at infinity and satisfies a one-sided Lipschitz condition. It has been shown to converge with mean-square order 1/2 (also in the case of additive noise). For similar variants of the Euler-Maruyama method, see, e.g, [START_REF] Sabanis | Euler approximations with varying coefficients: the case of superlinearly growing diffusion coefficients[END_REF][START_REF] Zhang | Order-preserving strong schemes for SDEs with locally Lipschitz coefficients[END_REF].

The strong convergence (without order) of a different class of variants, based on space truncation techniques, has been discussed in [START_REF] Hutzenthaler | Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF]. In particular, we recall the two methods

X TrEM (t i ) = X TrEM (t i-1 ) + F ( X TrEM (t i-1 ))∆ max 1, F ( X TrEM (t i-1 )) ∆ + Σ √ ∆ψ i-1 , (42)
X DTrEM (t i ) = X DTrEM (t i-1 ) + F ( X DTrEM (t i-1 ))∆ + Σ √ ∆ψ i-1 max 1, ∆ F ( X DTrEM (t i-1 ))∆ + Σ √ ∆ψ i-1 , ( 43
)
constructed to truncate the drift and the drift and diffusion term, respectively. Another type of truncated Euler-Maruyama method, with mean-square convergent rate arbitrarily close to 1, has been proposed in [START_REF] Mao | The truncated Euler-Maruyama method for stochastic differential equations[END_REF][START_REF] Mao | Convergence rates of the truncated Euler-Maruyama method for stochastic differential equations[END_REF]. Here, we recall the partially truncated variant discussed in [START_REF] Guo | The partially truncated Euler-Maruyama method and its stability and boundedness[END_REF]. This method assumes that the drift can be decomposed as

F (X(t)) = F 1 (X(t)) + F 2 (X(t)),
where F 1 is globally Lipschitz continuous and F 2 satisfies Assumption 1. It is given by

X PTrEM (t i ) = X PTrEM (t i-1 ) + F 1 ( X PTrEM (t i-1 )) + F ∆ 2 ( X PTrEM (t i-1 )) ∆ + Σ √ ∆ψ i-1 , (44)
where the function

F ∆ 2 is a truncated version of F 2 .
In particular, it is given by

F ∆ 2 (x) = F 2 min{ x , µ -1 h(∆) } x x ,
where µ : R + → R + such that µ(r) → ∞ as r → ∞ and sup

x ≤r

F 2 (x) ≤ µ(r), ∀ r ≥ 1,
and, for ∆ * ∈ (0, 1], h : (0, ∆ * ] → (0, ∞) such that

h(∆ * ) ≥ µ(1), lim ∆→0 h(∆) = ∞ and ∆ 1/4 h(∆) ≤ 1, ∀∆ ∈ (0, 1).
Thus, the method is not uniquely defined and depends on the choice of µ(•) and h(•). Following [START_REF] Guo | The partially truncated Euler-Maruyama method and its stability and boundedness[END_REF], for the cubic model problem ( 26), we consider F 1 ≡ 0, µ(r) = r 3 and h(∆) = ∆ -1/5 . For this choice, the above conditions on h are satisfied for ∆ * = 1, since h(∆ * ) = µ(1) = 1. Moreover, it holds that µ -1 (h(∆)) = ∆ -1/15 . Numerical experiments for the cubic model problem are reported in Appendix C. For the FHN model ( 2), we consider

F 1 (X(t)) =   1 V (t) -U (t) γV (t) -U (t) + β   , F 2 (X(t)) =   -1 V 3 (t) 0   , µ(r) = r 3 / and h(∆) = ∆ -1/5
. For this choice, the conditions on h are satisfied for ∆ * = 1/ -5 , since then h(∆ * ) = 1/ = µ(1). Therefore, when is small, this method requires very small time steps ∆, and is thus highly inefficient (see the subsequent sections).

In the following, we denote by tamed (TEM), diffusion tamed (DTEM), truncated (TrEM), diffusion truncated (DTrEM) and partially truncated (PTrEM) Euler-Maruyama method, the schemes ( 40), ( 41), ( 42), ( 43) and [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF], respectively. -10 -5 0 q q q q q q q q q q q q q q q q S LT TEM DTEM

TrEM/DTrEM PTrEM Order 1 Order 1/2 log 2 (∆) log 2 (RMSE(∆))
Figure 2: Illustration of the mean-square convergence order on the FHN model ( 2) via the RMSE [START_REF] Mclachlan | Splitting methods[END_REF]. All model parameters are set to 1, X 0 = (0, 0) and T = 5.

Convergence order

The mean-square convergence order can be illustrated by approximating the left-hand side of the inequality in Theorem 1 (for a fixed time T and p = 1) with the root mean-squared error (RMSE) defined by

RMSE(∆) := 1 M M l=1 X l (T ) -X l ∆ (T ) 2 1/2 , ( 45
)
where X l (T ) and X l ∆ (T ) denote the l-th simulated path at a fixed time T of the true process and the approximated process, respectively, for l = 1, . . . , M .

In Figure 2, we report the RMSEs of the different numerical methods in log2 scale as a function of the time step ∆ used to simulate X l ∆ (T ). Since the true process is not known, the reference values X l (T ) are simulated with the TEM method (40) using the small time step ∆ = 2 -14 . We verified that using a different scheme for the simulation of the reference paths does not affect the results of the experiments. The approximated trajectories X l ∆ (T ) are generated with the LT (15), S [START_REF] Hairer | Solving Ordinary Differential Equations: Nonstiff Problems[END_REF], TEM [START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF], DTEM [START_REF] Mao | Stochastic Differential Equations and Applications[END_REF], TrEM [START_REF] Mao | The truncated Euler-Maruyama method for stochastic differential equations[END_REF], DTrEM [START_REF] Mao | Convergence rates of the truncated Euler-Maruyama method for stochastic differential equations[END_REF] and PTrEM (44) method, respectively, under different choices of the time step, namely ∆ = 2 -k , k = 6, . . . , 12. We consider T = 5, M = 10 4 , X 0 = (0, 0) and set all model parameters to 1. All RMSEs are also reported in Table 1. The theoretical convergence order 1, established in Theorem 2, is confirmed numerically. The S splitting yields the smallest RMSEs among all considered numerical methods. The RMSEs of the LT method lie slightly above those obtained under the TrEM and DTrEM methods, which are identical (up to the reported precision). The RMSEs of the tamed Euler-Maruyama methods are larger than those obtained under the splitting, TrEM and DTrEM methods. For the DTEM method we only observe a convergence of order 1/2, in agreement with the observations in [START_REF] Kelly | Adaptive Euler methods for stochastic systems with non-globally Lipschitz coefficients[END_REF][START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF]. The PTrEM method yields the largest RMSEs. However, we observe that for smaller values of σ i , i = 1, 2, the method improves (see also Appendix C, where the impact of the noise intensity on the performance of that method is discussed).

Preservation of neuronal spiking dynamics: amplitudes and frequencies

In the following, we analyse the ability of the considered methods to preserve the qualitative neuronal spiking dynamics of the FHN model. In particular, we investigate whether the amplitudes and frequencies of the neuronal oscillations are kept when increasing the time step ∆. Throughout this and the next section, we omit the DTEM method (41) as it yields a performance comparable to that of the TEM method [START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. Moreover, we set β = σ 1 = 0.1 and σ 2 = 0.2, and consider different values for γ and . These parameters are of particular interest, because they regulate the spiking intensity of the neuron and separate the time scale of the two model components, respectively. When is small, both variables evolve on different time scales. This situation is often referred to as "stiff" case, while larger values of refer to the "nonstiff" case, see, e.g., [START_REF] Chen | Structure-preserving numerical integrators for Hodgkin-Huxley-type systems[END_REF]. Furthermore, these parameters determine the value of e A∆ , and thus the validity of Theorem 4. Since the true process is not available, all reference paths are obtained under the TEM method (40), using the small time step ∆ = 2 • 10 -5 . Also in this case, the choice of the scheme used to simulate the reference paths does not affect the results of the experiments. Moreover, note that all paths are generated using the same set of pseudo random numbers in each example.

In the following, the focus lies on the V -component of the process solving SDE (2), modelling the membrane voltage, which can be experimentally recorded with intracellular measurements. Similar results are obtained for the U -component.

In Figure 3, we report paths of the V -component of the FHN model generated under different values of the time step ∆. An increase in γ leads to an increase in the frequency of the oscillations, and thus in the number of released spikes. Both splitting methods yield almost overlapping paths as ∆ increases, preserving thus the qualitative dynamics of the model, independently of the choice of the intensity parameter γ. In contrast, the TEM method underestimates the frequency and overestimates the amplitude of the neuronal oscillations as ∆ increases, for both values of γ under consideration. For similar observations regarding tamed methods, we refer to [START_REF] Kelly | Adaptive time-stepping strategies for nonlinear stochastic systems[END_REF][START_REF] Kelly | Adaptive Euler methods for stochastic systems with non-globally Lipschitz coefficients[END_REF]. Note also that the paths btained under the TEM method already start deviating from the reference paths for ∆ = 2 • 10 -3 , performing thus worse than the TrEM and DTrEM methods. Since = 0.05, the quantity ∆ * required for the PTrEM method (see Section 7.1) equals 1/ -5 = 3.125•10 -7 . We then observe that this method produces the desired paths only for very small time steps (∆ = 2 • 10 -8 in Figure 3) and fails for the other values of ∆ under consideration.

For a deeper investigation of the neuronal spiking dynamics, we consider the spectral density of the V -component, which takes into account its autocovariance, and thus the dependence of the membrane voltage on previous epochs. It is given by

S V (ν) = F {r V } (ν) = ∞ -∞ r V (τ )e -i2πντ dτ, (46)
where F denotes the Fourier transformation, r V the autocovariance function of (V (t)) t∈[0,T ] and the frequency ν can be interpreted as the number of oscillations in one time unit. We estimate the spectral density S V (ν) with a smoothed periodogram estimator, see, e.g., [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF][START_REF] Quinn | On the periodogram estimators of periods from interleaved sparse, noisy timing data[END_REF], based on paths generated over the time interval [0, 10 3 ]. We use the R-function spectrum and set the required smoothing parameter to span= 0.3T . The estimated spectral densities obtained under different values of γ and different choices of the time step ∆ are reported in Figure 4. As desired, for a fixed γ, all spectral densities estimated from the paths generated under the splitting schemes are almost overlapping as ∆ increases. In contrast, the frequency ν estimated under the Euler-Maruyama type methods decreases as ∆ increases, and the height of the peaks, carrying information about the amplitude of the neuronal oscillations, increases with ∆. Their performance deteriorates as γ increases, the TrEM and DTrEM methods yielding better results than the TEM method. For the considered values of ∆, the spectral densities based on the PTrEM method cannot be derived, because this method produces "NaN" values. This is indicated by the horizontal lines in the bottom left panels. Note also that the estimated frequencies are in agreement with those deduced from Figure 3.

Moreover, the Euler-Maruyama type methods perform even worse in terms of second moment (amplitude) preservation when the parameter is increased, while the splitting methods preserve the qualitative behaviour of the model. This is illustrated in Figure 5, where we increase to 1 (the quantity ∆ * introduced for the PTrEM method in Section 7.1 thus equals 1), fix γ = 20 and report phase portraits of the system obtained under the different numerical methods for ∆ = 2 • 10 -4 and ∆ = 2 • 10 -2 . Again, the splitting methods preserve the behaviour of the process (X(t)) t∈[0,T ] as ∆ increases, while the Euler-Maruyama type methods produce larger orbits, overshooting the second moment of the process.

In addition, we investigate the ability of the considered numerical methods to approximate the underlying invariant density of the process (X(t)) t∈[0,T ] . In particular, we estimate the marginal invariant density of the V -component of the FHN model with a standard kernel density estimator given by

π V (v) = 1 nH n i=1 K v -V (t i ) H , ( 47
)
where H is a smoothing bandwidth and K a kernel function [START_REF] Pons | Functional Estimation for Density, Regression Models and Processes[END_REF]. Taking advantage of the ergodicity of the FHN model, the sample V (t i ), i = 1, . . . , n, in ( 47) is obtained from a long-time simulation of a single path. We use the R-function density, a kernel estimator as described in [START_REF] Milstein | A theorem on the order of convergence of mean-square approximations of solutions of systems of stochastic differential equations[END_REF].

In Figure 6, we report the marginal invariant densities of the process (V (t)) t∈[0,T ] estimated via [START_REF] Milstein | A theorem on the order of convergence of mean-square approximations of solutions of systems of stochastic differential equations[END_REF] based on paths generated over the time interval [0, 10 4 ], for = 1, γ = 20 and different values of ∆. Both splitting methods yield reliable estimates for all values of ∆ under consideration. In contrast, the densities obtained under the Euler-Maruyama type methods already deviate from the desired ones for ∆ = 2 • 10 -3 , and suggest a transition from a unimodal to a bimodal density when ∆ is further increased to 2 • 10 -2 . It is again visible that the Euler-Maruyama type methods overestimate the second moment, and thus the amplitudes of the process. Similar results are also obtained for the U -component.

Impact of the initial condition: preservation of phases

Finally, we compare the considered numerical methods regarding their sensitivity to changes in the initial condition X 0 . In particular, we illustrate that, when V 0 is large, i.e., when the process starts far away from the mean of the invariant distribution, the considered Euler-Maruyama type methods do not correctly reproduce the phases of the underlying oscillations, even when the time step ∆ is very small. In contrast, the splitting methods are less sensitive to changes in the initial condition. Similar observations are made when U 0 is large (figures not shown).

The impact of V 0 on the performance of the different numerical methods is shown in Figure 7 and Figure 8, where we report paths of the V -component, simulated under ∆ = 2 • 10 -4 , U 0 = 0 and different values of V 0 . The grey reference path is simulated under ∆ = 2 • 10 -7 using the TEM method [START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. As before, the results are not influenced by the choice of the numerical method used to generate the reference paths. The underlying parameter values are the same as in Section 7.3, choosing γ = 5 and = 0.05 in Figure 7, and γ = 20 and = 1 in Figure 8. As desired, the splitting methods are barely influenced by V 0 , even when it is very large, with paths overlapping with the reference paths for all t under consideration. In contrast, when V 0 is large, the Euler-Maruyama type methods introduce a delay in when the generated paths reach the oscillatory dynamics, this behaviour deteriorating as V 0 increases. Moreover, they also do not preserve the phases of the oscillations, introducing a shift. In Figure 7, the DTrEM method reaches the correct oscillatory dynamics, though shifted, almost as fast as the splitting methods for V 0 = 10 4 , but fails to reach the invariant regime for V 0 = 10 3 . In Figure 8, it does not enter the invariant regime for both V 0 = 10 3 and V 0 = 10 4 . Moreover, spurious oscillations produced by the DTrEM method were obtained for other parameter combinations, as also observed in [START_REF] Kelly | Adaptive Euler methods for stochastic systems with non-globally Lipschitz coefficients[END_REF][START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF]. For = 0.05 (see Figure 7), the PTrEM method does not produce the desired paths, even when V 0 is close to the invariant mean. For = 1 (see Figure 8), it yields the correct path when V 0 = 1, a path which initially deviates from the others when V 0 = 3, and produces high-amplitude oscillations, not entering the invariant regime, when V 0 = 10 3 and V 0 = 10 4 . Therefore, the PTrEM method reacts very sensitively to the choice of X 0 , this undesired behaviour being also observed for the cubic model problem [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF] introduced in Section 5, see Appendix C. Remark 6. Note that the only considered combination of γ and in this section for which Assumption 4, and thus Theorem 4 holds is γ = 1/ = 20. However, we do not observe a difference in the quality of the splitting methods depending on the combination of these parameters. Intuitively, this is because the underlying linear SDE (9) with matrix A as in (36), i.e., the damped stochastic oscillator, is geometrically ergodic for all γ > 0 and > 0.

Conclusion and discussion

We propose a splitting strategy to approximate the solutions of semi-linear SDEs with additive noise and globally one-sided Lipschitz continuous drift coefficients which are allowed to grow polynomially. We prove that the resulting explicit Lie-Trotter splitting method is mean-square convergent of order 1. In contrast to existing explicit mean-square convergent Euler-Maruyama type methods, which may also achieve a convergence rate of order 1, the constructed method preserves important structural properties of the model. First, it provides a more accurate approximation of the noise structure of the SDE through the covariance matrix of the exact solution of the stochastic subequation. In particular, while the conditional covariance matrix of Euler-Maruyama type methods only contains the information of the diffusion matrix Σ, the splitting method also relies on the matrix A in the semi-linear drift. This is particularly beneficial when the SDE is hypoelliptic. Indeed, while the conditional covariance matrix of the existing methods is degenerate in that case, we establish the desired 1-step hypoellipticity of the constructed splitting method, meaning that it admits a smooth transition density in every iteration step. In particular, the method yields non-degenerate Gaussian transition densities, a feature which is advantageous within likelihood-based estimation techniques, where the existing numerical methods cannot be applied [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF][START_REF] Melnykova | Parametric inference for hypoelliptic ergodic diffusions with full observations[END_REF][START_REF] Pokern | Parameter estimation for partially observed hypoelliptic diffusions[END_REF].

Second, Euler-Maruyama type methods do not preserve the geometric ergodicity of the process. As a consequence, they are not robust to changes in the initial condition, yield poor approximations of the underlying invariant distribution, or do not preserve the moments of the process. In contrast, the proposed splitting method is proved to preserve the Lyapunov structure of SDE, as long as an assumption on the solution f of the deterministic subequation is satisfied and it holds that e A∆ < 1 for all ∆ ∈ (0, ∆ 0 ]. If, in addition, the logarithmic norm µ(A) < 0, the method is proved to have an asymptotically bounded second moment. In the one-dimensional case, a precise bound of the second moment of the splitting method is derived in closed-form and illustrated on a cubic model problem. We also consider the FHN model, a well known equation used to describe the firing activity of single neurons. The geometric ergodicity of the proposed splitting method applied to this equation is established under a restricted parameter space.

Third, we illustrate on the FHN model that, in contrast to Euler-Maruyama type methods, the proposed splitting method preserves the amplitudes, frequencies and phases of neuronal oscillations, even for large time steps. This may make the method particularly beneficial when used, e.g., to simulate large networks of neurons, or when embedded within simulation-based inference procedures. Besides the Lie-Trotter splitting method, we also consider a method which is based on a Strang composition in our numerical experiments. Both splitting methods perform comparably good throughout, the Strang splitting behaving slightly better in some scenarios. As the considered Euler-Maruyama type methods do converge, their lack of structure preservation becomes less visible when using very small time steps. However, the use of significantly smaller time steps results in drastically higher computational costs, making these methods highly inefficient and, consequently, computationally infeasible within simulation-based inference algorithms, as previously illustrated in [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF].

Several generalisations of the considered approach are possible. The proposed splitting strategy can be, e.g., applied to the stochastic Van der Pol oscillator [START_REF] Van Der Pol | A theory of the amplitude of free and forced triode vibrations[END_REF][START_REF] Van Der Pol | On "relaxation-oscillations[END_REF], whose investigation leads to similar numerical results. The presented approach may be also applied to SDEs (1) with other types of nonlinearity. In particular, as long as the ODE determined by the function N is exactly solvable (see [START_REF] Kamke | Differentialgleichungen Lösungsmethoden und Lösungen: Gewöhnliche Differentialgleichungen[END_REF] for diverse solution methods) and satisfies some useful conditions, the method presented in this article may be used. Moreover, the proposed method may be extended to SDEs with multiplicative noise, e.g., to Σ(X(t)) = σX(t), σ > 0, where the stochastic subequation of the splitting framework corresponds to the geometric Brownian motion. This may be relevant, e.g., for the stochastic Ginzburg-Landau equation arising from the theory of superconductivity [START_REF] Ginzburg | On the theory of superconductivity[END_REF][START_REF] Hutzenthaler | Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF]. Furthermore, the investigation of conditions under which the presented results are still valid when the solution of ODE [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] is not available exactly, constitutes another topic for future research.

≤ 2(x -y) 2 + 2(y 3 -x 3 ) 2 = 2(x -y) 2 + 2(x -y) 2 (y 2 + xy + x 2 ) 2 ≤ 2(x -y) 2 + 9(x -y) 2 (x 4 + y 4 ),
where we used that 2xy ≤ x 2 + y 2 and that (3x 2 /2 + 3y 2 /2) 2 ≤ 9(x 4 + y 4 )/2 in the last inequality. Thus, we obtain that

N (x) -N (y) 2 ≤ 9(x -y) 2 1 + x 4 + y 4 ,
which proves that the assumption holds for c 2 = 9 and χ = 3. where denotes the derivative with respect to t. Consider two cases. First, let y / ∈ (0, 1). Then it holds that g (y; t) ≤ 0 for all t ≥ 0. Second, let y ∈ (0, 1). To prove that g (y; t) ≤ 1/2, we determine the global maximum of g (y; t) with respect to t. Solving g (y; t) = 0 with respect to t, gives that

t max = 1 2 log 1 y -1 .
Noting that t max exists and that g (y; t max ) = 1/2 for any y ∈ (0, 1) proves the result.

Remark 7. If y ∈ (0, 1), it also holds that g (y; t) ≤ 2e 2t y(1 -y) ≤ 1 2 e 2t ≤ e 2∆0 2 .

Thus, a simpler argument suffices to prove the statement for c 3 (∆ 0 ) = e 2∆ 0 2 > 1 2 . Assumptions 4 and 5: These statements are satisfied because A = -1, and thus the matrix norm e A∆ = e -∆ < 1 and the logarithmic norm µ(A) = -1 < 0.

B Proof of Proposition 3

Proof. Assumption 1: Denote x = (v 1 , u 1 ) and y = (v 2 , u 2 ) . We start with Assumption (A1) and have that (N (x) -N (y),

x -y) = 1 (v 1 -v 2 ) 2 1 -(v 2 1 + v 1 v 2 + v 2 2 ) ≤ 1 (v 1 -v 2 ) 2 ≤ 1 x -y 2 .
Thus, the assumption holds for c 1 = 1/ . Now, consider Assumption (A2). Applying similar arguments as in Appendix A, we have that

N (x) -N (y) 2 = 1 (v 1 -v 2 ) + 1 (v 3 2 -v 3 1 ) 2 ≤ 2 2 (v 1 -v 2 ) 2 + 9 2 (v 1 -v 2 ) 2 (v 4 1 + v 4 2 ).
Using that (v 1 -v 2 ) 2 ≤ x -y 2 and that v 4 1 + v 4 2 ≤ x 4 + y 4 , we finally obtain that

N (x) -N (y) 2 ≤ 9 2 x -y 2 1 + x 4 + y 4 .
Thus, the assumption holds for c 2 = 9/ 2 and χ = 3. Assumption 2: Condition (6) holds for the linear SDE (9), since

∂ u (Ax) 1 σ 2 = - σ 2 = 0.
Thus, the equation is hypoelliptic. It holds that λ 2 (∆) < 0 for all ∆ > 0, where denotes the derivative with respect to ∆. Thus, λ 2 (∆) is strictly decreasing in ∆. Noting that λ 2 (0) = 1 implies the statement. When κ < 0, γ < 1/2. In this case, the eigenvalues of B are given by

λ 1 (∆, γ) = e -∆ κ 4γ 2 -cosh( √ -κ∆) -2 1 -8γ 2 + cosh( √ -κ∆) sinh 2 ( √ -κ∆/2) , λ 2 (∆, γ) = e -∆ κ 4γ 2 -cosh( √ -κ∆) + 2 1 -8γ 2 + cosh( √ -κ∆) sinh 2 ( √ -κ∆)/2 ,
where λ 1 (∆, γ) ≤ λ 2 (∆, γ) for all ∆ > 0 and γ < 1/2. For γ < 1/2 arbitrary, but fixed, the partial derivative of λ 2 (∆, γ) with respect to ∆ exists and satisfies ∂ ∂∆ λ 2 (∆, γ) < 0, ∀ ∆ ∈ (0, ∆ 0 ].

Thus, the function λ 2 (∆, γ) is strictly decreasing in ∆. Moreover, we have that λ 2 (0, γ) = 1 for any γ, which implies the statement. When κ > 0, γ > 1/2. In this case, the eigenvalues of B are given by Again, we observe that λ 1 (∆, γ) ≤ λ 2 (∆, γ) for all ∆ > 0 and γ > 1/2. Consider γ > 1/2 arbitrary, but fixed. Moreover, define I ∆ := {2πl/ √ κ, l ∈ N}. Since cos(2πl) = 1 and sin(πl) = 0, we have that λ 1 (γ, ∆) = λ 2 (γ, ∆) = e -∆ < 1, ∀ ∆ ∈ I ∆ .

λ 1 (∆, γ) = e -∆ κ 4γ 2 -cos( √ κ∆) -2 -1 + 8γ
Let now ∆ ∈ (0, ∞)\I ∆ . For those values of ∆, the partial derivative of λ 2 (∆, γ) with respect to ∆ exists. In particular, we have that

∂ ∂∆ λ 2 (∆, γ) < 0, ∀ ∆ ∈ (0, ∞)\I ∆ .
Thus, for a fixed γ, the function λ 2 (∆, γ) is strictly decreasing in ∆. Noting that λ 2 (0, γ) = 1 for any γ implies the statement.

C Numerical experiments for the cubic model problem

Consider the cubic model problem [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF] introduced in Section 5. We now illustrate how the choice of X 0 influences the behaviour of paths of the ergodic process X(t) simulated under the different numerical methods. If X 0 is large compared to the invariant mean, the standard Euler-Maruyama method ( 16) produces paths which are computationally pushed to +/-infinity within a few iteration steps, even for very small values of ∆. This is not the case for the tamed/truncated variants of this method. However, they may also react sensitively to X 0 , even for small ∆. This is illustrated in Figure 9, where we report paths of SDE (26) generated for different values of X 0 , using ∆ = 10 -4 and σ = 1/2. The grey reference paths are simulated under ∆ = 10 -7 using the TEM method [START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. The choice of the reference method does not change the reported results, and all paths are generated under the same underlying pseudo random numbers. Note that, the DTEM method ( 41) is not reported in Figure 9, because it shows a performance comparable to that of the TEM method [START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. As desired, the paths obtained under the splitting methods ( 15) and ( 23) are not deterred by large values of X 0 , and overlap with the reference path for all values of X 0 under consideration. In contrast, for large values of X 0 , the Euler-Maruyama type methods introduce a delay in when the respective paths reach the reference path. This behaviour deteriorates as X 0 increases. The path obtained under the PTrEM method (44) initially deviates from the desired one, even when X 0 = 5, not reaching the reference path for the values of t under consideration for X 0 = 10 4 and X 0 = 3 • 10 4 . Note also that, for some values of X 0 , we observe that the DTrEM method (43) may produce spurious oscillations (figures not shown). See [START_REF] Kelly | Adaptive Euler methods for stochastic systems with non-globally Lipschitz coefficients[END_REF][START_REF] Tretyakov | A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications[END_REF], where such a behaviour has also been observed. and different values of σ. The grey reference paths are obtained under ∆ = 10 -7 using the TEM method [START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. All paths correspond to the same random realisation.

In addition, we observe that the PTrEM method (44) may also produce paths which deviate from the desired ones for larger values of the noise parameter σ. This is illustrated in Figure 10, where we report paths of SDE (26) generated under X 0 = 0, ∆ = 10 -4 and different values of σ. While the paths obtained under the LT splitting (15) (the same is observed for all other schemes except for the PTrEM method) overlap with the reference paths for both values of σ under consideration, the PTrEM method produces a path which deviates from the desired one when σ = 4 (right panel). This behaviour deteriorates as σ increases.

Figure 1 :

 1 Figure 1: Bound K LT (t i , X 0 ) (as a function of time) and asymptotic bound K LT ∞ of Corollary 2 for SDE (26) (a = 1 and c 3 = 1/2) with σ = 1/2 and X 0 = 2, and estimate of E[X 2 (t i )] obtained from 10 4 paths generated under the LT splitting.

Figure 3 :

 3 Figure 3: Paths of the V -component of the FHN model (2) simulated under the considered numerical methods for X 0 = (-1, 0) , β = σ 1 = 0.1, σ 2 = 0.2, = 0.05, two different values of γ and increasing time step ∆. Paths of the PTrEM method are also generated under a smaller time step than used for the other methods, i.e., ∆ = 2 • 10 -8 . All paths correspond to the same random realisation.

Figure 4 :

 4 Figure 4: Estimates of the spectral density (46) of the V -component of the FHN model (2) obtained under the considered numerical methods for X 0 = (0, 0) , β = σ 1 = 0.1, σ 2 = 0.2, = 0.05, two different values of γ and increasing time step ∆.

Figure 5 :

 5 Figure 5: Phase portraits of the FHN model (2) simulated under the considered numerical methods for X 0 = (0, 0) , β = σ 1 = 0.1, σ 2 = 0.2, = 1, γ = 20 and increasing time step ∆. All paths correspond to the same random realisation.

Figure 6 :

 6 Figure 6: Estimates of the invariant density (47) of the V -component of the FHN model (2) obtained under the considered numerical methods for X 0 = (0, 0) , β = σ 1 = 0.1, σ 2 = 0.2, = 1, γ = 20 and increasing time step ∆.

Figure 7 :Figure 8 :

 78 Figure 7: Paths of the V -component of the FHN model (2) simulated under the considered numerical methods for different values of V 0 (U 0 = ∆ = 2 • 10 -4 , β = σ 1 = 0.1, σ 2 = 0.2, γ = 5 and = 0.05. The grey reference paths are obtained under ∆ = 2 • 10 -7 using the TEM method[START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. All paths correspond to the same random realisation.

Assumption 2 : 3 : 1 2 =

 2312 Since d = 1, this is clear. Assumption We prove the statement for c 3 = 1/2. Setting y = x 2 , it has to be shown thatf 2 (x; t) = g(y; t) := y e -2t + y(1 -e -2t ) ≤ y + 1 2 t =: h(y; t), ∀t ∈ (0, ∆ 0 ].Since g(y; 0) = h(y; 0) = y, it suffices to prove that for any y ∈ R + 0 it holds that g (y; t) = -2e 2t (y -1)y(1 + (e 2t -1)y) 2 ≤ h (y; t), ∀t ∈ (0, ∆ 0 ],

Assumption 3 : 2 ∆

 32 We have that f (x; ∆) = (f 1 (v; ∆), f 2 (u; ∆)) .Consider the V -component. The fact that, for any v ∈ R it holds thatf 2 1 (v; ∆) ≤ v 2 + 1 ∀ ∆ ≥ 0,can be proved in the same way as in Appendix A. Regarding the U -component, by assumption, we have thatf 2 2 (u; ∆) = u 2 . Thus, f (x; ∆) 2 = f 2 1 (v; ∆) + f 2 2 (v; ∆) ≤ v 2 + u 2 statement for c 3 = 1/(2 ).Assumption 4: Recall that e A∆ = λ max ((e A∆ ) (e A∆ )), and define B := (e A∆ ) (e A∆ ). It suffices to prove that λ max (B) < 1 for all ∆ ∈ (0, ∆ 0 ].Since by assumption γ = 1/ , κ defined in (37) becomes κ = 4γ 2 -1. When κ = 0, this condition is equivalent to γ = 1/2. In this case, the eigenvalues of B are given byλ 1 (∆) = 1 2 e -∆ 2 + ∆ 2 -∆ 2 (4 + ∆ 2 ) ≤ λ 2 (∆) = 1 2 e -∆ 2 + ∆ 2 + ∆ 2 (4 + ∆ 2 ) .

Figure 9 :Figure 10 :

 910 Figure 9: Paths of SDE (26) simulated under the considered numerical methods for different values of X 0 , ∆ = 10 -4 and σ = 1/2. The grey reference paths are obtained under ∆ = 10 -7 using the TEM method[START_REF] León | Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate[END_REF]. All paths correspond to the same random realisation.

Table 1 :

 1 RMSE (45) for different values of ∆. All model parameters are set to 1, X 0 = (0, 0) and T = 5.

	∆	S	LT	TEM	DTEM	TrEM	DTrEM PTrEM
	2 -6	0.01320 0.01751 0.03628 0.27208 0.01667 0.01667 2.23688
	2 -7	0.00659 0.00871 0.01795 0.20075 0.00824 0.00824 1.41733
	2 -8	0.00323 0.00431 0.00889 0.14726 0.00410 0.00410 0.83296
	2 -9	0.00161 0.00213 0.00438 0.10639 0.00204 0.00204 0.45305
	2 -10 0.00079 0.00106 0.00213 0.07664 0.00101 0.00101 0.22307
	2 -11 0.00039 0.00053 0.00099 0.05446 0.00050 0.00050 0.10771
	2 -12 0.00020 0.00027 0.00043 0.03877 0.00027 0.00027 0.05273

E.B. was supported by the LCM -K2

Center within the framework of the Austrian COMET-K2 program. A.S. was supported by MIAI@Grenoble Alpes, (ANR-19-P3IA-0003). E.B., M.T. and I.T. were supported by the Austrian Science Fund (FWF): W1214-N15, project DK14. All authors were supported by the Austrian Exchange Service (OeAD), bilateral project FR 03/2017.

Appendix

A Proof of Proposition 1

Proof. Assumption 1: We start with Assumption (A1) and have that

Thus, the assumption holds for c 1 = 1, see also Example 1.2.16 in [START_REF] Humphries | Deterministic and random dynamical systems: Theory and numerics[END_REF]. Now, consider Assumption (A2). We have that N (x) -N (y)