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

Pierre Monmarché111Sorbonne Université, 4 place Jussieu 75005 Paris (France),
Édouard Strickler222Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy
(France)

###### 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 $\varepsilon^{-1}$, an averaging phenomenon occurs and the
solution of the ODE converges to a deterministic ODE as $\varepsilon$
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 $\varepsilon$ 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\geqslant 0}$ solution to

$x^{\prime}(t)=F_{\sigma(t)}\left(x(t)\right)$ (1)

where $(\sigma(t))_{t\geqslant 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/\varepsilon)$ for some small $\varepsilon$.
Provided $\sigma$ is ergodic with respect to some probability measure $\pi$,
the modulated process $x(t)$ is known to converge, as $\varepsilon$ vanishes,
to the solution of the averaged ODE

$x^{\prime}(t)=\bar{F}(x(t))\,,\quad\text{where}\quad\bar{F}(x)=\int_{\mathcal{S}}F_{s}(x)\pi(\text{d}s)\,.$
(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\geqslant 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\geq 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 exemples of
applications. Finally, an appendix gathers the proofs of some intermediary
results.

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

$\mathcal{Q}_{t}f(\sigma)=\mathbb{E}_{\sigma}\left(f(\sigma_{t})\right)\,,$

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

$\left\vvvert f\right\vvvert_{j}=\sum_{|\alpha|\leqslant
j}\|\partial^{\alpha}_{x}f\|_{\infty}\,.$

######  Assumption 1.

We assume that:

  1. 1.

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

  2. 2.

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

$\|\mathcal{Q}_{t}f-\pi f\|_{\infty}\leqslant Ce^{-\gamma t}\|f\|_{\infty}$
(3)

for all bounded measurable $f$ on $\mathcal{S}$ and all $t\geqslant 0$.

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

$Q^{-1}f:=\int_{0}^{\infty}\mathcal{Q}_{t}\left(\pi f-f\right)\text{d}t$ (4)

is well-defined for all $f\in\mathcal{A}$, with $Q^{-1}f\in\mathcal{A}$ and
$\left\vvvert Q^{-1}f\right\vvvert_{j}\leqslant C/\gamma\left\vvvert
f\right\vvvert_{j}$ for all $j\in\mathbb{N}$. Moreover, using that
$\partial_{t}\mathcal{Q}_{t}f=Q\mathcal{Q}_{t}f=\mathcal{Q}_{t}Qf$ for all
$t\geqslant 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\geqslant 0}$ is a standard Poisson
process with intensity $1$, $Y_{0}=\sigma(0)$ and $(Y_{k})_{k\geqslant 1}$ is
an i.i.d. sequence of random variables distributed according to $\pi$. It is
readily checked that, in that case, $\mathcal{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 $\mathcal{S}$ is a finite
state space and $(\sigma_{t})_{t\geq 0}$ is an irreducible continuous-time
Markov chain on $\mathcal{S}$, with transition rate matrix $Q$. In that case,
$Q^{-1}$ is the _group inverse_ of $Q$, defined as the unique matrix $X$
solution333such 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. to

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

In the particular case when the cardinal of $\mathcal{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$.

##  2 High-frequency expansion

The process $(X(t),\sigma(t/\varepsilon))_{t\geqslant 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 a 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\geqslant 0$, denote by $P_{t}^{\varepsilon}$
the associated Markov semigroup, namely

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

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\geq 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\vvvert P_{t}^{(k)}f\right\vvvert_{j}\leq C\left\vvvert
f\right\vvvert_{j+2k}\,.$ (6)

  * •

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

$\left\vvvert S_{t}^{(k)}f\right\vvvert_{j}\leq Ce^{-\gamma t}\left\vvvert
f\right\vvvert_{j+2k}\,.$ (7)

  * •

For all $n\in\mathbb{N}$, $T\geq 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_{\frac{t}{\varepsilon}}^{(k)}f,$
(8)

satisfies

$\left\vvvert R_{n,t}^{\varepsilon}f\right\vvvert_{j}\leq
C\varepsilon^{n+1}\left\vvvert f\right\vvvert_{j+3+2n}\,.$ (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)=\mathcal{Q}_{t}f(x,s)-\pi f(x)$ (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,$
(11)

where

$\mathbf{b}_{1}(t,x,s)=Q^{-1}(\partial_{t}P_{t}^{(0)}f-L_{c}P_{t}^{(0)}f)(x,s).$
(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\geqslant 0$ and $f\in\mathcal{C}^{\infty}(\mathcal{M})$,

$\left\vvvert
f\circ\varphi_{t}-f\left(\bar{x}\right)\right\vvvert_{j}\leqslant
Ce^{-at}\left\vvvert f\right\vvvert_{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\geqslant 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\geq 1}$ and $(M_{k})_{k\geq 1}$ such that for all
$n\in\mathbb{N}$ and any invariant measure $\mu_{\varepsilon}$ of
$(P_{t}^{\varepsilon})_{t\geqslant 0}$,

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

Moreover,

$c_{1}=\pi L_{c}Q^{-1}\left(L_{c}h-f\right)(\bar{x})$ (14)

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

$h(x,s)=\int_{0}^{\infty}\pi\left(f(\bar{x})-f\left(\varphi_{r}(x)\right)\right)dr\,.$
(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}\left(\pi f(\bar{x})-\pi
f(\varphi_{r}(x))\right)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\geqslant 0}$. By definition of
the flow $\varphi_{t}$,

$(P_{t}^{(0)})^{n}f(x,i)=\pi f\left(\varphi_{t}^{\circ n}(x)\right)=\pi
f\left(\varphi_{tn}(x)\right).$

Hence,

$\left\vvvert(P_{t}^{(0)})^{n}(\mu_{0}f-f)\right\vvvert_{j}=\left\vvvert\pi
f(\bar{x})-\pi(f\circ\varphi_{tn})\right\vvvert_{j}\leq Ce^{-atn}\left\vvvert
f\right\vvvert_{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)$ (16)

is well-defined, in $\mathcal{A}$ and such that, for all $j\geq 0$,
$\left\vvvert\Psi_{t}\right\vvvert_{j}\leq C_{j}\left\vvvert
f\right\vvvert_{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|\leqslant\|(P_{t}^{(0)}-P_{t}^{\varepsilon})\Psi_{t}\|_{\infty}\
\leqslant C\varepsilon\left\vvvert\Psi_{t}\right\vvvert_{3}\,.$

To get the higher order expansion, write

$\displaystyle\mu_{\varepsilon}f-\mu_{0}f$
$\displaystyle=\mu_{\varepsilon}(P_{t}^{(0)}-P_{t}^{\varepsilon})\Psi_{t}$
$\displaystyle=-\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}.$
(17)

From Proposition 1 and the bounds on $\left\vvvert\Psi_{t}\right\vvvert_{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|\leqslant
C_{n}\varepsilon^{n+1}$

for some $C_{n}$. Using that, from the convergence at order $0$,

$\left|\mu_{\varepsilon}P_{t}^{(k)}\Psi_{t}-\mu_{0}P_{t}^{(k)}\Psi_{t}\right|\leqslant
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

$\left|\mu_{\varepsilon}f-\mu_{0}f+\varepsilon\mu_{0}P_{t}^{(1)}\Psi_{t}\right|\leqslant
C_{1,0}\varepsilon^{2}$ (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

$\left|\mu_{\varepsilon}P_{t}^{(k)}\Psi_{t}-c_{0,k}-c_{1,k}\varepsilon\right|\leqslant
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

$\left|\mu_{\varepsilon}f-\mu_{0}f-\varepsilon
c_{0,1}-\varepsilon^{2}c_{1,1}-\varepsilon^{2}c_{0,2}\right|\leqslant
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\,,\qquad\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=\mathcal{Q}_{r}(f-\pi f)$, from

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

integrating with respect to $\pi$, we get

$\nabla_{x}S_{r}^{(0)}\Psi_{t}(\bar{x},s)=\nabla_{x}\mathcal{Q}_{r}\left(\Psi_{t}-\pi\Psi_{t}\right)(\bar{x},s)=\mathcal{\nabla}_{x}\mathcal{Q}_{r}\left(\pi
f-f\right)(\bar{x},s).$

As a consequence,

$\displaystyle\mu_{0}(\alpha)=\alpha(\bar{x})$
$\displaystyle=\int_{0}^{\infty}\pi L_{c}S_{r}^{(0)}\Psi_{t}(\bar{x})dr$
$\displaystyle=\int_{\mathcal{S}}F_{s}(\bar{x})\cdot\left[\int_{0}^{\infty}\nabla_{x}Q_{r}\left(\pi
f-f\right)dr\right](\bar{x},s)\pi(ds)$
$\displaystyle=\int_{\mathcal{S}}F_{s}(\bar{x})\cdot\nabla_{x}Q^{-1}f(\bar{x},s)\pi(ds)$
$\displaystyle=\pi L_{c}Q^{-1}f(\bar{x}).$

Besides,

$\displaystyle\mu_{0}(\beta)$ $\displaystyle=\int_{0}^{t}\pi
L_{c}\mathbf{b}_{1}^{\Psi}(r,\bar{x})dr$
$\displaystyle=\pi\left[L_{c}\int_{0}^{t}\mathbf{b}_{1}^{\Psi}(r,\cdot)dr\right](\bar{x})$
$\displaystyle=\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})$
$\displaystyle=-\pi\left[L_{c}Q^{-1}\int_{0}^{t}L_{c}P_{r}^{(0)}\Psi_{t}dr\right](\bar{x})\,,$

where we used that $\partial_{r}P_{r}^{(0)}\Psi_{t}$ does not depend on $s$
and is thus in the kernel of $Q^{-1}$. Hence,

$\displaystyle\int_{0}^{t}\pi L_{c}\mathbf{b}_{1}^{\Psi}(s,\bar{x})ds$
$\displaystyle=-\pi\left[L_{c}Q^{-1}\frac{1}{t}\int_{0}^{t}L_{c}P_{r}^{(0)}(t\Psi_{t})dr\right](\bar{x})$
$\displaystyle\underset{t\rightarrow
0}{\longrightarrow}-\pi\left[L_{c}Q^{-1}L_{c}(\pi h)\right](\bar{x}),$

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

$\displaystyle t\Psi_{t}(x)$ $\displaystyle=$ $\displaystyle t\pi
f(\bar{x})-tf(x,s)+t\sum_{n=1}^{\infty}\pi\left(f(\bar{x})-f(\varphi_{tn}(x))\right)$
$\displaystyle\underset{t\rightarrow 0}{\longrightarrow}$
$\displaystyle\int_{0}^{\infty}\pi\left(f(\bar{x})-f\left(\varphi_{r}(x)\right)\right)dr\
=\ h(x)\,,$

and the convergence $t\Psi\rightarrow h$ holds in all norms
$\left\vvvert\cdot\right\vvvert_{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^{\prime}(t)=A\left(\sigma(t/\varepsilon)\right)z(t)$ (19)

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

######  Assumption 3.

  1. 1.

The Markov process $(\sigma(t))_{t\geqslant 0}$ is Feller.

  2. 2.

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

  3. 3.

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

The fact that the matrices are cooperative implies that $\mathbb{R}_{+}^{d}$
is fixed by (19) for $t\geqslant 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)=\mathbbm{1}\cdot z(t)>0$ with
$\mathbbm{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^{\prime}(t)=(\mathbbm{1}\cdot
A\left(\sigma(t/\varepsilon)\right)\theta(t))\rho(t)\,,\qquad\theta^{\prime}(t)=F_{\sigma(t/\varepsilon)}(\theta(t))$

with

$F_{s}(\theta)=A(s)\theta-(\mathbbm{1}\cdot A(s)\theta)\theta\,.$

It is proven in [7] that, under Assumption 3, the Markov process
$(\theta(t),\sigma(t/\varepsilon))_{t\geqslant 0}$ admits a unique invariant
measure $\mu_{\varepsilon}$ on $\Delta\times\mathcal{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_{\varepsilon},$ (20)

where

$\Lambda_{\varepsilon}:=\int_{\Delta\times\mathcal{S}}\mathbbm{1}\cdot
A(s)\theta\,\mu_{\varepsilon}(\text{d}\theta,\text{d}s)$ (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_{\varepsilon\to 0}\Lambda_{\varepsilon}=\lambda_{\max}(\bar{A})\,.$

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

######  Proposition 2.

Under Assumptions 1 and 3, there exists a sequence $(c_{k})_{k\geqslant 1}$ of
real numbers such that, for all $n\geqslant 1$,

$\Lambda_{\varepsilon}=\lambda_{\max}(\bar{A})+\sum_{k=1}^{n}c_{k}\varepsilon^{k}+\underset{\varepsilon\rightarrow
0}{o}(\varepsilon^{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_{\mathcal{S}}\bar{y}^{\intercal}Q^{-1}(A)(s)\left(\bar{x}\bar{y}^{\intercal}-I\right)A(s)\bar{x}\pi(\text{d}s)\,.$

###### Proof.

First, considering $h$ given by (15), we claim that for all
$(x,s)\in\Delta\times\mathcal{S},$

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

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

$\displaystyle\pi f\left(\varphi_{r}(x)\right)$
$\displaystyle=\frac{\mathbbm{1}\cdot\bar{A}e^{r\bar{A}}x}{\mathbbm{1}_{d}\cdot
e^{r\bar{A}}x}$ $\displaystyle=\frac{d}{dr}\ln(\mathbbm{1}\cdot
e^{r\bar{A}}x)\,.$

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

$\displaystyle\int_{0}^{\infty}\left(\pi f\left(\bar{x}\right)-\pi
f\left(\varphi_{r}(x)\right)\right)\text{d}r$
$\displaystyle=\int_{0}^{\infty}\frac{d}{dr}\left(\ln\left(\frac{\mathbbm{1}\cdot
e^{r\bar{A}}\bar{x}}{\mathbbm{1}\cdot e^{r\bar{A}}x}\right)\right)\text{d}r$
$\displaystyle=\lim_{r\to\infty}\ln\left(\frac{\mathbbm{1}\cdot
e^{r\bar{A}}\bar{x}}{\mathbbm{1}\cdot e^{r\bar{A}}x}\right)$
$\displaystyle=-\ln(\bar{y}\cdot x),$

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

$\displaystyle\nabla_{x}Q^{-1}(L_{c}h-f)(\bar{x},s)$
$\displaystyle=\frac{Q^{-1}(A)(s)\bar{x}\cdot\bar{y}}{(\bar{y}\cdot\bar{x})^{2}}\bar{y}-\frac{[Q^{-1}(A)(s)]^{\intercal}\bar{y}}{\bar{y}\cdot\bar{x}}$
$\displaystyle=(Q^{-1}(A)(s)\bar{x}\cdot\bar{y})\bar{y}-[Q^{-1}(A)(s)]^{\intercal}\bar{y}\,,$

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

$\displaystyle F_{s}(\bar{x})\cdot\nabla_{x}Q^{-1}(L_{c}h-f)(\bar{x},s)$
$\displaystyle=$ $\displaystyle\left(A(s)\bar{x}-(\mathbbm{1}\cdot
A(s)\bar{x})\bar{x}\right)\cdot\left((Q^{-1}(A)(s)\bar{x}\cdot\bar{y})\bar{y}-[Q^{-1}(A)(s)]^{\intercal}\bar{y}\right)$
$\displaystyle=$
$\displaystyle(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)]^{\intercal}\bar{y}$
$\displaystyle-(\mathbbm{1}\cdot
A(s)\bar{x})(Q^{-1}(A)(s)\bar{x}\cdot\bar{y})\bar{x}\cdot\bar{y}+(\mathbbm{1}\cdot
A(s)\bar{x})\bar{x}\cdot[Q^{-1}(A)(s)]^{\intercal}\bar{y}\,.$

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

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

We end up with

$\displaystyle F_{s}(\bar{x})\cdot\nabla_{x}Q^{-1}(L_{c}h-f)(\bar{x},s)$
$\displaystyle=$
$\displaystyle(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}$
$\displaystyle=$
$\displaystyle\bar{y}^{\intercal}Q^{-1}(A)(s)\left[\bar{x}\bar{y}^{\intercal}-I\right]A(s)\bar{x}\,.$

Integrating with respect to $\pi$ leads to

$c_{1}=\int_{\mathcal{S}}\bar{y}^{\intercal}Q^{-1}(A)(s)\left(\bar{x}\bar{y}^{\intercal}-I\right)A(s)\bar{x}\pi(\text{d}s)\,.$

∎

######  Example 3.

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

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

we get the following simple expression for $c_{1}$:

$c_{1}=p(1-p)\left[\bar{y}^{\intercal}(A_{0}-A_{1})^{2}\bar{x}-\left(\bar{y}^{\intercal}(A_{0}-A_{1})\bar{x}\right)^{2}\right].$
(24)

Indeed, recalling that $Q^{-1}=Q$ (see Example 2) and using that
$\pi_{s}Q_{s,s^{\prime}}=p(1-p)$ if $s\neq s^{\prime}$ and $-p(1-p)$ if
$s=s^{\prime}$, we get

$c_{1}=\bar{y}^{\intercal}\left(A_{0}PA_{0}+A_{1}PA_{1}-A_{0}PA_{1}-A_{1}PA_{0}\right)\bar{x},$

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

$\displaystyle\bar{y}^{\intercal}A_{0}\bar{x}\bar{y}^{\intercal}A_{1}\bar{x}+$
$\displaystyle\bar{y}^{\intercal}A_{1}\bar{x}\bar{y}^{\intercal}A_{0}\bar{x}$
$\displaystyle-\bar{y}^{\intercal}A_{0}\bar{x}\bar{y}^{\intercal}A_{0}\bar{x}-\bar{y}^{\intercal}A_{1}\bar{x}\bar{y}^{\intercal}A_{1}\bar{x}$
$\displaystyle=-\left(\bar{y}^{\intercal}(A_{0}-A_{1})\bar{x}\right)^{2}$

and

$\bar{y}^{\intercal}A_{0}A_{1}\bar{x}+\bar{y}^{\intercal}A_{1}A_{0}\bar{x}-\bar{y}^{\intercal}A_{0}A_{0}\bar{x}-\bar{y}^{\intercal}A_{1}A_{1}\bar{x}=-\bar{y}^{\intercal}(A_{0}-A_{1})^{2}\bar{x},$

which induces formula (24).

When considering the switching between matrices $(A_{i})_{i\in\mathcal{S}}$,
the following natural question arises: if all the matrices are stable (i.e.
$\lambda_{\max}(A_{i})<0$ for all $i\in\mathcal{S}$), is the switched system
(19) stable, in the sense that $\Lambda_{\varepsilon}<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 a 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},\;A_{1}=\begin{pmatrix}-10&0&10\\\ 0&-10&0\\\
0&10&-1\end{pmatrix}.$ (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:

  1. 1.

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

  2. 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 1(a) 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 1(b), where is plotted
$p\mapsto c_{1}(p)$, that on the interval $[0.3,0.5]$, we have $c_{1}(p)\geq
15$, so that in particular $c_{1}(p^{*})\geq 15>0$.

(a) $p\mapsto\lambda_{\max}(p)$

(b) $p\mapsto c_{1}(p)$

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)=\mathbbm{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)\geq\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^{*})\geq\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. 1.

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

  2. 2.

For all $p\in[0,1]$,
$\lambda_{\max}(A_{p}^{\delta})\leq\lambda_{\max}(A_{p^{*}}^{\delta})=\lambda_{\max}(p^{*})+\delta=-\frac{c_{1}^{2}(p^{*})}{8M(p^{*})}<0.$

∎

######  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}^{\intercal}\bar{A}=\lambda(\bar{A})\bar{y}^{\intercal}$ and
$\bar{y}^{\intercal}\bar{x}=1$, we get that
$\bar{y}^{\intercal}\bar{A}(\bar{x}\bar{y}^{\intercal}-I)=0$, which yields

$c_{1}=\int_{\mathcal{S}}\bar{y}^{\intercal}A^{2}(s)\bar{x}-\left(\bar{y}^{T}A(s)\bar{x}\right)^{2}\pi(ds)\,.$

Notice that, from Jensen’s inequality,

$\int_{\mathcal{S}}\left(\bar{y}^{T}A(s)\bar{x}\right)^{2}\pi(ds)\geq\left(\int_{\mathcal{S}}\bar{y}^{T}A(s)\bar{x}\pi(ds)\right)^{2}=\lambda_{\max}(\bar{A})^{2}\,,$

so that

$c_{1}\leq\int_{\mathcal{S}}\bar{y}^{\intercal}A^{2}(s)\bar{x}\pi(ds)-\lambda_{\max}(\bar{A})^{2}=\bar{y}^{\intercal}\overline{A^{2}}\bar{x}-\lambda_{\max}(\overline{A}^{2})=\bar{y}^{\intercal}\left(\overline{A^{2}}-\overline{A}^{2}\right)\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))\,,$ (26)

the vector field $\bar{F}$ is splited 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^{\varepsilon}_{n+1}=\Phi^{1}_{\tau_{1}/\varepsilon}\circ\dots\circ\Phi^{N}_{\tau_{N}/\varepsilon}(x^{\varepsilon}_{n})$

where $\Phi^{i}_{t}$ 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 $\llbracket 1,N\rrbracket$ 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}\left(\varepsilon
S_{n}/N\right)|\leqslant\max_{k\in\llbracket
1,N\rrbracket}\|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\geqslant 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\rightarrow
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}$,

$\left|f\left(\varphi_{t}(x)\right)-2\mathbb{E}\left(f(x_{\varepsilon/2}(t))\right)+\mathbb{E}\left(f(x_{\varepsilon}(t))\right)\right|\leqslant
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:

$\begin{cases}\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{cases}$
(27)

We assume that $a^{s}_{11},a^{s}_{22}>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)\,.$
(28)

We prove in the Annex the following result:

######  Lemma 1.

Assume that $p_{0}:=\inf_{s\in\mathcal{S}}\frac{a_{10}^{s}}{a_{11}^{s}}>0$,
$p_{1}:=\sup_{s\in\mathcal{S}}\frac{a_{10}^{s}}{a_{11}^{s}}<+\infty$ and
$\inf_{s\in\mathcal{S}}a_{s}^{11}>0$. Then the process
$(X(t),\sigma(t/\varepsilon))_{t\geqslant 0}$ with $X$ given by (28) and
$X(0)\neq 0$ will ultimately lie in $E:=[p_{0},p_{1}]\times\mathcal{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_{\mathcal{M}\times\mathcal{S}}(a_{20}^{s}+a_{21}^{s}x)\mu_{\varepsilon}(dx,ds)\,.$
(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))\simeq\frac{1}{t}\int_{0}^{t}\left[a^{\sigma(u/\varepsilon)}_{20}+a^{\sigma(u/\varepsilon)}_{21}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_{\mathcal{S}}\alpha^{s}\pi(ds)$ the mean with respect to
$\pi$ of a function $\alpha$ defined on $\mathcal{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\mathcal{S}}a_{10}^{s}>0$ and
$\sup_{s\in\mathcal{S}}a_{11}^{s}<+\infty$. Then, for any invariant
probability measure $\mu_{\varepsilon}$ of
$(X_{t},\sigma(t/\varepsilon))_{t\geq 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_{\mathcal{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)+\underset{\varepsilon\rightarrow
0}{o}(\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}+\underset{\varepsilon\rightarrow 0}{o}(\varepsilon),$ (30)

with

$c_{1}=\pi L_{c}Q^{-1}\left(L_{c}h-f\right)(\bar{x}),$

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

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

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

$\displaystyle L_{c}Q^{-1}\left(L_{c}h-f\right)(\bar{x},s)$
$\displaystyle=\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_{\mathcal{S}}\left(\frac{a_{10}^{s}}{\bar{a}_{10}}-\frac{a^{s}_{11}}{\bar{a}_{11}}\right)Q^{-1}\left(\frac{a_{11}}{\bar{a}_{11}}-\frac{a_{21}}{\bar{a}_{21}}\right)(s)\pi(\text{d}s)\,,$
(31)

which concludes. ∎

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

  1. 1.

If $a_{11}^{s}=\bar{a}_{11}$ and $a_{10}^{s}=\bar{a}_{10}$ for all $s$, then
$\Lambda_{y}^{\prime}(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_{\frac{\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}^{\prime}(0)=0$.

  2. 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

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle-\frac{\bar{a}_{21}\bar{a}_{10}^{2}}{\bar{a}_{11}}\int_{\mathcal{S}}\frac{a_{11}^{s}}{\bar{a}_{11}}Q^{-1}\left(\frac{a_{11}}{\bar{a}_{11}}\right)(s)\pi(\text{d}s)$
$\displaystyle=$
$\displaystyle-\frac{\bar{a}_{21}\bar{a}_{10}^{2}}{\bar{a}_{11}^{3}}\int_{\mathcal{S}}(a_{11}^{s}-\bar{a}_{11})Q^{-1}\left(a_{11}-\bar{a}_{11}\right)(s)\pi(\text{d}s)\,.$

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\geqslant 0$,

$\displaystyle
t\mathbb{E}_{\pi}\left|\frac{1}{t}\int_{0}^{t}g(\sigma(s))\text{d}s\right|^{2}$
$\displaystyle=$
$\displaystyle\frac{2}{t}\mathbb{E}_{\pi}\int_{0}^{t}\int_{s}^{t}g(\sigma(u))g(\sigma(s))\text{d}s\text{d}u$
$\displaystyle=$
$\displaystyle\frac{2}{t}\int_{0}^{t}\int_{s}^{t}\pi\left(g\mathcal{Q}_{u-s}g\right)\text{d}s\text{d}u$
$\displaystyle=$
$\displaystyle\frac{2}{t}\int_{0}^{t}\int_{0}^{t}\mathbbm{1}_{v<t-s}\pi\left(g\mathcal{Q}_{v}g\right)\text{d}s\text{d}v$
$\displaystyle=$
$\displaystyle\frac{2}{t}\int_{0}^{t}(t-v)\pi\left(g\mathcal{Q}_{v}g\right)\text{d}v$
$\displaystyle\underset{t\rightarrow\infty}{\longrightarrow}$
$\displaystyle-2\pi\left(gQ^{-1}g\right)\,,$

where we used that $\int_{0}^{t}\mathcal{Q}_{v}g$ converges to $-Q^{-1}g$ and
bounded

$\left|\int_{0}^{t}v\pi\left(g\mathcal{Q}_{v}g\right)\text{d}v\right|\leqslant\int_{0}^{\infty}Cve^{-\gamma
v}\|g\|_{\infty}^{2}\text{d}v\,.$

######  Example 5.

We consider the two-state case, i.e., $\sigma$ is Markov chain on
$\mathcal{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}}(1+\frac{a_{21}^{1}}{a_{20}^{1}}x)(1-\frac{a_{11}^{0}}{a_{10}^{0}}x)-\frac{a_{20}^{0}}{a_{10}^{0}}(1+\frac{a_{21}^{0}}{a_{20}^{0}}x)(1-\frac{a_{11}^{1}}{a_{10}^{1}}x)\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}}\geqslant\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^{\prime}}=p(1-p)$ if $s\neq s^{\prime}$ and $-p(1-p)$ if
$s=s^{\prime}$, from which we compute

$\displaystyle c_{1}$
$\displaystyle=\sum_{s,s^{\prime}\in{0,1}}\pi_{s}Q_{s,s^{\prime}}\left(\frac{a_{10}^{s}}{\bar{a}_{10}}-\frac{a_{11}^{s}}{\bar{a}_{11}}\right)\left(\frac{a_{11}^{s^{\prime}}}{\bar{a}_{11}}-\frac{a_{21}^{s^{\prime}}}{\bar{a}_{21}}\right)$
$\displaystyle=p(1-p)\Big{[}\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)$
$\displaystyle\qquad-\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)\Big{]}$
$\displaystyle=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)$
$\displaystyle=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)$
$\displaystyle=\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}}{a_{10}^{1}}-\frac{a_{11}^{0}}{a_{10}^{0}}\right)\left(a_{11}^{0}a_{21}^{1}-a_{11}^{1}a_{21}^{0}\right)\,.$

Since we have assumed that
$\frac{a_{11}^{1}}{a_{10}^{1}}\geqslant\frac{a_{11}^{0}}{a_{10}^{0}}$, the
sign of $\Lambda_{y}^{\prime}(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^{s}_{21}}{\bar{a}_{21}}-\frac{a^{s}_{11}}{\bar{a}_{11}}$

which entails

$c_{1}=\frac{\overline{a_{10}a_{21}}}{\bar{a}_{10}\bar{a}_{21}}+\frac{\overline{a_{11}^{2}}}{(\bar{a}_{11})^{2}}-\frac{\overline{a_{10}a_{11}}}{\bar{a}_{10}\bar{a}_{11}}-\frac{\overline{a_{11}a_{21}}}{\bar{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\mathcal{C}^{\infty}(\mathcal{M}\times\mathcal{S})$. For all
$t\geq 0$,

$\partial_{t}P_{t}^{\varepsilon}f=L^{\varepsilon}P_{t}^{\varepsilon}f.$ (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\geq 0$,

$QP_{t}^{(0)}f=0,$ (33)

$\partial_{t}P_{t}^{(k)}f=L_{c}\left(P_{t}^{(k)}f\right)+QP_{t}^{(k+1)}f$ (34)

and

$\partial_{t}S_{t}^{(0)}f=QS_{t}^{(0)}f,$ (35)

$\partial_{t}S_{t}^{(k)}f=QS_{t}^{(k)}f+L_{c}(S_{t}^{(k-1)}f).$ (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_{\mathcal{S}}g(t,x,s)\pi(\text{d}s)$. Equation
(33) is equivalent to say that there exists $\theta_{0}(t,x)$ such that for
all $s\in\mathcal{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}.$ (37)

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

$\pi\left(\partial_{t}u_{0}-L_{c}u_{0}\right)=0.$

This is equivalent to

$\partial_{t}\theta_{0}(t,x)-\bar{F}(x)\cdot\nabla\theta_{0}(t,x)=0,$ (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)).$ (39)

Now, reasoning by induction, we assume that for some $k\geq 1$, we have
constructed $u_{k-1}$ satisfying
$\pi\left(\partial_{t}u_{k-1}-L_{c}u_{k-1}\right)=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),$
(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}.$ (41)

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

$\pi\left(\partial_{t}u_{k}-L_{c}u_{k}\right)(t,x)=0.$ (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}.$ (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),$ (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.$ (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}=Qv_{0}$ (46)

and

$\partial_{t}v_{k}=Qv_{k}+L_{c}v_{k-1}.$ (47)

Equation (46) is simply solved as

$v_{0}(t,x,s)=\mathcal{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)=\mathcal{Q}_{t}f(x,s)-\pi f(x),$ (48)

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

$v_{k}(t,x,s)=\mathcal{Q}_{t}v_{k}(0,x,s)+\int_{0}^{t}\mathcal{Q}_{t-r}L_{c}v_{k-1}(r,x,s)dr.$

It remains to chose 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}\mathcal{Q}_{t}=\pi$, (7) can only hold if

$0=\pi\left(v_{k}(0,x)+\int_{0}^{\infty}\mathcal{Q}_{r}L_{c}v_{k-1}(r,x)dr\right)$

and therefore,

$\pi
v_{k}(0,x)=-\int_{0}^{\infty}\pi\mathcal{Q}_{r}L_{c}v_{k-1}(r,x)dr=-\int_{0}^{\infty}\pi
L_{c}v_{k-1}(r,x)dr.$ (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\geq 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\geq 0$ and all
$f\in\mathcal{C}^{\infty}(\mathcal{M}\times\mathcal{S})$,

$\sup_{t\in[0,T]}\left\vvvert\partial_{t}^{p}u_{0}\right\vvvert_{j}\leq
C_{0}(T,j,p)\left\vvvert f\right\vvvert_{j+p}.$

Let $j\geq 0$, $\alpha\in\llbracket 1,N\rrbracket^{j}$ and $p\in\mathbb{N}$.
Then, by Equation (48),

$\partial^{\alpha}_{x}v_{0}=\mathcal{Q}_{t}(\partial^{\alpha}_{x}f-\pi(\partial^{\alpha}_{x}f)).$

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

$\left\vvvert v_{0}(t,\cdot)\right\vvvert_{j}\leq Ce^{-\gamma t}\left\vvvert
f\right\vvvert_{j}.$

For $k\geq 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]}\left\vvvert\partial_{t}^{p}u_{k-1}(t,\cdot)\right\vvvert_{j}\leq
C_{k-1}(T,j,p)\left\vvvert f\right\vvvert_{j+p+2(k-1)}$

and for all $t\geq 0$,

$\left\vvvert v_{k-1}(t,\cdot)\right\vvvert_{j}\leq C_{k-1}(j)e^{-\gamma
t}(1+t^{k-1})\left\vvvert f\right\vvvert_{j+2(k-1)}.$

Since $s,x\mapsto F_{s}(x)$ and its derivatives in $x$ are uniformly bounded
over $\mathcal{M}\times\mathcal{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\mathcal{C}^{\infty}(\mathcal{M})$,

$\displaystyle\left\vvvert L_{c}g\right\vvvert_{j}$ $\displaystyle\leq$
$\displaystyle C(j)\left\vvvert g\right\vvvert_{j+1}$
$\displaystyle\sup_{t\in[0,T]}\left\vvvert\partial_{t}^{p}(g\circ\varphi_{t})\right\vvvert_{j}$
$\displaystyle\leq$ $\displaystyle C(T,j,p)\left\vvvert g\right\vvvert_{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\leq T$

$\displaystyle\left\vvvert\partial_{t}^{p}(\theta_{k}(0,\cdot)\circ\varphi_{t})\right\vvvert_{j}$
$\displaystyle\leq C\left\vvvert\theta_{k}(0,\cdot)\right\vvvert_{j+p}$
$\displaystyle\leq C\int_{0}^{\infty}\left\vvvert
L_{c}v_{k-1}(r,\cdot)\right\vvvert_{j+p}dr$ $\displaystyle\leq
C\int_{0}^{\infty}\left\vvvert v_{k-1}(r,\cdot)\right\vvvert_{j+p+1}dr$
$\displaystyle\leq C\left\vvvert f\right\vvvert_{j+p+1+2(k-1)}$

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\leq T$

$\displaystyle\left\vvvert\partial_{t}^{p}\mathbf{b}_{k}\right\vvvert_{j}$
$\displaystyle\leq
C\left(\left\vvvert\partial_{t}^{p+1}u_{k-1}\right\vvvert_{j}+\left\vvvert
L_{c}(\partial_{t}^{p}u_{k-1})\right\vvvert_{j}\right)$ $\displaystyle\leq
C\left(\left\vvvert\partial_{t}^{p+1}u_{k-1}\right\vvvert_{j}+\left\vvvert\partial_{t}^{p}u_{k-1}\right\vvvert_{j+1}\right)$
$\displaystyle\leq C\left\vvvert f\right\vvvert_{j+p+1+2(k-1)}\,,$

by induction. Now, from (45), for all $t\leq T$,

$\displaystyle\left\vvvert\theta_{k}(t,\cdot)\right\vvvert_{j}$
$\displaystyle\leq\left\vvvert\theta_{k}(0,\cdot)\circ\varphi_{t}\right\vvvert_{j}+\int_{0}^{t}\left\vvvert
L_{c}\mathbf{b}_{k}(s,\cdot)\circ\varphi_{t-s}\right\vvvert_{j}ds$
$\displaystyle\leq
C\left\vvvert\theta_{k}(0,\cdot)\right\vvvert_{j}+C\int_{0}^{t}\left\vvvert\mathbf{b}_{k}(s,\cdot)\right\vvvert_{j+1}ds$
$\displaystyle\leq C\left\vvvert f\right\vvvert_{j+2+2(k-1)}.$

To bound the derivatives in time of $\theta_{k}$, instead of using (45), it is
simpler to work with (44), which gives for $p\geq 1$

$\displaystyle\left\vvvert\partial_{t}^{p}\theta_{k}(t,.)\right\vvvert_{j}$
$\displaystyle\leqslant\left\vvvert\pi
L_{c}\partial_{t}^{p-1}\theta_{k}(t,\cdot)\right\vvvert_{j}+\left\vvvert\pi
L_{c}\partial_{t}^{p-1}\mathbf{b}_{k}(t,\cdot)\right\vvvert_{j}$
$\displaystyle\leq
C\left\vvvert\partial_{t}^{p-1}\theta_{k}(t,\cdot)\right\vvvert_{j+1}+C\left\vvvert\partial_{t}^{p-1}\mathbf{b}_{k}(t,\cdot)\right\vvvert_{j+1}$
$\displaystyle\leq
C\left\vvvert\partial_{t}^{p-1}\theta_{k}(t,\cdot)\right\vvvert_{j+1}+C\left\vvvert
f\right\vvvert_{j+p+1+2(k-1)}.$

Reasoning by induction on $p$ yields
$\left\vvvert\partial_{t}^{p}\theta_{k}(t,.)\right\vvvert_{j}\leq
C\left\vvvert f\right\vvvert_{j+p+2+2(k-1)}$. Gathering the previous bounds in
(40) yields

$\left\vvvert\partial_{t}^{p}u_{k}(t,\cdot)\right\vvvert_{j}\leq\left\vvvert\partial_{t}^{p}\theta_{k}(t,.)\right\vvvert_{j,0}+\left\vvvert\partial_{t}^{p}\mathbf{b}_{k}(t,\cdot)\right\vvvert_{j}\leq
C\left\vvvert f\right\vvvert_{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

$\displaystyle v_{k}(t)$ $\displaystyle=$
$\displaystyle\mathcal{Q}_{t}v_{k}(0)+\int_{0}^{t}\mathcal{Q}_{t-r}L_{c}v_{k-1}(r)dr$
$\displaystyle=$ $\displaystyle\mathcal{Q}_{t}(v_{k}(0)-\pi
v_{k}(0))+\int_{0}^{t}\mathcal{Q}_{t-r}(L_{c}v_{k-1}(r)-\pi L_{c}v_{k-1}(r))$
$\displaystyle-\int_{t}^{\infty}\pi L_{c}v_{k-1}(r)dr.$

We have

$\left\vvvert\mathcal{Q}_{t}(v_{k}(0)-\pi v_{k}(0))\right\vvvert_{j}\leq
Ce^{-\gamma t}\left\vvvert v_{k}(0)\right\vvvert_{j}\leq Ce^{-\gamma
t}\left\vvvert f\right\vvvert_{j+2k},$

where we used that $v_{k}(0)=-u_{k}(0)$ and the previous bound on $u_{k}(0)$.
Similarly,

$\displaystyle\left\vvvert\int_{0}^{t}\mathcal{Q}_{t-r}(L_{c}v_{k-1}(r)-\pi
L_{c}v_{k-1}(r))dr\right\vvvert_{j}$ $\displaystyle\leq$ $\displaystyle
C\int_{0}^{t}e^{-\gamma(t-r)}\left\vvvert L_{c}v_{k-1}(r)\right\vvvert_{j}dr$
$\displaystyle\leq$ $\displaystyle C\int_{0}^{t}e^{-\gamma(t-r)}\left\vvvert
v_{k-1}(r)\right\vvvert_{j+1}dr$ $\displaystyle\leq$ $\displaystyle
Ce^{-\gamma t}(1+t^{k})\left\vvvert f\right\vvvert_{j+1+2(k-1)},$

where we have used the induction hypothesis on $v_{k-1}$. Finally,

$\displaystyle\int_{t}^{\infty}\left\vvvert
L_{c}v_{k-1}(r)\right\vvvert_{j}dr$ $\displaystyle\leq
C\int_{t}^{\infty}\left\vvvert v_{k-1}(r)\right\vvvert_{j+1}dr$
$\displaystyle\leq C\int_{t}^{\infty}(1+r^{k-1})e^{-\gamma r}dr\left\vvvert
f\right\vvvert_{j+1+2(k-1)}$ $\displaystyle\leq C(1+t^{k-1})e^{-\gamma
t}\left\vvvert f\right\vvvert_{j+1+2(k-1)}.$

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\geqslant 1$. Introducing
$\eta_{t}=(\partial_{t}-L_{\varepsilon})R_{n,t}^{\varepsilon}f$ for a given
$f\in\mathcal{C}^{\infty}$, we see that, thanks to (33), (34), (35) and (36),
all the low order terms in $\varepsilon$ vanishes (as designed) so that

$\displaystyle\eta_{t}$
$\displaystyle=(\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)$
$\displaystyle=\varepsilon^{n}\left(L_{c}P_{t}^{(n)}f-\partial_{t}P_{t}^{(n)}f+L_{c}S_{t/\varepsilon}^{(n)}f\right).$

The bounds obtained in the previous step yield

$\sup_{t\in[0,T]}\left\vvvert\eta_{t}\right\vvvert_{j}\leq
C\varepsilon^{n}\left\vvvert f\right\vvvert_{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}\leq
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]}\left\vvvert
R_{n,t}^{\varepsilon}f\right\vvvert_{j^{\prime}}\leq
C\varepsilon^{n}\left\vvvert f\right\vvvert_{j^{\prime}+1+2n}$ (50)

for all $j^{\prime}<j$ for some $j\geq 1$ and some $C>0$. Indeed, then,
considering a multi-index $\alpha\in\llbracket 1,d\rrbracket^{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^{\prime}|\leq
j}B_{\alpha^{\prime}}\partial_{x}^{\alpha^{\prime}}R_{n,t}^{\varepsilon}f,$
(51)

where, for all multi-index $\alpha^{\prime}$, $B_{\alpha^{\prime}}(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}^{\varepsilon}f$ of order $j$, then $G_{t}$
solves an equation of the form

$\partial_{t}G_{t}=L_{\varepsilon}G_{t}+B^{\prime}G_{t}+\eta_{t}^{\prime},$
(52)

where $L_{\varepsilon}$ acts coefficient-wise on $G_{t}$, the coefficients of
the matrix $B^{\prime}(x,s)$ are derivatives in $x$ of $F_{s}$ and
$\eta_{t}^{\prime}$ gathers $\partial_{x}^{\alpha}\eta_{t}$ and all the terms
in (51) involving derivatives of $R_{n,t}^{\varepsilon}f$ of order
$|\alpha^{\prime}|<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]}\left\vvvert\eta_{t}^{\prime}\right\vvvert_{\infty}\leq
C\varepsilon^{n}\left\vvvert f\right\vvvert_{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}^{\prime}(X(r),\sigma(r))dr\right)$

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

$\sup_{t\in[0,T]}\left\vvvert G_{t}\right\vvvert_{\infty}\leq
C\varepsilon^{n}\left\vvvert f\right\vvvert_{j+1+2n}$

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

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

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

$\displaystyle\left\vvvert R_{n,t}^{\varepsilon}f\right\vvvert_{j}$
$\displaystyle\leq\left\vvvert
R_{n+1,t}^{\varepsilon}f\right\vvvert_{j}+\varepsilon^{n+1}\left\vvvert
P_{t}^{(n+1)}f\right\vvvert_{j}+\varepsilon^{n+1}\left\vvvert
S_{\frac{t}{\varepsilon}}^{(n+1)}f\right\vvvert_{j}$ $\displaystyle\leq
C\varepsilon^{n+1}\left\vvvert f\right\vvvert_{j+3+2n},$

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

∎

###  4.2 Feynman Kac formula

######  Lemma 2.

Let $L$ be a Markov generator on a set $\mathcal{S}^{\prime}$. For $n\geqslant
1$, consider $f_{t}(z)=(f_{1,t}(z),\dots,f_{n,t}(z))\in\mathbb{R}^{n}$ for
$z\in\mathcal{S}^{\prime}$, $t\geqslant 0$ solution of

$\partial_{t}f_{t}(z)=Lf_{t}(z)+A_{t}(z)f_{t}(z)+c_{t}(z),\quad
z\in\mathcal{S}^{\prime},\ t\geqslant 0\,,$ (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\geqslant 0$ and $z\in\mathcal{S}^{\prime}$,

$f_{t}(z)=\mathbb{E}\left(E_{t}f_{0}(Z_{t})\right)+\int_{0}^{t}\mathbb{E}\left(E_{r}c_{t-r}(Z_{r})\right)\text{d}r\,,$

where $(Z_{t})_{t\geqslant 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\mathcal{S}^{\prime}$. Then $(Z_{r},E_{r})_{r\geqslant 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)=ef_{t-r}(z)$, we
end up with

$\displaystyle\partial_{r}\mathbb{E}\left(E_{r}f_{t-r}(Z_{r})\right)$
$\displaystyle=$
$\displaystyle\partial_{r}\mathbb{E}\left(h(r,Z_{r},E_{r})\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left(\left(\partial_{r}+\mathcal{L}_{r}\right)h(r,Z_{r},E_{r})\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left(E_{s}\left(-\partial_{r}+A_{t-r}+L\right)f_{t-r}(Z_{r})\right)$
$\displaystyle=$ $\displaystyle-\mathbb{E}\left(E_{r}c_{t-r}(Z_{r})\right).$

Using that $E_{0}=I$, conclusion follows from

$f_{t}(z)=\mathbb{E}\left(E_{0}f_{t}(Z_{0})\right)=\mathbb{E}\left(E_{t}f_{0}(Z_{t})\right)-\int_{0}^{t}\partial_{r}\mathbb{E}\left(E_{r}f_{t-r}(Z_{r})\right)\text{d}r\,.$

∎

###  4.3 Proof of Lemma 1

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

For $s\in\mathcal{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\mathcal{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}\neq 0$, there exists a time $T_{0}$ such that, for all $t\geq T_{0}$,
$X_{t}\in[p_{0},p_{1}]$. Now, on the state space
$E=[p_{0},p_{1}]\times\mathcal{S}$, we introduce the metric

$d((x,s),(y,s^{\prime}))=1_{s\neq
s^{\prime}}+1_{s=s^{\prime}}\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

$\mathcal{W}_{d}(\mu,\nu)=\inf_{\Gamma}\int
d((x,s),(y,s^{\prime}))d\,\Gamma((x,s),(x,s^{\prime})),$

where the infimum runs overs 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^{\prime})\in E$ and all $t\geq 0$,

$\mathcal{W}_{d}(\delta_{(x,s)}P_{t},\delta_{(y,s^{\prime})}P_{t})\leq
Ce^{-\delta t}d((x,s),(y,s^{\prime})),$ (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
$\mathcal{W}_{d}$ is complete.

We first prove (54) for point $(x,s)$ and $(y,s^{\prime})$ such that
$s^{\prime}=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^{s}_{t}),(X_{t}^{y,s},\sigma^{s}_{t}))$, where for all
$T\geq 0$,

$X_{T}^{x,s}=x+\int_{0}^{T}F_{\sigma^{s}(u)}(X_{u}^{x,s})du,\quad
X_{T}^{x,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\geq 0}$ starting from $s$.

######  Lemma 3.

Almost surely, for all $x,y\in[p_{0},p_{1}]$, $s\in\mathcal{S}$ and
$t\geqslant 0$,

$|X_{t}^{x,s}-X_{t}^{y,s}|\leq\frac{p_{1}}{p_{0}}e^{-\eta t}|x-y|,$ (55)

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

$\mathcal{W}_{d}(\delta_{(x,s)}P_{t}^{\varepsilon},\delta_{(y,s)}P_{t}^{\varepsilon})\leq\frac{p_{1}}{p_{0}}e^{-\eta
t}d((x,s),(y,s)).$

###### Proof.

Let $s\in\mathcal{S}$ and write $\alpha_{t}=a^{\sigma^{s}(t)}_{10}$ and
$\beta_{t}=a^{\sigma^{s}(t)}_{11}$. Let $\varphi_{t}(z)$ the solution of the
differential equation

$\dot{z}_{t}=\alpha_{t}-\beta_{t}e^{z_{t}},\quad z_{0}=z$ (56)

Then, for all $z,z^{\prime}\geq\ln(p_{0})$, and all $t\geq 0$, we have

$|\varphi_{t}(z)-\varphi_{t}(z^{\prime})|\leq
e^{-p_{0}\int_{0}^{t}\beta_{u}du}|z-z^{\prime}|.$ (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

$\displaystyle\partial_{z}\varphi_{t}(z)$
$\displaystyle=\exp\left(-\int_{0}^{t}\beta_{u}e^{\varphi_{u}(z)}du\right)$
$\displaystyle\leq\exp\left(-p_{0}\int_{0}^{t}\beta_{u}du\right),$

where the inequality comes from the fact that $z\geq\ln(p_{0})$ implies
$\varphi_{t}(z)\geq\ln(p_{0})$ for all $t\geq 0$.

Now, consider the differential equation

$\dot{x}_{t}=x_{t}(\alpha_{t}-\beta_{t}x_{t}),\quad x_{0}=x$ (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\geq 0$,

$\displaystyle|\psi_{t}(x)-\psi_{t}(y)|$
$\displaystyle=|e^{\varphi_{t}(\ln(x))}-e^{\varphi_{t}(\ln(y))}|$
$\displaystyle\leq p_{1}|\varphi_{t}(x)-\varphi_{t}(y)|$ $\displaystyle\leq
p_{1}e^{-p_{0}\int_{0}^{t}\beta_{u}du}|\ln(x)-\ln(y)|$
$\displaystyle\leq\frac{p_{1}}{p_{0}}e^{-p_{0}\int_{0}^{t}\beta_{u}du}|x-y|.$

This entails the result since $X_{t}^{x,s}$ and $X_{t}^{y,s}$ solve (58) with
respective initial conditions $x$ and $y$. ∎

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

$\mathbb{P}\left(\sigma_{t/2}^{s}\neq\sigma_{t/2}^{s^{\prime}}\right)\leqslant
Ce^{-\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}^{\varepsilon}$ (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^{\prime}}$. Then, we define
$(\sigma_{r}^{s},X_{r}^{x,s},\sigma_{r}^{s^{\prime}},X_{r}^{y,s^{\prime}})_{r\geqslant
t/2}$ as follows: under the event
$\\{\sigma_{t/2}^{s}=\sigma_{t/2}^{s^{\prime}}\\}$, we let
$(\sigma_{r}^{s})_{r\geqslant t/2}$ be a Markov process associated to
$(\mathcal{Q}_{r})_{r\geqslant 0}$ with initial condition $\sigma_{t/2}^{s}$,
and we set $\sigma_{r}^{s^{\prime}}=\sigma_{r}^{s}$ for all $r\geqslant t/2$;
otherwise, we define these processes as two independent processes associated
to $(\mathcal{Q}_{r})_{r\geqslant 0}$ with respective initial conditions
$\sigma_{t/2}^{s}$ and $\sigma_{t/2}^{s^{\prime}}$. In both cases, we define
$(X_{r}^{x,s},X_{r}^{y,s^{\prime}})_{r\geqslant 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^{\prime}}=X_{t/2}^{y,s^{\prime}}+\int_{t/2}^{r}F_{\sigma^{s^{\prime}}(u)}(X_{u}^{y,s^{\prime}})du.$

Then
$(\sigma_{t}^{s},X_{t}^{x,s},\sigma_{t}^{s^{\prime}},X_{t}^{y,s^{\prime}})$ is
a coupling of $\delta_{(x,s)}P_{t}^{\varepsilon}$ and
$\delta_{(y,s)}P_{t}^{\varepsilon}$ with, using (55) and that
$|X_{t/2}^{x,s}-X_{t/2}^{x,s}|\leqslant p_{1}-p_{0}$ almost surely,

$\displaystyle\mathbb{E}\left(d((X_{t}^{x,s},\sigma_{t}^{s}),(X_{t}^{y,s^{\prime}},\sigma_{t}^{s^{\prime}}))\right)$
$\displaystyle=$
$\displaystyle\mathbb{P}\left(\sigma_{t/2}^{s}\neq\sigma_{t/2}^{s^{\prime}}\right)+\mathbb{E}\left(\frac{|X_{t}^{x,s}-X_{t}^{y,s^{\prime}}|}{p_{1}-p_{0}}1_{\sigma_{t/2}^{s}=\sigma_{t/2}^{s^{\prime}}}\right)$
$\displaystyle\leq$ $\displaystyle Ce^{-\gamma
t/2}+\frac{p_{1}}{p_{0}}e^{-\eta t/2}$ $\displaystyle\leq$
$\displaystyle\left(C+\frac{p_{1}}{p_{0}}\right)e^{-(\gamma\wedge\eta)t/2}d\left((x,s),(y,s^{\prime})\right)\,.$

As a conclusion, we have obtained that (54) holds for all
$(x,s),(y,s^{\prime})\in[p_{0},p_{1}]\times\mathcal{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. Automat. 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.l 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.