# Asymptotic expansion of the invariant measure for Markov-modulated ODEs at high frequency

Pierre Monmarché, Edouard Strickler

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Monmarché, Edouard Strickler. Asymptotic expansion of the invariant measure for Markov-modulated ODEs at high frequency. 2024. hal-04779988

HAL Id: hal-04779988

https://hal.inrae.fr/hal-04779988v1

Preprint submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Asymptotic expansion of the invariant measure for Markov-modulated ODEs at high frequency

Pierre Monmarché*, Édouard Strickler†

September 29, 2023

## Abstract

We consider time-inhomogeneous ODEs whose parameters are governed by an underlying ergodic Markov process. When this underlying process is accelerated by a factor $\epsilon^{-1}$, an averaging phenomenon occurs and the solution of the ODE converges to a deterministic ODE as $\epsilon$ vanishes. We are interested in cases where this averaged flow is globally attracted to a point. In that case, the equilibrium distribution of the solution of the ODE converges to a Dirac mass at this point. We prove an asymptotic expansion in terms of $\epsilon$ for this convergence, with a somewhat explicit formula for the first order term. The results are applied in three contexts: linear Markov-modulated ODEs, randomized splitting schemes, and Lotka-Volterra models in random environment. In particular, as a corollary, we prove the existence of two matrices whose convex combinations are all stable but such that, for a suitable jump rate, the top Lyapunov exponent of a Markov-modulated linear ODE switching between these two matrices is positive.

## 1 Introduction

Let $\mathcal{M}$ be a $d$-dimensional compact $\mathcal{C}^\infty$ manifold. Markov-modulated ODEs on $\mathcal{M}$ are dynamical systems $(x(t))_{t \ge 0}$ solution to

$$x'(t) = F_{\sigma(t)}(x(t)) \quad (1)$$

where $(\sigma(t))_{t \ge 0}$ is a Markov process on some space $\mathcal{S}$ and, for all $s \in \mathcal{S}$, $F_s$ is a vector field on $\mathcal{M}$.

We are interested in the high frequency regime, namely when $\sigma(t)$ is replaced in (1) by $\sigma(t/\epsilon)$ for some small $\epsilon$. Provided $\sigma$ is ergodic with respect to some

<footnote>*Sorbonne Université, 4 place Jussieu 75005 Paris (France)</footnote>

<footnote>†Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy (France)</footnote>
probability measure $\pi$, the modulated process $x(t)$ is known to converge, as $\varepsilon$ vanishes, to the solution of the averaged ODE

$$x'(t) = \bar{F}(x(t)), \quad \text{where} \quad \bar{F}(x) = \int_S F_s(x)\pi(ds). \qquad (2)$$

In the case where this averaged ODE has a unique global attractor $\bar{x}$, we expect the solution of (1) to be close to $\bar{x}$ for $t$ large and $\varepsilon$ small. In other words, the invariant measure of the Markov process $(x(t), \sigma(t/\varepsilon))_{t \ge 0}$ converges, as $\varepsilon$ vanishes, to $\delta_{\bar{x}} \otimes \pi$. The goal of the present work is to provide an infinitesimal expansion in terms of $\varepsilon$ for this invariant measure. This is obtained by combining a similar expansion for the law of the process for a fixed $t > 0$ (following [22]) with a long-time convergence result for the limit process. The structure of the proof follows the work of Talay and Tubaro [23] who established a similar expansion for the invariant measure of Euler-Maruyam schemes for diffusion processes (the averaging estimates for a fixed $t > 0$ being in that case replaced by finite-time discretization error expansions).

One of our main motivation is to get the first-order term in the expansion of the top Lyapunov exponent of systems of switched linear ODE [5, 16, 9], as detailed in Section 3.1. Our result also has applications for instance in some population models [21, 6, 4] or random splitting numerical schemes [1], cf. Section 3. More generally, Markov-modulated ODEs form a flexible class of Markov processes which appear in a variety of models, often to describe systems evolving in a randomly fluctuating environment, as in finance [12], biology [17] or fiability [19]. It is related to random ODEs, see e.g. [2, 10] and references within.

In the specific context of Piecewise-Deterministic Markov processes (PDMP; namely, when $\sigma$ is a Markov chain on a finite set), the question of fast averaging has been addressed in [13] where a large deviation principle is proven, in [22] where the expansion of the law of the process at fixed time $t$ is given, or in [8] where the convergence of the invariant measure of the Markov process $(x(t), \sigma(t/\varepsilon))_{t \ge 0}$ is proven. In the recent paper [15], the authors deal with the case where $\sigma$ is not ergodic.

The paper is organized as follows. In the remaining of this introduction, we introduce our general notations and assumptions. Our main result, Theorem 1, is stated and proven in Section 2. Section 3 is devoted to examples of applications. Finally, an appendix gathers the proofs of some intermediary results.

To conclude this section, let us clarify our settings. Denote by $(Q_t)_{t \ge 0}$ the semi-group associated to $(\sigma_t)_{t \ge 0}$, namely

$$Q_t f(\sigma) = \mathbb{E}_{\sigma}(f(\sigma_t)),$$

and by $Q$ its infinitesimal generator.

Denote by $\mathcal{A}$ the set of bounded measurable functions from $\mathcal{M} \times \mathcal{S}$ to $\mathbb{R}$ which are $\mathcal{C}^{\infty}$ in their first variable $x$ and such that all their derivatives in $x$ are bounded over $\mathcal{M} \times \mathcal{S}$, and which moreover are such that for all $x \in \mathcal{M}$ and
all multi-index $\alpha$, $s \mapsto \partial_x^\alpha f(x, s)$ is in the domain of $Q$. We consider on $\mathcal{A}$ the norms

$$ \|f\|_j = \sum_{|\alpha| \le j} \|\partial_x^\alpha f\|_\infty. $$

**Assumption 1.** We assume that:

1. Almost surely, the ODE (1) is well-defined for all times, with values in $M$. Moreover, the coordinates of $(x, s) \mapsto F_s(x)$ are in $\mathcal{A}$.

2. The semigroup $(Q_t)_{t \ge 0}$ admits a unique invariant probability measure $\pi$. Moreover, there exist $C, \gamma > 0$ such that

$$ \|Q_t f - \pi f\|_{\infty} \le C e^{-\gamma t} \|f\|_{\infty} \quad (3) $$

for all bounded measurable $f$ on $S$ and all $t \ge 0$.

The uniform exponential ergodicity (3) implies that the operator $Q^{-1}$ given by

$$ Q^{-1}f := \int_{0}^{\infty} Q_t (\pi f - f) \, dt \quad (4) $$

is well-defined for all $f \in \mathcal{A}$, with $Q^{-1}f \in \mathcal{A}$ and $\|Q^{-1}f\|_j \le C/\gamma \|f\|_j$ for all $j \in \mathbb{N}$. Moreover, using that $\partial_t Q_t f = Q Q_t f = Q_t Q f$ for all $t \ge 0$, we get that $QQ^{-1}f = Q^{-1}Qf = f - \pi f$ for $f \in \mathcal{A}$, namely $Q^{-1}$ is the pseudo-inverse of $Q$.

**Example 1.** A case of interest is given by $Qf = \pi f - f$, which corresponds to $\sigma(t) = Y_{N_t}$ where $(N_t)_{t \ge 0}$ is a standard Poisson process with intensity 1, $Y_0 = \sigma(0)$ and $(Y_k)_{k \ge 1}$ is an i.i.d. sequence of random variables distributed according to $\pi$. It is readily checked that, in that case, $Q_t f = e^{-t}f + (1 - e^{-t})\pi f$ and then $Q^{-1} = Q$.

**Example 2.** We will be particularly concerned with the case when $S$ is a finite state space and $(\sigma_t)_{t \ge 0}$ is an irreducible continuous-time Markov chain on $S$, with transition rate matrix $Q$. In that case, $Q^{-1}$ is the group inverse of $Q$, defined as the unique matrix $X$ solution<fnref n="1" /> to

$$ QXQ = Q, \quad XQX = X, \quad XQ = QX. \quad (5) $$

In the particular case when the cardinal of $S$ is 2, then $Q$ can be written as

$$ Q = \begin{pmatrix} -p & p \\ q & -q \end{pmatrix}, $$

for some $p, q > 0$. Moreover, one can easily check, using $Q^2 = -(p+q)Q$ that $X = \frac{1}{(p+q)^2}Q$ satisfies (5), and therefore $Q^{-1} = \frac{1}{(p+q)^2}Q$.

<footnote n="1">¹such a matrix exists and is unique since $Q$ has index 1, meaning that $Q^2$ and $Q$ have the same rank, because 0 is a simple eigenvalue of $Q$. See Chapter 4 in [3] for more details.</footnote>
## 2 High-frequency expansion

The process $(X(t), \sigma(t/\varepsilon))_{t \ge 0}$ is a Markov process on $\mathcal{M} \times \mathcal{S}$ with generator $L = L_c + \frac{1}{\varepsilon}Q$, where $Q$ is seen as an operator on functions on $\mathcal{M} \times \mathcal{S}$ acting only on the second variable and

$$L_c f(x, s) = F_s(x) \cdot \nabla_x f(x, s),$$

for all $f \in \mathcal{A}$. For $t \ge 0$, denote by $P_t^\varepsilon$ the associated Markov semigroup, namely

$$P_t^\varepsilon f(x, s) = \mathbb{E}_{x,s} (f(X(t), \sigma(t/\varepsilon))).$$

We start by stating an expansion in $\varepsilon$ of $P_t^\varepsilon$ for a fixed $t$ (the proof is postponed to the Appendix). This is essentially the result (and proof) of [22] but written in a dual form and with more explicit functional settings. Moreover, [22] only considers the case where $\mathcal{S}$ is a finite set, but it allows the jump rates of $\sigma$ to depend on $x$.

**Proposition 1.** There exist two families of operators $P_t^{(k)}$, $S_t^{(k)}$ indexed by $t \ge 0$, $k \in \mathbb{N}$, acting on $\mathcal{A}$, with the following properties.

* For all $k, j \in \mathbb{N}$ and $T > 0$, there exists $C > 0$ such that for all $f \in \mathcal{A}$,

$$\sup_{t \in [0, T]} \left\| P_t^{(k)} f \right\|_j \le C \|f\|_{j+2k} \quad (6)$$

* For all $k, j \in \mathbb{N}$, there exist $C, \gamma > 0$ such that for all $t \ge 0$ and all $f \in \mathcal{A}$,

$$\left\| S_t^{(k)} f \right\|_j \le C e^{-\gamma t} \|f\|_{j+2k} \quad (7)$$

* For all $n \in \mathbb{N}$, $T \ge 0$, $j \in \mathbb{N}$, there exists $C > 0$ such that for all $\varepsilon > 0$ and $f \in \mathcal{A}$, the remainder $R_{n,t}^\varepsilon f$ defined by

$$R_{n,t}^\varepsilon f = P_t^\varepsilon f - \sum_{k=0}^n \varepsilon^k P_t^{(k)} f - \sum_{k=0}^n \varepsilon^k S_t^{(k)} f, \quad (8)$$

satisfies

$$\left\| R_{n,t}^{\varepsilon} f \right\|_{j} \le C \varepsilon^{n+1} \left\| f \right\|_{j+3+2n} . \quad (9)$$

* The first operators are given by

$$P_t^{(0)} f(x, s) = \pi f(\bar{\varphi}_t(x)), \quad S_t^{(0)} f(x, s) = Q_t f(x, s) - \pi f(x) \quad (10)$$

and

$$P_t^{(1)} f(x, s) = \mathbf{b}_1(t, x, s) + \int_0^{+\infty} \pi(L_c S_r^{(0)} f)(x) dr + \int_0^t \pi(L_c \mathbf{b}_1)(r, \bar{\varphi}_{t-r}(x)) dr, \quad (11)$$

where

$$\mathbf{b}_1(t, x, s) = Q^{-1}(\partial_t P_t^{(0)} f - L_c P_t^{(0)} f)(x, s). \quad (12)$$
Next, to address the question of the long-time behaviour of the process, we work under the following condition:

**Assumption 2.** There is a globally attractive $\bar{x} \in \mathcal{M}$ for the averaged flow $\bar{F}$ given in (2), and for all $j \in \mathbb{N}$, there exist $C, a > 0$ such that for all $t \ge 0$ and $f \in C^\infty(\mathcal{M})$,

$$ \|f \circ \varphi_t - f(\bar{x})\|_j \le C e^{-at} \|f\|_{j+1}, $$

where $\varphi$ is the flow associated to $\bar{F}$.

Under Assumption 2, $\mu_0 = \delta_{\bar{x}} \otimes \pi$ is the unique invariant measure of $P_t^{(0)}$ given in (10). Our main result is the following expansion in terms of $\varepsilon$ of any invariant measure of the process $(P_t^\varepsilon)_{t \ge 0}$, in the spirit of Talay-Tubaro expansions in terms of the step size for discretization schemes [23].

**Theorem 1.** Under Assumptions 1 and 2, for all $f \in \mathcal{A}$ there exists real sequences $(c_k)_{k \ge 1}$ and $(M_k)_{k \ge 1}$ such that for all $n \in \mathbb{N}$ and any invariant measure $\mu_\varepsilon$ of $(P_t^\varepsilon)_{t \ge 0}$,

$$ \left| \mu_{\varepsilon} f - \mu_{0} f - \sum_{k=1}^{n} c_{k} \varepsilon^{k} \right| \le M_{n} \varepsilon^{n+1}. \qquad (13) $$

Moreover,

$$ c_1 = \pi L_c Q^{-1} (L_c h - f)(\bar{x}) \qquad (14) $$

where, for $(x, s) \in \mathcal{M} \times \mathcal{S}$,

$$ h(x, s) = \int_{0}^{\infty} \pi (f(\bar{x}) - f(\varphi_{r}(x))) dr. \qquad (15) $$

**Remark 1.** In dimension 1, it is possible to get an alternative expression for the function $h$ appearing in formula (14). Indeed, for fixed $x \in \mathbb{R}$, with the change of variables $y = \varphi_r(x)$; we obtain, due to $\frac{d\varphi_r(x)}{dr} = \bar{F}(\varphi_r(x))$,

$$ h(x, s) = \int_{0}^{+\infty} (\pi f(\bar{x}) - \pi f(\varphi_r(x))) dr = \int_{\bar{x}}^{x} \frac{\pi f(y) - \pi f(\bar{x})}{\bar{F}(y)} dy. $$

We will use this expression in Section 3.3

*Proof of Theorem 1.* The proof is divided in two parts. In the first one, we prove (13), and we give an expression of $c_1$ which depends on an arbitrary time $t > 0$. In the second part, we obtain the announced expression of $c_1$ by letting $t$ vanish in this first expression.

**Step 1.** Fix $t > 0$, $f \in \mathcal{A}$ and let $\mu_\varepsilon$ be an invariant measure of $(P_s^\varepsilon)_{s \ge 0}$. By definition of the flow $\varphi_t$,

$$ (P_t^{(0)})^n f(x, i) = \pi f(\varphi_t^{\circ n}(x)) = \pi f(\varphi_{tn}(x)). $$

Hence,

$$ \|(P_t^{(0)})^n (\mu_0 f - f)\|_j = \|\pi f(\bar{x}) - \pi(f \circ \varphi_{tn})\|_j \le C e^{-atn} \|f\|_{j+1}. $$
As a consequence, the function $\Psi_t$ given by

$$ \Psi_t = \sum_{n=0}^{\infty} (P_t^{(0)})^n (\mu_0 f - f) \qquad (16) $$

is well-defined, in $\mathcal{A}$ and such that, for all $j \ge 0$, $\|\Psi_t\|_j \le C_j \|f\|_{j+1}$ for some $C_j > 0$ independent from $f$. Moreover,

$$ (P_t^{(0)} - \mathrm{Id})\Psi_t = \sum_{n=1}^{\infty} (P_t^{(0)})^n (\mu_0 f - f) - \sum_{n=0}^{\infty} (P_t^{(0)})^n (\mu_0 f - f) = f - \mu_0 f. $$

Then, using this Poisson equation and that $\mu_\varepsilon$ is invariant for $P_t^\varepsilon$,

$$ \mu_\varepsilon f - \mu_0 f = \mu_\varepsilon (P_t^{(0)} - \mathrm{Id})\Psi_t = \mu_\varepsilon (P_t^{(0)} - P_t^\varepsilon)\Psi_t. $$

To get the convergence of $\mu_\varepsilon f$ to $\mu_0 f$ at zeroth order, we simply bound, thanks to Proposition 1,

$$ |\mu_\varepsilon f - \mu_0 f| \le \|(P_t^{(0)} - P_t^\varepsilon)\Psi_t\|_\infty \le C\varepsilon \|\Psi_t\|_3. $$

To get the higher order expansion, write

$$ \begin{aligned} \mu_\varepsilon f - \mu_0 f &= \mu_\varepsilon (P_t^{(0)} - P_t^\varepsilon)\Psi_t \\ &= -\mu_\varepsilon \left( \sum_{k=1}^n \varepsilon^k P_t^{(k)} + \sum_{k=0}^n \varepsilon^k S_{\frac{t}{\varepsilon}}^{(k)} + R_{n,t}^\varepsilon \right) \Psi_t. \end{aligned} \qquad (17) $$

From Proposition 1 and the bounds on $\|\Psi_t\|_j$ for all $j \in \mathbb{N}$,

$$ \left| \mu_\varepsilon \left( \sum_{k=0}^{n} \varepsilon^k S_{\frac{t}{\varepsilon}}^{(k)} + R_{n,t}^\varepsilon \right) \Psi_t \right| \le C_n \varepsilon^{n+1} $$

for some $C_n$. Using that, from the convergence at order 0,

$$ |\mu_\varepsilon P_t^{(k)} \Psi_t - \mu_0 P_t^{(k)} \Psi_t| \le C_{0,k} \varepsilon $$

for all $k \in \mathbb{N}$ for some $C_{0,k} > 0$, we get first with (17) at order $n = 1$ that

$$ |\mu_\varepsilon f - \mu_0 f + \varepsilon \mu_0 P_t^{(1)} \Psi_t| \le C_{1,0} \varepsilon^2 \qquad (18) $$

for some $C_{1,0} > 0$. We can thus apply this first order expansion with $f$ replaced by $P_t^{(k)}\Psi_t$ to get

$$ |\mu_\varepsilon P_t^{(k)} \Psi_t - c_{0,k} - c_{1,k} \varepsilon| \le C_{1,k} \varepsilon^2 $$

for some $c_{0,k}, c_{1,k}, C_{1,k}$. Hence, considering the expansion (17) at order $n = 2$, we get

$$ |\mu_\varepsilon f - \mu_0 f - \varepsilon c_{0,1} - \varepsilon^2 c_{1,1} - \varepsilon^2 c_{0,2}| \le C_{2,0} \varepsilon^3, $$
for some $C_{2,0} > 0$. Again, this can be used with $f$ replaced by $P_t^{(k)}\Psi_t$, and a straightforward induction on $n$ concludes the proof of (13). In particular, from (18), we see that for all $t > 0$, $c_1$ can be written as

$$c_1 = -\mu_0 P_t^{(1)} \Psi_t.$$

**Step 2.** The goal is now to let $t$ vanish in this last expression. For a fixed $t > 0$, denote by $\mathbf{b}_1^\Psi$ the function given by (12) but with $f$ replaced by $\Psi_t$. Then

$$P_t^{(1)} \Psi_t(x, s) = \mathbf{b}_1^\Psi(t, x, s) + \alpha(x) + \beta(t, x)$$

with

$$\alpha(x) = \int_{0}^{+\infty} \pi(L_c S_r^{(0)} \Psi_t)(x) dr, \quad \beta(t, x) = \int_{0}^{t} \pi(L_c \mathbf{b}_1^\Psi)(r, \bar{\varphi}_{t-r}(x)) dr.$$

First, by definition of $Q^{-1}$, the average of $\mathbf{b}_1^\Psi$ with respect to $\pi$, hence with respect to $\mu_0$, is zero. At this stage, we have obtained that $c_1 = -\mu_0\alpha - \mu_0\beta$. Recalling that $S_r^{(0)}f = Q_r(f - \pi f)$, from

$$\Psi_t(x, s) = \pi f(\bar{x}) - f(x, s) + \sum_{n=1}^{\infty} \pi(f(\bar{x}) - f(\varphi_{tn}(x))),$$

integrating with respect to $\pi$, we get

$$\nabla_x S_r^{(0)} \Psi_t(\bar{x}, s) = \nabla_x Q_r(\Psi_t - \pi \Psi_t)(\bar{x}, s) = \nabla_x Q_r(\pi f - f)(\bar{x}, s).$$

As a consequence,

$$\begin{aligned}
\mu_0(\alpha) &= \alpha(\bar{x}) = \int_0^\infty \pi L_c S_r^{(0)} \Psi_t(\bar{x}) dr \\
&= \int_S F_s(\bar{x}) \cdot \left[ \int_0^\infty \nabla_x Q_r(\pi f - f) dr \right] (\bar{x}, s) \pi(ds) \\
&= \int_S F_s(\bar{x}) \cdot \nabla_x Q^{-1} f(\bar{x}, s) \pi(ds) \\
&= \pi L_c Q^{-1} f(\bar{x}).
\end{aligned}$$

Besides,

$$\begin{aligned}
\mu_0(\beta) &= \int_0^t \pi L_c \mathbf{b}_1^\Psi(r, \bar{x}) dr \\
&= \pi \left[ L_c \int_0^t \mathbf{b}_1^\Psi(r, \cdot) dr \right] (\bar{x}) \\
&= \pi \left[ L_c \int_0^t Q^{-1} \left( \partial_r P_r^{(0)} \Psi_t - L_c P_r^{(0)} \Psi_t \right) dr \right] (\bar{x}) \\
&= -\pi \left[ L_c Q^{-1} \int_0^t L_c P_r^{(0)} \Psi_t dr \right] (\bar{x}),
\end{aligned}$$
where we used that $\hat{\partial}_r P_r^{(0)} \Psi_t$ does not depend on $s$ and is thus in the kernel of $Q^{-1}$. Hence,

$$ \int_0^t \pi L_c \mathbf{b}_1^\Psi(s, \bar{x}) ds = -\pi \left[ L_c Q^{-1} \frac{1}{t} \int_0^t L_c P_r^{(0)}(t \Psi_t) dr \right] (\bar{x}) \\ \xrightarrow{t \to 0} -\pi [L_c Q^{-1} L_c(\pi h)](\bar{x}), $$

where we used that $P_0^{(0)} f = \pi f$ and that $t\Psi_t \to h$. Indeed, thanks to Assumption 2,

$$ t\Psi_t(x) = t\pi f(\bar{x}) - tf(x, s) + t \sum_{n=1}^{\infty} \pi(f(\bar{x}) - f(\varphi_{tn}(x))) \\ \xrightarrow{t \to 0} \int_0^{\infty} \pi(f(\bar{x}) - f(\varphi_r(x))) dr = h(x), $$

and the convergence $t\Psi \to h$ holds in all norms $\|\cdot\|_j$. The proof of (14) is concluded by letting $t$ vanish in the equality $c_1 = -\mu_0\alpha - \mu_0\beta$ (since $c_1$ is independent from $t$). □

# 3 Applications

## 3.1 Top Lyapunov exponent for cooperative linear Markov-modulated ODEs

In this section, we consider linear Markov-modulated ODEs on $\mathbb{R}^d$ of the form

$$ z'(t) = A(\sigma(t/\epsilon)) z(t) \qquad (19) $$

where, for all $\sigma \in S$, $A(\sigma)$ is a $d \times d$ matrix. We work under Assumption 1 and write $\bar{A} = \int_S A(s)\pi(ds)$. Moreover, we focus on the settings of [7], given as follows:

**Assumption 3.** 1. The Markov process $(\sigma(t))_{t \ge 0}$ is Feller.

2. For all $s \in S$, $A(s)$ is a cooperative matrix, in the sense that its off-diagonal coefficients are non-negative.

3. The averaged matrix $\bar{A}$ is irreducible in the sense that for all $i, j \in [1, d]$ there exists a path $i = i_0, \dots, i_q = j$ with $\bar{A}_{i_{k-1}, i_k} > 0$ for all $k \in [1, q]$.

The fact that the matrices are cooperative implies that $\mathbb{R}_+^d$ is fixed by (19) for $t \ge 0$ and that Assumption 2 holds. We decompose solutions of (19) on $\mathbb{R}_+^d$ as $z(t) = \rho(t)\theta(t)$ where $\rho(t) = 1 \cdot z(t) > 0$ with $1 = (1, \dots, 1) \in \mathbb{R}^d$ and $\theta(t) = z(t)/\rho(t) \in \Delta = \{x \in \mathbb{R}_+^d, x_1 + \dots + x_d = 1\}$. The ODE (19) is then equivalent to

$$ \rho'(t) = (1 \cdot A(\sigma(t/\epsilon)) \theta(t)) \rho(t), \quad \theta'(t) = F_{\sigma(t/\epsilon)}(\theta(t)) $$
with

$$F_s(\theta) = A(s)\theta - (1 \cdot A(s)\theta)\theta.$$

It is proven in [7] that, under Assumption 3, the Markov process $(\theta(t), \sigma(t/\epsilon))_{t \ge 0}$ admits a unique invariant measure $\mu_\epsilon$ on $\Delta \times S$ and that for all initial conditions $z \in \mathbb{R}_+^d \setminus \{0\}$, almost surely,

$$\lim_{t \to \infty} \frac{\rho(t)}{t} = \Lambda_{\epsilon}, \qquad (20)$$

where

$$\Lambda_{\epsilon} := \int_{\Delta \times S} 1 \cdot A(s) \theta \, \mu_{\epsilon}(d\theta, ds) \qquad (21)$$

is called the top Lyapunov exponent of the process. Denoting by $\lambda_{\max}(\bar{A})$ the principal eigenvalue of $\bar{A}$ (i.e., with maximal real part), Proposition 4 in [7] entails that

$$\lim_{\epsilon \to 0} \Lambda_{\epsilon} = \lambda_{\max}(\bar{A}).$$

From Theorem 1 applied to the function $f : (\theta, s) \mapsto 1 \cdot A(s)\theta$, we get an expansion of $\Lambda_\epsilon$ for small $\epsilon$.

**Proposition 2.** Under Assumptions 1 and 3, there exists a sequence $(c_k)_{k \ge 1}$ of real numbers such that, for all $n \ge 1$,

$$\Lambda_{\epsilon} = \lambda_{\max}(\bar{A}) + \sum_{k=1}^{n} c_k \epsilon^k + \mathcal{O}_{\epsilon \to 0}(\epsilon^n).$$

Moreover, denoting by $\bar{x}$ and $\bar{y}$, respectively, the right and left eigenvectors of $\bar{A}$ associated with $\lambda_{\max}(\bar{A})$ and such that $\bar{x} \in \Delta$ and $\bar{x} \cdot \bar{y} = 1$, it holds

$$c_1 = \int_S \bar{y}^\top Q^{-1}(A)(s) (\bar{x}\bar{y}^\top - I) A(s) \bar{x} \pi(ds).$$

*Proof.* First, considering $h$ given by (15), we claim that for all $(x, s) \in \Delta \times S$,

$$h(x, s) = -\ln(\bar{y} \cdot x). \qquad (22)$$

Indeed, note that, for all $x \in \mathcal{M}$,

$$\begin{aligned} \pi f(\varphi_r(x)) &= \frac{1 \cdot \bar{A} e^{r \bar{A}} x}{1_d \cdot e^{r \bar{A}} x} \\ &= \frac{d}{dr} \ln(1 \cdot e^{r \bar{A}} x). \end{aligned}$$

Theorefore, using that $\bar{x} = \varphi_r(\bar{x})$, one has

$$\begin{aligned} \int_0^\infty (\pi f(\bar{x}) - \pi f(\varphi_r(x))) dr &= \int_0^\infty \frac{d}{dr} \left( \ln \left( \frac{1 \cdot e^{r \bar{A}} \bar{x}}{1 \cdot e^{r \bar{A}} x} \right) \right) dr \\ &= \lim_{r \to \infty} \ln \left( \frac{1 \cdot e^{r \bar{A}} \bar{x}}{1 \cdot e^{r \bar{A}} x} \right) \\ &= -\ln(\bar{y} \cdot x), \end{aligned}$$
where we have used Perron-Frobenius Theorem. From (22), we deduce that
$L_c h(x, s) = -F_s(x) \cdot \frac{\bar{y}}{\bar{y} \cdot x}$. Now, using that $F_s(x) = A(s)x - f_s(x)x$, we get that

$$
L_c h(x, s) - f(x, s) = - \frac{A(s)x \cdot \bar{y}}{\bar{y} \cdot x},
$$

and then

$$
Q^{-1}(L_c h - f)(x, s) = - \frac{Q^{-1}(A)(s)x \cdot \bar{y}}{\bar{y} \cdot x},
$$

from which we deduce that

$$
\begin{align*}
\nabla_x Q^{-1}(L_c h - f)(\bar{x}, s) &= \frac{Q^{-1}(A)(s)\bar{x} \cdot \bar{y}}{(\bar{y} \cdot \bar{x})^2} \bar{y} - \frac{[Q^{-1}(A)(s)]^\top \bar{y}}{\bar{y} \cdot \bar{x}} \\
&= (Q^{-1}(A)(s)\bar{x} \cdot \bar{y})\bar{y} - [Q^{-1}(A)(s)]^\top \bar{y},
\end{align*}
$$

where we have used that $\bar{x} \cdot \bar{y} = 1$ As a consequence,

$$
\begin{align*}
F_s(\bar{x}) \cdot \nabla_x Q^{-1}(L_c h - f)(\bar{x}, s) & \\
&= (A(s)\bar{x} - (\mathbf{1} \cdot A(s)\bar{x})\bar{x}) \cdot ((Q^{-1}(A)(s)\bar{x} \cdot \bar{y})\bar{y} - [Q^{-1}(A)(s)]^\top \bar{y}) \\
&= (Q^{-1}(A)(s)\bar{x} \cdot \bar{y})A(s)\bar{x} \cdot \bar{y} - A(s)\bar{x} \cdot [Q^{-1}(A)(s)]^\top \bar{y} \\
&\quad - (\mathbf{1} \cdot A(s)\bar{x})(Q^{-1}(A)(s)\bar{x} \cdot \bar{y})\bar{x} \cdot \bar{y} + (\mathbf{1} \cdot A(s)\bar{x})\bar{x} \cdot [Q^{-1}(A)(s)]^\top \bar{y}.
\end{align*}
$$

Using that $\bar{x} \cdot \bar{y} = 1$, we see that the last line is equal to

$$
-(\mathbf{1} \cdot A(s)\bar{x})Q^{-1}(A)(s)\bar{x} \cdot \bar{y} + (\mathbf{1} \cdot A(s)\bar{x})Q^{-1}(A)(s)\bar{x} \cdot \bar{y} = 0.
$$

We end up with

$$
\begin{align*}
F_s(\bar{x}) \cdot \nabla_x Q^{-1}(L_c h - f)(\bar{x}, s) & \\
&= (Q^{-1}(A)(s)\bar{x} \cdot \bar{y})A(s)\bar{x} \cdot \bar{y} - Q^{-1}(A)(s)A(s)\bar{x} \cdot \bar{y} \\
&= \bar{y}^\top Q^{-1}(A)(s) [\bar{x}\bar{y}^\top - I] A(s)\bar{x}.
\end{align*}
$$

Integrating with respect to $\pi$ leads to

$$
c_1 = \int_S \bar{y}^\top Q^{-1}(A)(s) (\bar{x} \bar{y}^\top - I) A(s) \bar{x} \pi(ds).
$$

<figure><img src="image_3.png" /></figure>

**Example 3.** If $\sigma$ is a Markov chain on $S = \{0,1\}$, with matrices rates given by

$$
Q = \begin{pmatrix} -p & p \\ 1-p & -(1-p) \end{pmatrix}, \qquad (23)
$$

we get the following simple expression for $c_1$:

$$
c_1 = p(1-p) \left[ \bar{y}^{\top} (A_0 - A_1)^2 \bar{x} - (\bar{y}^{\top} (A_0 - A_1) \bar{x})^2 \right]. \quad (24)
$$
Indeed, recalling that $Q^{-1} = Q$ (see Example 2) and using that $\pi_s Q_{s,s'} = p(1-p)$ if $s \neq s'$ and $-p(1-p)$ if $s = s'$, we get

$$c_1 = \bar{y}^T (A_0 P A_0 + A_1 P A_1 - A_0 P A_1 - A_1 P A_0) \bar{x},$$

where $P = I - \bar{x}\bar{y}^T$. Now,

$$\begin{aligned} & \bar{y}^T A_0 \bar{x} \bar{y}^T A_1 \bar{x} + \bar{y}^T A_1 \bar{x} \bar{y}^T A_0 \bar{x} \\ & \qquad - \bar{y}^T A_0 \bar{x} \bar{y}^T A_0 \bar{x} - \bar{y}^T A_1 \bar{x} \bar{y}^T A_1 \bar{x} \\ & = -(\bar{y}^T (A_0 - A_1) \bar{x})^2 \end{aligned}$$

and

$$\bar{y}^T A_0 A_1 \bar{x} + \bar{y}^T A_1 A_0 \bar{x} - \bar{y}^T A_0 A_0 \bar{x} - \bar{y}^T A_1 A_1 \bar{x} = - \bar{y}^T (A_0 - A_1)^2 \bar{x},$$

which induces formula (24).

When considering the switching between matrices $(A_i)_{i \in S}$, the following natural question arises: if all the matrices are stable (i.e. $\lambda_{\max}(A_i) < 0$ for all $i \in S$), is the switched system (19) stable, in the sense that $\Lambda_{\epsilon} < 0$? It is now known that this is not true. In the case of random switching between two $2 \times 2$ matrices, examples of stable matrices giving an unstable system can be founded in [5] and [20]. In the first reference, for $p = 1/2$, $\lambda_{\max}(\bar{A}) < 0$ so that fast switching leads to an unstable system. In the second reference, it is a bit more complicated, since $\lambda_{\max}(\bar{A}) > 0$ for all $p \in (0, 1)$, meaning that both slow and fast switching lead to a stable system. However, the authors prove in [20] that switching not too fast, nor too slowly, the system is unstable. Yet, in [20], the matrices are not cooperative in the sense of Assumption 3. This is not surprising, since it is proven in [16] that switching between cooperative matrices of size $2 \times 2$ such that every matrices in the convex hull of the given matrices are stable, will always lead to a stable system. However, it is also shown in [16] that it is possible in some higher dimension to construct an example where all the matrices in the convex hull are stable, and for which there exists a periodic switching such that the linear system explodes. Later, an explicit example in dimension 3 was given by Fainshil, Margaliot and Chiganski [14]. Precisely, consider the matrices

$$A_0 = \begin{pmatrix} -1 & 0 & 0 \\ 10 & -1 & 0 \\ 0 & 0 & -10 \end{pmatrix}, \quad A_1 = \begin{pmatrix} -10 & 0 & 10 \\ 0 & -10 & 0 \\ 0 & 10 & -1 \end{pmatrix}. \qquad (25)$$

It is shown in [14] that every convex combination of $A_0$ and $A_1$ is stable, and yet a switch of period 1 between $A_0$ and $A_1$ yields an explosion. In [8], the authors asked the question whether the same system but with a Markovian switching can lead to explosion. Using numerical simulations, they suggest that this is true and the Lyapunov exponent is positive for not too fast nor too slow switching. Now, using Formula (24), we prove rigorously the following assertion:

**Proposition 3.** There exists two cooperative matrices $B_0, B_1$ such that:
<figure><img src="image_4.png" /><figcaption>(a) $p \mapsto \lambda_{\max}(p)$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) $p \mapsto c_1(p)$</figcaption></figure>

1. For all $p \in [0, 1]$, $\lambda_{\max}((1-p)B_0 + pB_1) < 0$,

2. For some $\varepsilon > 0$ and $p \in (0, 1)$, the top Lyapunov exponent $\Lambda_{\varepsilon}$ of the system (19) with $\sigma$ generated by Q given in (23) is positive.

*Proof.* To emphasize the dependency on $p$, we write $\bar{A}(p)$, $\Lambda_{\varepsilon}(p)$ and $c_1(p)$ for the mean matrix, the top Lyapunov exponent and the first order derivative, respectively, when the switching rates are given by the matrix $Q$ in (23), and the matrices are (25). We also write $\lambda_{\max}(p)$ for $\lambda_{\max}(\bar{A}(p))$. In Figure 1a is plotted $p \mapsto \lambda_{\max}(p)$. One sees that the maximal value of this function is attained at some $p^* \in [0.3, 0.5]$ and worth $\lambda_{\max}^* \in [-0.5, -0.45]$. One sees on Figure 1b, where is plotted $p \mapsto c_1(p)$, that on the interval $[0.3, 0.5]$, we have $c_1(p) \ge 15$, so that in particular $c_1(p^*) \ge 15 > 0$.

Let $\delta > 0$ and $A_i^\delta = A_i + \delta I$. Replacing $A_i$ by $A_i^\delta$ in (19), we easily check that the Lyapunov exponent of the system is $\Lambda_{\varepsilon}^\delta(p) = \Lambda_{\varepsilon}(p) + \delta$. Moreover, $f_i^\delta(x) = 1 \cdot (A_i + \delta I)x$, so that $F_i^\delta = F_i$, and the dynamic of the angular part of (19) is independent of $\delta$. Hence, Proposition 2 entails

$$ \Lambda_{\varepsilon}^{\delta}(p) \ge \lambda_{\max}(p) + \delta + c_1(p)\varepsilon - M(p)\varepsilon^2, $$

for some constant $M(p)$ depending on $p$. Choosing $p = p^*$ and $\varepsilon = \frac{c_1(p^*)}{2M(p^*)} > 0$ yields

$$ \Lambda_{\varepsilon}^{\delta}(p^*) \ge \lambda_{\max}(p^*) + \delta + \frac{c_1(p^*)^2}{4M(p^*)}. $$

Letting $\delta = -\frac{c_1(p^*)^2}{8M(p^*)} - \lambda_{\max}(p^*)$ concludes the proof, since

1. $\Lambda_{\varepsilon}^{\delta}(p^*) \ge \frac{c_1^2(p^*)}{8M(p^*)} > 0$;

2. For all $p \in [0, 1]$, $\lambda_{\max}(A_p^\delta) \le \lambda_{\max}(A_{p^*}^\delta) = \lambda_{\max}(p^*) + \delta = -\frac{c_1^2(p^*)}{8M(p^*)} < 0$. $\square$
**Example 4.** Assume as in Example 1 that $Qf = Q^{-1}f = \pi f - f$. Then, $Q^{-1}(A)(s) = \bar{A} - A(s)$. Using $\bar{y}^\top \bar{A} = \lambda(\bar{A}) \bar{y}^\top$ and $\bar{y}^\top \bar{x} = 1$, we get that $\bar{y}^\top \bar{A}(\bar{x} \bar{y}^\top - I) = 0$, which yields

$$c_1 = \int_S \bar{y}^\top A^2(s) \bar{x} - (\bar{y}^\top A(s) \bar{x})^2 \pi(ds).$$

Notice that, from Jensen's inequality,

$$\int_S (\bar{y}^\top A(s) \bar{x})^2 \pi(ds) \geq \left( \int_S \bar{y}^\top A(s) \bar{x} \pi(ds) \right)^2 = \lambda_{\max}(\bar{A})^2,$$

so that

$$c_1 \leq \int_S \bar{y}^\top A^2(s) \bar{x} \pi(ds) - \lambda_{\max}(\bar{A})^2 = \bar{y}^\top \overline{A^2} \bar{x} - \lambda_{\max}(\overline{A^2}) = \bar{y}^\top (\overline{A^2} - \overline{A^2}) \bar{x}.$$

## 3.2 Richardson extrapolation for randomized splitting schemes

In [1], in order to approximate the solution of an ODE of the form

$$\dot{y}(t) = \bar{F}(y(t)), \qquad (26)$$

the vector field $\bar{F}$ is split as $\bar{F}(x) = \sum_{k=1}^N F_k(x)$, where each flow $\dot{x} = F_k(x)$ can be solved exactly. More precisely, the two main examples of interest in [1] (to which we refer for details and motivations) are the so-called Lorenz-96 model and a finite-dimensional Galerkin projection of the vorticity formulation of 2D Navier-Stokes. The trajectory $y$ is then approximated by a discrete time scheme of the form

$$x_{n+1}^{\varepsilon} = \Phi_{\tau_1/\varepsilon}^{1} \circ \cdots \circ \Phi_{\tau_N/\varepsilon}^{N}(x_n^{\varepsilon})$$

where $\Phi_t^i$ is the flow associated to $F_i$ at time $t$ and $\tau_1, \dots, \tau_N$ are independent random variables distributed according to the standard exponential law. More precisely, for a fixed $t > 0$, $y(t)$ is approximated by $x_n^\varepsilon$ with $n\varepsilon = t$, for a small $\varepsilon$. Notice that this does not enter directly our framework. However, consider the Markov-modulated ODE

$$\dot{x}_{\varepsilon}(t) = F_{\sigma(Nt/\varepsilon)}(x_{\varepsilon}(t)),$$

where $\sigma$ is a cyclic chain on $[1, N]$ with rate 1 (i.e. it is a standard Poisson Process modulo $N$). We can take the times $\tau_i$ to define the jump times of $\sigma$, and thus $x_n^\varepsilon$ is exactly $x_\varepsilon(\varepsilon S_n/N)$ where $S_n$ is the $(nN)^{th}$ jump time of $\sigma$. In particular, if $n\varepsilon = t$, assuming that the vector fields $F_k$ are bounded,

$$|x_{\varepsilon}(t) - x_{n}^{\varepsilon}| = |x_{\varepsilon}(t) - x_{\varepsilon}(\varepsilon S_{n}/N)| \leq \max_{k \in [1, N]} \|F_{k}\|_{\infty} |\varepsilon S_{n}/N - t|,$$

and

$$\mathbb{E} \left| \frac{\varepsilon S_n}{N} - t \right|^2 = \frac{n\varepsilon^2}{N} = \frac{t\varepsilon}{N}.$$
In other words, $(x_{\varepsilon}(t))_{t \ge 0}$ (which can be sampled exactly if $(x_n^\varepsilon)_{n \in \mathbb{N}}$ can) can be seen as an alternative variation of the scheme of [1].

The convergence of $x_n^\varepsilon$ to $y(t)$ as $\varepsilon \to 0$ for a fixed $t = n\varepsilon$, or more precisely of the corresponding Markov transition operators, is stated in [1, Theorem 4.1], which is thus similar to our Proposition 1 at order $n = 0$. Notice that, a priori, the ODE (26) is in $\mathbb{R}^d$ but [1, Theorem 4.1] is proven by restricting the study on the orbit of the process, which is assumed to be bounded (this is Assumption 1 of [1]), so that our results apply in this context.

The higher order expansion in Proposition 1 enables the use of a Richardson extrapolation to get better convergence rates: for instance we obtain, for a fixed $t > 0$ and $f \in \mathcal{A}$,

$$|f(\varphi_t(x)) - 2\mathbb{E}(f(x_{\varepsilon/2}(t))) + \mathbb{E}(f(x_{\varepsilon}(t)))| \le C \|f\|_3 \varepsilon^2,$$

for some constant $C > 0$.

## 3.3 Species invasion rate in a two species Lotka-Volterra model in random environment

We consider a model of two species in interaction, in a random environment, described by the following system of Lotka - Volterra equations:

$$\left\{ \begin{aligned} \frac{dX(t)}{dt} &= X(t) \left( a_{10}^{\sigma(t)} - a_{11}^{\sigma(t)} X(t) + a_{12}^{\sigma(t)} Y(t) \right) \\ \frac{dY(t)}{dt} &= Y(t) \left( a_{20}^{\sigma(t)} + a_{21}^{\sigma(t)} X(t) - a_{22}^{\sigma(t)} Y(t) \right) \end{aligned} \right. \quad (27)$$

We assume that $a_{11}^s, a_{22}^s > 0$ (intraspecific competition), the other parameters can be positive or negative. The coefficients $a_{10}^s$ and $a_{20}^s$ represents the growth rates of species 1 and 2 respectively, while the coefficients $a_{12}^s$ and $a_{21}^s$ models the interaction between the two species. If for example $a_{10}^s, a_{21}^s > 0$ while $a_{20}^s, a_{12}^s < 0$, we get a prey - predator system (species 1 being the prey, and species 2 the predator). Assumption 1 is also enforced in this whole section.

In the absence of species 2 and after accelerating the dynamics of $\sigma$ by a factor $1/\varepsilon$, species 1 evolves according to a logistic equation in random environment:

$$\frac{dX(t)}{dt} = X(t) \left( a_{10}^{\sigma(t/\varepsilon)} - a_{11}^{\sigma(t/\varepsilon)} X(t) \right). \quad (28)$$

We prove in the Annex the following result:

**Lemma 1.** Assume that $p_0 := \inf_{s \in S} \frac{a_{10}^s}{a_{11}^s} > 0$, $p_1 := \sup_{s \in S} \frac{a_{10}^s}{a_{11}^s} < +\infty$ and $\inf_{s \in S} a_s^{11} > 0$. Then the process $(X(t), \sigma(t/\varepsilon))_{t \ge 0}$ with $X$ given by (28) and $X(0) \ne 0$ will ultimately lie in $E := [p_0, p_1] \times S$ and admits a unique stationary distribution $\mu_\varepsilon$ on $E$. Moreover, $(X(t), \sigma(t/\varepsilon))$ converges in law to $\mu_\varepsilon$ when $t$ goes to infinity.

The stationary distribution $\mu_\varepsilon$ represents the population of species 1 at equilibrium in the absence of species 2. The invasion growth rate of species 2 when
species 1 is at equilibrium is then given by

$$ \Lambda_y(\varepsilon) = \int_{M \times S} (a_{20}^s + a_{21}^s x) \mu_\varepsilon(dx, ds). \qquad (29) $$

Intuitively, we are interested in the sign of this quantity, since when the abundance of the second species is close to 0, we have roughly

$$ \frac{1}{t} \ln(Y(t)) \approx \frac{1}{t} \int_{0}^{t} \left[ a_{20}^{\sigma(u/\varepsilon)} + a_{21}^{\sigma(u/\varepsilon)} X(u) \right] du, $$

where $X$ solves (28) and the right hand side converges to $\Lambda_y(\varepsilon)$ when $t$ goes to infinity, thanks to Lemma 1 and the ergodic theorem. It is made rigorous in [4] that, indeed, the sign of $\Lambda_y(\varepsilon)$ gives information on the local behaviour of (27) near the boundary $\{y=0\}$. Denoting by $\bar{\alpha} = \int_S \alpha^s \pi(ds)$ the mean with respect to $\pi$ of a function $\alpha$ defined on $S$, it is proven in [6], in the specific case when $\sigma$ is a two-states Markov chain, that $\Lambda_y(\varepsilon)$ converges as $\varepsilon$ goes to 0 to $\bar{a}_{20} + \bar{a}_{21}\bar{x}$. The proposition below extends this result to the general case of a process $\sigma$ satisfying Assumption 1 and with a first order expansion:

**Proposition 4.** Assume that $\inf_{s \in S} a_{10}^s > 0$ and $\sup_{s \in S} a_{11}^s < +\infty$. Then, for any invariant probability measure $\mu_\varepsilon$ of $(X_t, \sigma(t/\varepsilon))_{t \ge 0}$,

$$ \Lambda_y(\varepsilon) = \bar{a}_{20} + \bar{a}_{21}\bar{x} + \varepsilon \frac{\bar{a}_{21}\bar{a}_{10}^2}{\bar{a}_{11}} \int_S \left( \frac{a_{10}^s}{\bar{a}_{10}} - \frac{a_{11}^s}{\bar{a}_{11}} \right) Q^{-1} \left( \frac{a_{11}}{\bar{a}_{11}} - \frac{a_{21}}{\bar{a}_{21}} \right) (s) \pi(ds) + o_{\varepsilon \to 0}(\varepsilon). $$

*Proof.* First, note that the unique equilibrium of $\bar{F}(x) = x(\bar{a}_{10} - \bar{a}_{11}x)$ is $\bar{x} := \bar{a}_{10}/\bar{a}_{11}$ and Assumption 2 holds (since thanks to Lemma 1, it is sufficient to consider initial conditions in the compact set $[p_0, p_1]$). By Theorem 1, and since the marginal of $\mu_\varepsilon$ in $s$ is independent from $\varepsilon$ (so that $\mu_\varepsilon(a_{20}) = \bar{a}_{20}$), one has

$$ \Lambda_y(\varepsilon) = \bar{a}_{20} + \bar{a}_{21}\bar{x} + \varepsilon c_1 + o_{\varepsilon \to 0}(\varepsilon), \qquad (30) $$

with

$$ c_1 = \pi L_c Q^{-1} (L_c h - f)(\bar{x}), $$

for $f(x, s) = a_{21}^s x$. Using the formula for $h$ derived in Remark 1, we have

$$ \begin{aligned} h(x, s) &= \int_{\bar{x}}^{x} \frac{\bar{f}(y) - \bar{f}(\bar{x})}{\bar{F}(y)} dy \\ &= \int_{\bar{x}}^{x} \bar{a}_{21} \frac{y - \bar{x}}{y(\bar{a}_{10} - \bar{a}_{11}y)} dy \\ &= \frac{\bar{a}_{21}}{\bar{a}_{11}} \ln \left( \frac{\bar{x}}{x} \right), \end{aligned} $$

which yields

$$ L_c h(x, s) = -(a_{10}^s - a_{11}^s x) \frac{\bar{a}_{21}}{\bar{a}_{11}}. $$
Then, $\nabla_x(L_c h - f) = a_{11}\bar{a}_{21}/\bar{a}_{11} - a_{21}$, so that

$$L_c Q^{-1} (L_c h - f) (\bar{x}, s) = \bar{x}(a_{10}^s - a_{11}^s \bar{x}) Q^{-1} \left( \frac{\bar{a}_{21}}{\bar{a}_{11}} a_{11} - a_{21} \right) (s)$$

and finally, integrating this with respect to $\pi$ and using that $\bar{x} = \bar{a}_{10}/\bar{a}_{11}$,

$$c_1 = \frac{\bar{a}_{21}\bar{a}_{10}^2}{\bar{a}_{11}} \int_S \left( \frac{a_{10}^s}{\bar{a}_{10}} - \frac{a_{11}^s}{\bar{a}_{11}} \right) Q^{-1} \left( \frac{a_{11}}{\bar{a}_{11}} - \frac{a_{21}}{\bar{a}_{21}} \right) (s) \pi(ds), \quad (31)$$

which concludes.

<figure><img src="image_6.png" /></figure>

Let us discuss the sign of the expression (31) in some particular cases.

1. If $a_{11}^s = \bar{a}_{11}$ and $a_{10}^s = \bar{a}_{10}$ for all $s$, then $\Lambda'_y(0) = 0$. This is expected, since in that case, the process $X$ solution to (28) is deterministic and does not depend on $\varepsilon$, and $\mu_\varepsilon = \delta_{\bar{a}_{10}/\bar{a}_{11}} \otimes \pi$ for all $\varepsilon$. More interestingly, if $a_{11}^s = \bar{a}_{11}$ and $a_{21}^s = \bar{a}_{21}$ but $a_{10}$ is varying, we still have $\Lambda'_y(0) = 0$.

2. If $a_{10}^s = \bar{a}_{10}$ and $a_{21}^s = \bar{a}_{21}$ for all $s$, using that $\pi Q^{-1} = Q^{-1}1 = 0$, we get

$$\begin{aligned} c_1 &= -\frac{\bar{a}_{21}\bar{a}_{10}^2}{\bar{a}_{11}} \int_S \frac{a_{11}^s}{\bar{a}_{11}} Q^{-1} \left( \frac{a_{11}}{\bar{a}_{11}} \right) (s) \pi(ds) \\ &= -\frac{\bar{a}_{21}\bar{a}_{10}^2}{\bar{a}_{11}^3} \int_S (a_{11}^s - \bar{a}_{11}) Q^{-1} (a_{11} - \bar{a}_{11}) (s) \pi(ds). \end{aligned}$$

This is always non-negative, as it can be interpreted as an asymptotic variance. Indeed, considering $g \in \mathcal{A}$ with $\pi g = 0$, for $t \ge 0$,

$$\begin{aligned} t \mathbb{E}_{\pi} \left| \frac{1}{t} \int_{0}^{t} g(\sigma(s)) ds \right|^{2} &= \frac{2}{t} \mathbb{E}_{\pi} \int_{0}^{t} \int_{s}^{t} g(\sigma(u)) g(\sigma(s)) ds du \\ &= \frac{2}{t} \int_{0}^{t} \int_{s}^{t} \pi(g Q_{u-s} g) ds du \\ &= \frac{2}{t} \int_{0}^{t} \int_{0}^{t} \mathbb{1}_{v<t-s} \pi(g Q_v g) ds dv \\ &= \frac{2}{t} \int_{0}^{t} (t-v) \pi(g Q_v g) dv \\ &\xrightarrow{t \to \infty} -2\pi (g Q^{-1} g), \end{aligned}$$

where we used that $\int_0^t Q_v g$ converges to $-Q^{-1}g$ and bounded

$$\left| \int_{0}^{t} v \pi(g Q_v g) dv \right| \le \int_{0}^{\infty} C v e^{-\gamma v} \|g\|_{\infty}^{2} dv.$$
**Example 5.** We consider the two-state case, i.e., $\sigma$ is Markov chain on $S = \{0, 1\}$, with matrices rates given by

$$Q = \begin{pmatrix} -p & p \\ 1-p & -(1-p) \end{pmatrix}$$

for some $p \in (0, 1)$. The behaviour of the solutions to (27) was studied by Benaïm and Lobry in [6], through the signs of the invasion growth rates of species 1 and 2, in the competitive case (ie, when $a_{21}^s$ and $a_{12}^s$ are negative). Their study was complemented by Malrieu and Zitt in [21], who proposed an alternative formula for the invasion growth rate which made possible for them to understand the monotonicity of $\varepsilon \mapsto \Lambda_y(\varepsilon)$. In particular, it is a consequence of the proof of Lemma 4.1 in [21] that $\varepsilon \mapsto \Lambda_y(\varepsilon)$ is increasing if $A_2 > 0$ while it is decreasing if $A_2 < 0$, where $A_2$ is the coefficient of the second order term of the polynomial

$$P(x) = \left[ \frac{a_{20}^1}{a_{10}^1} \left(1 + \frac{a_{21}^1}{a_{20}^1} x\right) \left(1 - \frac{a_{10}^0}{a_{10}^0} x\right) - \frac{a_{20}^0}{a_{10}^0} \left(1 + \frac{a_{21}^0}{a_{20}^0} x\right) \left(1 - \frac{a_{11}^1}{a_{10}^1} x\right) \right] \frac{\frac{a_{11}^1}{a_{10}^1} - \frac{a_{11}^0}{a_{10}^0}}{|\frac{a_{11}^1}{a_{10}^1} - \frac{a_{11}^0}{a_{10}^0}|},$$

that is $A_2 = (a_{11}^1 a_{21}^0 - a_{11}^0 a_{21}^1)/(a_{10}^0 a_{10}^1)$, if we assume without loss of generality that $\frac{a_{11}^1}{a_{10}^1} \ge \frac{a_{11}^0}{a_{10}^0}$. Let us prove that our results are in accordance with the results of Malrieu and Zitt, by studying the sign of the first order term of the expansion of $\varepsilon \mapsto \Lambda_y(\varepsilon)$ when $\varepsilon \to 0$. Using Example 2, we have $\pi = (1-p, p)$, $Q^{-1} = Q$ and $\pi_s Q_{s,s'} = p(1-p)$ if $s \ne s'$ and $-p(1-p)$ if $s = s'$, from which we compute

$$\begin{aligned}
c_1 &= \sum_{s,s' \in \{0,1\}} \pi_s Q_{s,s'} \left( \frac{a_{10}^s}{\bar{a}_{10}} - \frac{a_{11}^s}{\bar{a}_{11}} \right) \left( \frac{a_{11}^{s'}}{\bar{a}_{11}} - \frac{a_{21}^{s'}}{\bar{a}_{21}} \right) \\
&= p(1-p) \left[ \left( \frac{a_{10}^0}{\bar{a}_{10}} - \frac{a_{11}^0}{\bar{a}_{11}} \right) \left( \frac{a_{11}^1}{\bar{a}_{11}} - \frac{a_{21}^1}{\bar{a}_{21}} \right) + \left( \frac{a_{10}^1}{\bar{a}_{10}} - \frac{a_{11}^1}{\bar{a}_{11}} \right) \left( \frac{a_{11}^0}{\bar{a}_{11}} - \frac{a_{21}^0}{\bar{a}_{21}} \right) \right. \\
&\qquad \left. - \left( \frac{a_{10}^0}{\bar{a}_{10}} - \frac{a_{11}^0}{\bar{a}_{11}} \right) \left( \frac{a_{11}^0}{\bar{a}_{11}} - \frac{a_{21}^0}{\bar{a}_{21}} \right) - \left( \frac{a_{10}^1}{\bar{a}_{10}} - \frac{a_{11}^1}{\bar{a}_{11}} \right) \left( \frac{a_{11}^1}{\bar{a}_{11}} - \frac{a_{21}^1}{\bar{a}_{21}} \right) \right] \\
&= p(1-p) \left( \frac{a_{10}^0}{\bar{a}_{10}} - \frac{a_{11}^0}{\bar{a}_{11}} - \left( \frac{a_{10}^1}{\bar{a}_{10}} - \frac{a_{11}^1}{\bar{a}_{11}} \right) \right) \left( \frac{a_{11}^1}{\bar{a}_{11}} - \frac{a_{21}^1}{\bar{a}_{21}} - \left( \frac{a_{11}^0}{\bar{a}_{11}} - \frac{a_{21}^0}{\bar{a}_{21}} \right) \right) \\
&= p(1-p) \left( \frac{a_{10}^0 - a_{10}^1}{\bar{a}_{10}} - \frac{a_{11}^0 - a_{11}^1}{\bar{a}_{11}} \right) \left( \frac{a_{11}^1 - a_{11}^0}{\bar{a}_{11}} - \frac{a_{21}^1 - a_{21}^0}{\bar{a}_{21}} \right) \\
&= \frac{p(1-p)a_{10}^0 a_{10}^1}{\bar{a}_{10}(\bar{a}_{11})^2 \bar{a}_{21}} \left( \frac{a_{11}^1}{\bar{a}_{10}} - \frac{a_{11}^0}{\bar{a}_{10}} \right) (a_{11}^0 a_{21}^1 - a_{11}^1 a_{21}^0).
\end{aligned}$$

Since we have assumed that $\frac{a_{11}^1}{a_{10}^1} \ge \frac{a_{11}^0}{a_{10}^0}$, the sign of $\Lambda'_y(0)$ is indeed the same as that of $a_{11}^1 a_{21}^0 - a_{11}^0 a_{21}^1$ (recall that $\bar{a}_{21} < 0$).

**Example 6.** We now consider the case described in Example 1, where $Qf = Q^{-1}f = \pi f - f$. In that case, we get

$$Q^{-1} \left( \frac{a_{11}}{\bar{a}_{11}} - \frac{a_{21}}{\bar{a}_{21}} \right) (s) = \frac{a_{21}^s}{\bar{a}_{21}} - \frac{a_{11}^s}{\bar{a}_{11}}$$
which entails

$$c_1 = \frac{\overline{a_{10}a_{21}}}{\overline{a_{10}\bar{a}_{21}}} + \frac{\overline{a_{11}^2}}{(\overline{a_{11}})^2} - \frac{\overline{a_{10}a_{11}}}{\overline{a_{10}\bar{a}_{11}}} - \frac{\overline{a_{11}a_{21}}}{\overline{a_{11}\bar{a}_{21}}}.$$

# 4 Annex

## 4.1 Proof of Proposition 1

*Proof.* This is essentially the proof of [22], in dual form since we work with $P_t f$ rather than the law of the process. It is organized in three steps. First, assuming formally that (8) gives an expansion in $\varepsilon$ of $P_t^\varepsilon$, we deduce the expression of $P_t^{(k)}$ and $S_t^{(k)}$. With these definitions, in a second step, we establish the bounds (6) and (7) on $P_t^{(k)}$ and $S_t^{(k)}$. Finally, from this, we obtain the bound (9) on $R_{n,t}^\varepsilon$.

**Step 1.** Fix $f \in C^\infty(M \times S)$. For all $t \ge 0$,

$$\partial_t P_t^\varepsilon f = L^\varepsilon P_t^\varepsilon f. \qquad (32)$$

Having in mind the ansatz that (8) gives an expansion in $\varepsilon$ of $P_t^\varepsilon f$, developing in $\varepsilon$ each side of the above equality, and equating the terms of the same order (treating separately the terms $P_t^{(k)} f$ and $S_t^{(k)} f$), we end up with the following equations: for all $k \ge 0$,

$$QP_t^{(0)} f = 0, \qquad (33)$$

$$\partial_t P_t^{(k)} f = L_c(P_t^{(k)} f) + QP_t^{(k+1)} f \qquad (34)$$

and

$$\partial_t S_t^{(0)} f = QS_t^{(0)} f, \qquad (35)$$

$$\partial_t S_t^{(k)} f = QS_t^{(k)} f + L_c(S_t^{(k-1)} f). \qquad (36)$$

To shorten the notation, for $k \in \mathbb{N}$, we write $u_k(t, x, s) = P_t^{(k)} f(x, s)$. Moreover, for a function $g(t, x, s)$, we use the notation $\pi g(t, x) = \int_S g(t, x, s)\pi(ds)$. Equation (33) is equivalent to say that there exists $\theta_0(t, x)$ such that for all $s \in S$, $u_0(t, x, s) = \theta_0(t, x)$. Now, if we want Equation (34) for $k = 0$ to have a solution $u_1$, this will induce a constraint, fixing the value of $\theta_0$. Indeed, this equation can be rewritten as

$$Qu_1 = \partial_t u_0 - L_c u_0. \qquad (37)$$

The left hand side averages to 0 with respect to $\pi$ (for all fixed $x \in M$ $t \ge 0$), and thus we have to impose that

$$\pi(\partial_t u_0 - L_c u_0) = 0.$$

This is equivalent to

$$\partial_t \theta_0(t, x) - \bar{F}(x) \cdot \nabla \theta_0(t, x) = 0, \qquad (38)$$
which is a transport equation with solution given by

$$ \theta_0(t, x) = \theta_0(0, \varphi_t(x)). $$

In order to identify a suitable initial condition $\theta_0(0, x) = \pi u_0(0, x)$, we formally let $\varepsilon$ and then $t$ go to 0 in

$$ \pi P_t^\varepsilon f(x) = \pi f(x) + \int_0^t \pi L_c P_r^\varepsilon f(x) dr, $$

which is obtained by integrating (32) and using that $\pi Q = 0$. This yields $\theta_0(0, x) = \pi f(x)$ and as announced,

$$ P_t^{(0)} f(x, s) = u_0(t, x, s) = \pi f(\varphi_t(x)). \quad (39) $$

Now, reasoning by induction, we assume that for some $k \ge 1$, we have constructed $u_{k-1}$ satisfying $\pi(\partial_t u_{k-1} - L_c u_{k-1}) = 0$, from which we want to define $u_k$. Recalling the definition (37) of the pseudo-inverse $Q^{-1}$, Equation (34) is then equivalent to the existence of a function $\theta_k : \mathbb{R}_+ \times \mathbb{R}^d \to \mathbb{R}$ such that

$$ u_k(t, x, s) = Q^{-1}(\partial_t u_{k-1} - L_c u_{k-1})(t, x, s) + \theta_k(t, x), \quad (40) $$

Once again, to find a suitable $\theta_k$, we use Equation (34) at order $k + 1$ that we rewrite as

$$ Qu_{k+1} = \partial_t u_k - L_c u_k. \quad (41) $$

As before, for fixed $(t, x)$, since the left-hand side averages to 0 with respect to $\pi$, we impose

$$ \pi (\partial_t u_k - L_c u_k)(t, x) = 0. \quad (42) $$

Notice that $\pi Q^{-1} = 0$, so that, averaging Equation (40) with respect to $\pi$ and differentiating with respect to time,

$$ \pi \partial_t u_k = \partial_t \theta_k. \quad (43) $$

Injecting Equations (43) and (40) in (42) yields

$$ \partial_t \theta_k(t, x) - \bar{F}(x) \cdot \nabla \theta_k(t, x) = \pi L_c \mathbf{b}_k(t, x), \quad (44) $$

where we set

$$ \mathbf{b}_k = Q^{-1}(\partial_t u_{k-1} - L_c u_{k-1}). $$

This is again a transport equation, similar to (38), but with a source term. The solution is given by

$$ \theta_k(t, x) = \theta_k(0, \varphi_t(x)) + \int_0^t \pi L_c \mathbf{b}_k(r, \varphi_{t-r}(x)) dr. \quad (45) $$

For now, by construction, for any choice of the function $\theta_k(0, \cdot)$, the function $u_k$ defined by (40) where $\theta_k$ is given by (45) satisfies (34) and (42). It remains to
specify a suitable initial condition $\theta_k(0, x)$. To do so, we need first to look at the so-called boundary layer correction terms $S^{(j)}$.

Set $v_k(t, x, s) = S_t^{(k)} f(x, s)$. Then, Equations (35) and (36) read

$$ \partial_t v_0 = Q v_0 \qquad (46) $$

and

$$ \partial_t v_k = Q v_k + L_c v_{k-1}. \qquad (47) $$

Equation (46) is simply solved as

$$ v_0(t, x, s) = Q_t v_0(0, x, s). $$

Moreover, taking $t = 0$ and letting $\varepsilon \to 0$ in (8) yields $v_0(0, x, s) + u_0(0, x, s) = f(x, s)$, hence we set

$$ v_0(t, x, s) = Q_t f(x, s) - \pi f(x), \qquad (48) $$

which indeed solves (46). Next, assuming by induction that $v_{k-1}$ has been defined for some $k \ge 1$, we can solve equation (47) to find

$$ v_k(t, x, s) = Q_t v_k(0, x, s) + \int_0^t Q_{t-r} L_c v_{k-1}(r, x, s) dr. $$

It remains to choose a suitable initial condition $v_k(0, x)$, which will be determined by the requirement of the long-time decay (7). Indeed, since $\lim_{t \to \infty} Q_t = \pi$, (7) can only hold if

$$ 0 = \pi \left( v_k(0, x) + \int_0^\infty Q_r L_c v_{k-1}(r, x) dr \right) $$

and therefore,

$$ \pi v_k(0, x) = - \int_0^\infty \pi Q_r L_c v_{k-1}(r, x) dr = - \int_0^\infty \pi L_c v_{k-1}(r, x) dr. \qquad (49) $$

Besides, from (8) applied with $t = 0$, we get that $v_k(0, x, s) + u_k(0, x, s) = 0$ for all $k \ge 1$ and, integrating (40) at $t = 0$ with respect to $\pi$, $\theta_k(0, x) = \pi u_k(0, x)$. Hence, the only choice of $\theta_k(0, \cdot)$ that is compatible with (7) is

$$ \theta_k(0, x) = \int_0^\infty \pi L_c v_{k-1}(r, x) dr. $$

Finally, $v_k(0, x, s) = -u_k(0, x, s)$ is determined using again (40) at $t = 0$. At this point, we have completely determined a definition of $P_t^{(k)}$ and $S_t^{(k)}$, which is the one that we use in the rest on the proof. It is such that (33), (34), (35) and (36) hold.

**Step 2** We prove that bounds (6) and (7) on $P_t^{(k)}$ and $S_t^{(k)}$ hold true. We work by induction on $k$, starting with $k = 0$. By compactness, recalling
the definition (39) of $u_0$, it is straightforward that for all $T > 0$, there exists $C_0(T, j, p)$ such that for all $j \ge 0$ and all $f \in C^\infty(M \times S)$,

$$ \sup_{t \in [0, T]} \| \partial_t^p u_0 \|_j \le C_0(T, j, p) \| f \|_{j+p}. $$

Let $j \ge 0$, $\alpha \in [1, N]^j$ and $p \in \mathbb{N}$. Then, by Equation (48),

$$ \partial_x^\alpha v_0 = Q_t(\partial_x^\alpha f - \pi(\partial_x^\alpha f)). $$

Using the ergodicity condition (3), we get for all $f \in C^\infty$, $t \ge 0$ and $j \in \mathbb{N}$,

$$ \|v_0(t, \cdot)\|_j \le C e^{-\gamma t} \|f\|_j. $$

For $k \ge 1$, reasoning by induction, we assume that there exist constants $C_{k-1}(T, j, p)$ and $C_{k-1}(j)$ such that

$$ \sup_{t \in [0, T]} \| \partial_t^p u_{k-1}(t, \cdot) \|_j \le C_{k-1}(T, j, p) \| f \|_{j+p+2(k-1)} $$

and for all $t \ge 0$,

$$ \| v_{k-1}(t, \cdot) \|_j \le C_{k-1}(j) e^{-\gamma t} (1 + t^{k-1}) \| f \|_{j+2(k-1)}. $$

Since $s, x \mapsto F_s(x)$ and its derivatives in $x$ are uniformly bounded over $M \times S$, it is straightforward to check that that for any $j, p \in \mathbb{N}$, and $T > 0$, there exist $C(j)$ and $C(T, j, p)$ such that for any $g \in C^\infty(M)$,

$$ \|L_c g\|_j \le C(j) \|g\|_{j+1} $$

$$ \sup_{t \in [0, T]} \| \partial_t^p (g \circ \varphi_t) \|_j \le C(T, j, p) \| g \|_{j+p}, $$

which we use extensively in the rest of this step. In the following, $T, p, j, k$ are fixed, and we denote by $C$ a constant which may change from line to line and depends only on $T, p, j, k$. First, for all $t \le T$

$$ \begin{aligned} \| \partial_t^p (\theta_k(0, \cdot) \circ \varphi_t) \|_j &\le C \| \theta_k(0, \cdot) \|_{j+p} \\ &\le C \int_0^\infty \| L_c v_{k-1}(r, \cdot) \|_{j+p} dr \\ &\le C \int_0^\infty \| v_{k-1}(r, \cdot) \|_{j+p+1} dr \\ &\le C \| f \|_{j+p+1+2(k-1)}. \end{aligned} $$

by induction. Next, using that $\partial_t^p \mathbf{b}_k = Q^{-1}(\partial_t^{p+1} u_{k-1} - L_c(\partial_t^p u_{k-1}))$, we get that for all $t \le T$

$$ \begin{aligned} \| \partial_t^p \mathbf{b}_k \|_j &\le C \left( \| \partial_t^{p+1} u_{k-1} \|_j + \| L_c (\partial_t^p u_{k-1}) \|_j \right) \\ &\le C \left( \| \partial_t^{p+1} u_{k-1} \|_j + \| \partial_t^p u_{k-1} \|_{j+1} \right) \\ &\le C \| f \|_{j+p+1+2(k-1)}, \end{aligned} $$
by induction. Now, from (45), for all $t \le T$,

$$ \begin{aligned} \|\theta_k(t, \cdot)\|_j &\le \|\theta_k(0, \cdot) \circ \varphi_t\|_j + \int_0^t \|L_c \mathbf{b}_k(s, \cdot) \circ \varphi_{t-s}\|_j ds \\ &\le C \|\theta_k(0, \cdot)\|_j + C \int_0^t \|\mathbf{b}_k(s, \cdot)\|_{j+1} ds \\ &\le C \|f\|_{j+2+2(k-1)}. \end{aligned} $$

To bound the derivatives in time of $\theta_k$, instead of using (45), it is simpler to work with (44), which gives for $p \ge 1$

$$ \begin{aligned} \|\partial_t^p \theta_k(t, \cdot)\|_j &\le \|\pi L_c \partial_t^{p-1} \theta_k(t, \cdot)\|_j + \|\pi L_c \partial_t^{p-1} \mathbf{b}_k(t, \cdot)\|_j \\ &\le C \|\partial_t^{p-1} \theta_k(t, \cdot)\|_{j+1} + C \|\partial_t^{p-1} \mathbf{b}_k(t, \cdot)\|_{j+1} \\ &\le C \|\partial_t^{p-1} \theta_k(t, \cdot)\|_{j+1} + C \|f\|_{j+p+1+2(k-1)}. \end{aligned} $$

Reasoning by induction on $p$ yields $\|\partial_t^p \theta_k(t, \cdot)\|_j \le C \|f\|_{j+p+2+2(k-1)}$. Gathering the previous bounds in (40) yields

$$ \|\partial_t^p u_k(t, \cdot)\|_j \le \|\partial_t^p \theta_k(t, \cdot)\|_{j,0} + \|\partial_t^p \mathbf{b}_k(t, \cdot)\|_j \le C \|f\|_{j+p+2k} $$

This concludes the study by induction of $u_k$. We now turn to the study of $v_k$. Using Equation (49) and omitting the dependency in $x, s$ to alleviate notations, we decompose $v_k$ as

$$ \begin{aligned} v_k(t) &= Q_t v_k(0) + \int_0^t Q_{t-r} L_c v_{k-1}(r) dr \\ &= Q_t(v_k(0) - \pi v_k(0)) + \int_0^t Q_{t-r} (L_c v_{k-1}(r) - \pi L_c v_{k-1}(r)) dr \\ &\quad - \int_t^\infty \pi L_c v_{k-1}(r) dr. \end{aligned} $$

We have

$$ \|Q_t(v_k(0) - \pi v_k(0))\|_j \le C e^{-\gamma t} \|v_k(0)\|_j \le C e^{-\gamma t} \|f\|_{j+2k}, $$

where we used that $v_k(0) = -u_k(0)$ and the previous bound on $u_k(0)$. Similarly,

$$ \begin{aligned} \left\| \int_0^t Q_{t-r} (L_c v_{k-1}(r) - \pi L_c v_{k-1}(r)) dr \right\|_j &\le C \int_0^t e^{-\gamma(t-r)} \|L_c v_{k-1}(r)\|_j dr \\ &\le C \int_0^t e^{-\gamma(t-r)} \|v_{k-1}(r)\|_{j+1} dr \\ &\le C e^{-\gamma t} (1 + t^k) \|f\|_{j+1+2(k-1)}, \end{aligned} $$
where we have used the induction hypothesis on $v_{k-1}$. Finally,

$$ 
\begin{aligned}
\int_t^\infty \|L_c v_{k-1}(r)\|_j dr &\le C \int_t^\infty \|v_{k-1}(r)\|_{j+1} dr \\
&\le C \int_t^\infty (1+r^{k-1})e^{-\gamma r} dr \|f\|_{j+1+2(k-1)} \\
&\le C(1+t^{k-1})e^{-\gamma t} \|f\|_{j+1+2(k-1)}.
\end{aligned}
$$

This proves the desired bounds on $v_k$. In particular, at this stage, we have proven (6) and (7).

**Step 3.** We now bound the remainder $R_{n,t}^\varepsilon$ given by (8) for a given $n \ge 1$. Introducing $\eta_t = (\partial_t - L_\varepsilon)R_{n,t}^\varepsilon f$ for a given $f \in C^\infty$, we see that, thanks to (33), (34), (35) and (36), all the low order terms in $\varepsilon$ vanishes (as designed) so that

$$ 
\begin{aligned}
\eta_t &= (\partial_t - L_\varepsilon) \left( P_t^\varepsilon f - \sum_{k=0}^n \varepsilon^k P_t^{(k)} f + \sum_{k=0}^n \varepsilon^k S_{\frac{t}{\varepsilon}}^{(k)} f \right) \\
&= \varepsilon^n \left( L_c P_t^{(n)} f - \partial_t P_t^{(n)} f + L_c S_{t/\varepsilon}^{(n)} f \right).
\end{aligned}
$$

The bounds obtained in the previous step yield

$$ \sup_{t \in [0,T]} \|\eta_t\|_j \le C \varepsilon^n \|f\|_{j+1+2n} $$

for some $C$ depending only on $T, j, n$. Interpreting $R_{n,t}^\varepsilon f$ as the solution of the equation $\partial_t R_{n,t}^\varepsilon f = L_\varepsilon R_{n,t}^\varepsilon f + \eta_t$ with $R_{n,0}^\varepsilon f = 0$ yields the Feynman-Kac representation

$$ R_{n,t} f(x,s) = \mathbb{E}_{x,s} \left( \int_0^t \eta_r(X(r), \sigma(r)) dr \right), $$

see [11, Theorem 6.3]. The bounds on $\eta_t$ immediately yields

$$ \sup_{t \in [0,T]} \|R_{n,t}^\varepsilon f\|_\infty \le C \varepsilon^n \|f\|_{1+2n} $$

for some $C$ depending only on $T, n$.

Now, to get similar bounds for the derivatives in space of $R_{n,t}^\varepsilon f$, we work by induction on $j$, assuming that we have already obtained that

$$ \sup_{t \in [0,T]} \|R_{n,t}^\varepsilon f\|_{j'} \le C \varepsilon^n \|f\|_{j'+1+2n} \quad (50) $$

for all $j' < j$ for some $j \ge 1$ and some $C > 0$. Indeed, then, considering a multi-index $\alpha \in [1, d]^j$, differentiating the equation $\partial_t R_{n,t}^\varepsilon = L_\varepsilon R_{n,t}^\varepsilon f + \eta_t$ yields

$$ \partial_t \partial_x^\alpha R_{n,t}^\varepsilon f = L_\varepsilon \partial_x^\alpha R_{n,t}^\varepsilon f + \partial_x^\alpha \eta_t + \sum_{|\alpha'| \le j} B_{\alpha'} \partial_x^{\alpha'} R_{n,t}^\varepsilon f, \quad (51) $$
where, for all multi-index $\alpha'$, $B_{\alpha'}(x, s)$ is a matrix whose coefficients are some derivatives in $x$ of $F_s$. In other words, considering a vector $G_t$ whose coefficients are all the partial derivatives in $x$ of $R_{n,t}^\epsilon f$ of order $j$, then $G_t$ solves an equation of the form

$$ \partial_t G_t = L_\epsilon G_t + B' G_t + \eta'_t, \qquad (52) $$

where $L_\epsilon$ acts coefficient-wise on $G_t$, the coefficients of the matrix $B'(x, s)$ are derivatives in $x$ of $F_s$ and $\eta'_t$ gathers $\partial_x^\alpha \eta_t$ and all the terms in (51) involving derivatives of $R_{n,t}^\epsilon f$ of order $|\alpha'| < j$. In particular, thanks to the bounds (6) and (7) established in step 2 and to the induction assumption, we have

$$ \sup_{t \in [0, T]} \|\eta'_t\|_{\infty} \le C \epsilon^n \|f\|_{j+1+2n} $$

for some $C > 0$ depending on $T, n, j$. We use again the Feynman-Kac formula for $G_t$ based on Equation (52), which reads

$$ G_t(x, s) = \mathbb{E}_{x,s} \left( \int_0^t E_r \eta'_{t-r} (X(r), \sigma(r)) dr \right) $$

where $(E_t)_{t \ge 0}$ is the solution of the matrix-valued ODE $\partial_t E_t = E_t B'(X(t), \sigma(t))$ (for completeness, we provide a proof of this formula in appendix). Since $B'$ is uniformly bounded, so it $E_r$ over $r \in [0, T]$, and thus

$$ \sup_{t \in [0, T]} \|G_t\|_{\infty} \le C \epsilon^n \|f\|_{j+1+2n} $$

for some $C > 0$ depending on $T, n, j$. Since this holds for all $\alpha \in [1, d]^j$, this concludes the proof by induction of (50) for all $j'$. To conclude, using that

$$ R_{n,t}^\epsilon f = R_{n+1,t}^\epsilon f + \epsilon^{n+1} P_t^{(n+1)} f + \epsilon^{n+1} S_{\frac{t}{\epsilon}}^{(n+1)} f, $$

we get that, for all $j \ge 0$,

$$ \begin{aligned} \|R_{n,t}^\epsilon f\|_j &\le \|R_{n+1,t}^\epsilon f\|_j + \epsilon^{n+1} \|P_t^{(n+1)} f\|_j + \epsilon^{n+1} \|S_{\frac{t}{\epsilon}}^{(n+1)} f\|_j \\ &\le C \epsilon^{n+1} \|f\|_{j+3+2n}, \end{aligned} $$

where we used (6), (7) and (50) (for $n+1$). This concludes the proof of (9), and thus the proof of Proposition 1. $\square$

## 4.2 Feynman Kac formula

**Lemma 2.** Let $L$ be a Markov generator on a set $S'$. For $n \ge 1$, consider $f_t(z) = (f_{1,t}(z), \dots, f_{n,t}(z)) \in \mathbb{R}^n$ for $z \in S'$, $t \ge 0$ solution of

$$ \partial_t f_t(z) = L f_t(z) + A_t(z) f_t(z) + c_t(z), \quad z \in S', \quad t \ge 0, \qquad (53) $$
where $(t, z) \mapsto A$ is a bounded $n \times n$ matrix field, $(t, z) \mapsto c$ is a bounded vector field and $Lf_t$ is to be understood component-wise. Then, for all $t \ge 0$ and $z \in S'$,

$$f_t(z) = \mathbb{E}(E_t f_0(Z_t)) + \int_0^t \mathbb{E}(E_r c_{t-r}(Z_r)) dr,$$

where $(Z_t)_{t \ge 0}$ is a Markov process associated to $L$ with $Z_0 = z$ and $(E_r)_{r \in [0, t]}$ is the solution of the $n \times n$ matrix-valued random ODE $\partial_r E_r = E_r A_{t-r}(Z_r)$ with $E_0 = I_n$, provided this is well-defined almost surely.

*Proof.* Fix $t > 0$, $z \in S'$. Then $(Z_r, E_r)_{r \ge 0}$ given in the statement of the lemma is a time-inhomogeneous Markov process with generator

$$\mathcal{L}_r g(z, e) = eA_{t-r}(z) \cdot \nabla_e g(z, e) + Lg(z, e),$$

and thus, considering the vector-valued function $h(r, z, e) = e f_{t-r}(z)$, we end up with

$$\begin{align*}
\partial_r \mathbb{E}(E_r f_{t-r}(Z_r)) &= \partial_r \mathbb{E}(h(r, Z_r, E_r)) \\
&= \mathbb{E}((\partial_r + \mathcal{L}_r) h(r, Z_r, E_r)) \\
&= \mathbb{E}(E_s (-\partial_r + A_{t-r} + L) f_{t-r}(Z_r)) \\
&= -\mathbb{E}(E_r c_{t-r}(Z_r)).
\end{align*}$$

Using that $E_0 = I$, conclusion follows from

$$f_t(z) = \mathbb{E}(E_0 f_t(Z_0)) = \mathbb{E}(E_t f_0(Z_t)) - \int_0^t \partial_r \mathbb{E}(E_r f_{t-r}(Z_r)) dr.$$

□

## 4.3 Proof of Lemma 1

We can write $\sigma^\epsilon(t) = \sigma(t/\epsilon)$, and then drop the superscript $\epsilon$ to alleviate notations; in other words, without loss of generality it is sufficient to treat the case $\epsilon = 1$ (up to replacing $\gamma$ by $\gamma/\epsilon$ in Assumption 1).

For $s \in S$, let $F_s(x) = x(a_{10}^s - a_{11}^s)x$. By definition of $p_0$ and $p_1$, we have that for all $s \in S$, $F_s(x) > 0$ for $x \in (0, p_0)$ and $F_s(x) < 0$ for $x > p_1$, which implies that for all $X_0 \ne 0$, there exists a time $T_0$ such that, for all $t \ge T_0$, $X_t \in [p_0, p_1]$. Now, on the state space $E = [p_0, p_1] \times S$, we introduce the metric

$$d((x, s), (y, s')) = 1_{s \ne s'} + 1_{s=s'} \frac{|x - y|}{p_1 - p_0},$$

and we consider on the space of probability measures on $E$ the associated Wasserstein distance defined for all $\mu, \nu$ by

$$W_d(\mu, \nu) = \inf_{\Gamma} \int d((x, s), (y, s')) d\Gamma((x, s), (x, s')),$$
where the infimum runs over all the coupling measures of $\mu$ and $\nu$. We will show that there exists $C, \delta > 0$ such that, for all $(x, s), (y, s') \in E$ and all $t \ge 0$,

$$W_d(\delta_{(x,s)}P_t, \delta_{(y,s')}P_t) \le Ce^{-\delta t}d((x, s), (y, s')), \quad (54)$$

where $P_t$ stands for $P_t^1$, the semigroup of $(X, \sigma)$ at time $t$. The above inequality implies Lemma 1 by classical arguments, since the space of probability measures on $E$ endowed with the Wasserstein metric $W_d$ is complete.

We first prove (54) for point $(x, s)$ and $(y, s')$ such that $s' = s$. In that case, a coupling of $\delta_{(x,s)}P_t$ and $\delta_{(y,s)}P_t$ is given by the random variables $((X_t^{x,s}, \sigma_t^s), (X_t^{y,s}, \sigma_t^s))$, where for all $T \ge 0$,

$$X_T^{x,s} = x + \int_0^T F_{\sigma^s(u)}(X_u^{x,s})du, \quad X_T^{y,s} = y + \int_0^T F_{\sigma^s(u)}(X_u^{y,s})du.$$

In other words, $X_t^{x,s}$ and $X_t^{y,s}$ denote the solution to (28) with respective initial conditions $x$ and $y$ and constructed with the same realisation of the processus $(\sigma(u)_{u \ge 0})$ starting from $s$.

**Lemma 3.** Almost surely, for all $x, y \in [p_0, p_1]$, $s \in S$ and $t \ge 0$,

$$|X_t^{x,s} - X_t^{y,s}| \le \frac{p_1}{p_0}e^{-\eta t}|x - y|, \quad (55)$$

where $\eta = p_0 \inf_{s \in S} a_{11}^s > 0$. This immediately yields

$$W_d(\delta_{(x,s)}P_t^{\varepsilon}, \delta_{(y,s)}P_t^{\varepsilon}) \le \frac{p_1}{p_0}e^{-\eta t}d((x, s), (y, s)).$$

*Proof.* Let $s \in S$ and write $\alpha_t = a_{10}^{\sigma^s(t)}$ and $\beta_t = a_{11}^{\sigma^s(t)}$. Let $\varphi_t(z)$ the solution of the differential equation

$$\dot{z}_t = \alpha_t - \beta_t e^{z_t}, \quad z_0 = z \quad (56)$$

Then, for all $z, z' \ge \ln(p_0)$, and all $t \ge 0$, we have

$$|\varphi_t(z) - \varphi_t(z')| \le e^{-p_0 \int_0^t \beta_u du} |z - z'|. \quad (57)$$

Indeed, the associated variational equation is

$$\frac{d}{dt}(\partial_z \varphi_t(z)) = -\beta_u(\partial_z \varphi_t(z))e^{\varphi_t(z)},$$

thus

$$\begin{aligned} \partial_z \varphi_t(z) &= \exp \left( - \int_0^t \beta_u e^{\varphi_u(z)} du \right) \\ &\le \exp \left( -p_0 \int_0^t \beta_u du \right), \end{aligned}$$
where the inequality comes from the fact that $z \ge \ln(p_0)$ implies $\varphi_t(z) \ge \ln(p_0)$ for all $t \ge 0$.

Now, consider the differential equation

$$ \dot{x}_t = x_t(\alpha_t - \beta_t x_t), \quad x_0 = x \qquad (58) $$

and set $z_t = \ln(x_t)$. Then, $z_t$ solves (56) with initial condition $\ln(x)$. Denoting by $\psi_t(x)$ the solution to (58), it holds true, using (57), that for all $x, y \in [p_0, p_1]$ and all $t \ge 0$,

$$ \begin{aligned} |\psi_t(x) - \psi_t(y)| &= |e^{\varphi_t(\ln(x))} - e^{\varphi_t(\ln(y))}| \\ &\le p_1 |\varphi_t(x) - \varphi_t(y)| \\ &\le p_1 e^{-p_0 \int_0^t \beta_u du} |\ln(x) - \ln(y)| \\ &\le \frac{p_1}{p_0} e^{-p_0 \int_0^t \beta_u du} |x - y|. \end{aligned} $$

This entails the result since $X_t^{x,s}$ and $X_t^{y,s}$ solve (58) with respective initial conditions $x$ and $y$. $\square$

We now tackle the case where the initial conditions of $\sigma$ in the two chains are distinct $s \ne s' \in S$. Let $t > 0$. Thanks to Assumption 1, there exists a coupling $(\sigma_{t/2}^s, \sigma_{t/2}^{s'})$ of $\delta_s Q_{t/2}$ and $\delta_{s'} Q_{t/2}$ with

$$ \mathbb{P} \left( \sigma_{t/2}^{s} \neq \sigma_{t/2}^{s'} \right) \le C e^{-\gamma t}. $$

Fix $x, y \in [p_0, p_1]$, and let $X_{t/2}^{x,s}$ be such that $(\sigma_{t/2}^s, X_{t/2}^{x,s})$ is distributed according to $\delta_{x,s} P_{t/2}^\epsilon$ (which is possible since regular version of the conditional laws exists in Polish spaces, see e.g. [18, Theorem 5.3]), and similarly for $X_{t/2}^{y,s'}$. Then, we define $(\sigma_r^s, X_r^{x,s}, \sigma_r^{s'}, X_r^{y,s'})_{r \ge t/2}$ as follows: under the event $\{\sigma_{t/2}^s = \sigma_{t/2}^{s'}\}$, we let $(\sigma_r^s)_{r \ge t/2}$ be a Markov process associated to $(Q_r)_{r \ge 0}$ with initial condition $\sigma_{t/2}^s$, and we set $\sigma_r^{s'} = \sigma_r^s$ for all $r \ge t/2$; otherwise, we define these processes as two independent processes associated to $(Q_r)_{r \ge 0}$ with respective initial conditions $\sigma_{t/2}^s$ and $\sigma_{t/2}^{s'}$. In both cases, we define $(X_r^{x,s}, X_r^{y,s'})_{r \ge t/2}$ as solutions to

$$ X_r^{x,s} = X_{t/2}^{x,s} + \int_{t/2}^{r} F_{\sigma^s(u)}(X_u^{x,s}) du, \quad X_r^{y,s'} = X_{t/2}^{y,s'} + \int_{t/2}^{r} F_{\sigma^{s'}(u)}(X_u^{y,s'}) du. $$

Then $(\sigma_t^s, X_t^{x,s}, \sigma_t^{s'}, X_t^{y,s'})$ is a coupling of $\delta_{(x,s)} P_t^\epsilon$ and $\delta_{(y,s)} P_t^\epsilon$ with, using (55)
and that $|X_{t/2}^{x,s} - X_{t/2}^{x,s}| \le p_1 - p_0$ almost surely,

$$ \begin{aligned} & \mathbb{E} \left( d \left( (X_t^{x,s}, \sigma_t^s), (X_t^{y,s'}, \sigma_t^{s'}) \right) \right) \\ &= \mathbb{P} \left( \sigma_{t/2}^s \neq \sigma_{t/2}^{s'} \right) + \mathbb{E} \left( \frac{|X_t^{x,s} - X_t^{y,s'}|}{p_1 - p_0} 1_{\sigma_{t/2}^s = \sigma_{t/2}^{s'}} \right) \\ &\le Ce^{-\gamma t/2} + \frac{p_1}{p_0} e^{-\eta t/2} \\ &\le \left( C + \frac{p_1}{p_0} \right) e^{-(\gamma \wedge \eta)t/2} d \left( (x, s), (y, s') \right). \end{aligned} $$

As a conclusion, we have obtained that (54) holds for all $(x, s), (y, s') \in [p_0, p_1] \times S$ for some constants $C, \delta > 0$, which concludes the proof.

## Acknowledgements

The research of P. Monmarché is supported by the French ANR grant SWIDIMS (ANR-20-CE40-0022).

## References

[1] A. Agazzi, J. C. Mattingly, and O. Melikechi. Random Splitting of Fluid Models: Unique Ergodicity and Convergence. *Communications in Mathematical Physics*, March 2023.

[2] Ludwig Arnold, Hans Crauel, and Jean-Pierre Eckmann. *Lyapunov exponents: proceedings of a conference held in Oberwolfach, May 28-June 2, 1990*. Springer, 2006.

[3] A. Ben-Israel and T. N. E. Greville. *Generalized Inverses: Theory and Applications*. CMS Books in Mathematics. Springer, 2 edition, 2003.

[4] M. Benaïm. Stochastic persistence. *arXiv preprint arXiv:1806.08450*, 2018.

[5] M. Benaïm, S. Le Borgne, F. Malrieu, and P.-A. Zitt. On the stability of planar randomly switched systems. *Ann. Appl. Probab.*, 24(1):292–311, 2014.

[6] M. Benaïm and C. Lobry. Lotka Volterra in fluctuating environment or “how switching between beneficial environments can make survival harder”. *Ann. Appl. Probab.*, 26(6):3754–3785, 2016.

[7] M. Benaïm, C. Lobry, T. Sari, and E. Strickler. A note on the top Lyapunov exponent of linear cooperative systems. *arXiv e-prints*, page arXiv:2302.05874, February 2023.

[8] M. Benaïm and E. Strickler. Random switching between vector fields having a common zero. *Ann. Appl. Probab.*, 29(1):326–375, 02 2019.
[9] Y. Chitour, G. Mazanti, P. Monmarché, and M. Sigalotti. On the gap between deterministic and probabilistic Lyapunov exponents for continuous-time linear systems. *arXiv e-prints*, page arXiv:2112.07005, December 2021.

[10] I. Chueshov. *Monotone random systems theory and applications*, volume 1779 of *Lecture Notes in Mathematics*. Springer-Verlag, Berlin, 2002.

[11] M. H. A. Davis. Piecewise-deterministic markov processes: A general class of non-diffusion stochastic models. *Journal of the Royal Statistical Society. Series B (Methodological)*, 46(3):353–388, 1984.

[12] R. J. Elliott and T. K. Siu. On markov-modulated exponential-affine bond price formulae. *Applied Mathematical Finance*, 16(1):1–15, 2009.

[13] Alessandra Faggionato, Davide Gabrielli, and M Ribezzi Crivellari. Non-equilibrium thermodynamics of piecewise deterministic markov processes. *Journal of Statistical Physics*, 137:259–304, 2009.

[14] L. Fainshil, M. Margaliot, and P. Chigansky. On the stability of positive linear switched systems under arbitrary switching laws. *IEEE Trans. Automatic. Control*, 54(4):897–899, 2009.

[15] BD Goddard, M Ottobre, KJ Painter, and I Souttar. On the study of slow-fast dynamics, when the fast process has multiple invariant measures. *arXiv preprint arXiv:2305.04632*, 2023.

[16] L. Gurvits, R. Shorten, and O. Mason. On the stability of switched positive linear systems. *IEEE Trans. Automat. Control*, 52(6):1099–1103, 2007.

[17] Haralampos Hatzikirou, Nikos I Kavallaris, and Marta Leocata. A novel averaging principle provides insights in the impact of intratumoral heterogeneity on tumor progression. *Mathematics*, 9(20):2530, 2021.

[18] O. Kallenberg. *Foundations of modern probability*. Probability and its Applications (New York). Springer-Verlag, New York, second edition, 2002.

[19] JEFFREY P. KHAROUFEH and STEVEN M. COX. Stochastic models for degradation-based reliability. *IIE Transactions*, 37(6):533–542, 2005.

[20] S. D. Lawley, J. C. Mattingly, and M.I C. Reed. Sensitivity to switching rates in stochastically switched ODEs. *Commun. Math. Sci.*, 12(7):1343–1352, 2014.

[21] F. Malrieu and P.-A. Zitt. On the persistence regime for Lotka-Volterra in randomly fluctuating environments. *ALEA Lat. Am. J. Probab. Math. Stat.*, 14(2):733–749, 2017.

[22] K. Pakdaman, M. Thieullen, and G. Wainrib. Asymptotic expansion and central limit theorem for multiscale piecewise-deterministic markov processes. *Stochastic Processes and their Applications*, 122(6):2292–2318, 2012.
[23] D. Talay and L. Tubaro. Expansion of the global error for numerical schemes solving stochastic differential equations. *Stochastic Analysis and Applications*, 8(4):483–509, 1990.