Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?

Abstract

When studying the dynamics of trait distribution of populations in a heterogeneous environment, classical models from quantitative genetics choose to look at its system of moments, specifically the first two ones. Additionally, in order to close the resulting system of equations, they often assume the local trait distributions are Gaussian [see for instance Ronce and Kirkpatrick (Evolution 55(8):1520–1531, 2001. https://doi.org/10.1111/j.0014-3820.2001.tb00672.x.37)]. The aim of this paper is to introduce a mathematical framework that follows the whole trait distribution (without prior assumption) to study evolutionary dynamics of sexually reproducing populations. Specifically, it focuses on complex traits, whose inheritance can be encoded by the infinitesimal model of segregation (Fisher in Trans R Soc Edinb 52(2):399–433, 1919. https://doi.org/10.1017/S0080456800012163). We show that it allows us to derive a regime in which our model gives the same dynamics as when assuming Gaussian local trait distributions. To support that, we compare the stationary problems of the system of moments derived from our model with the one given in Ronce and Kirkpatrick (Evolution 55(8):1520–1531, 2001. https://doi.org/10.1111/j.0014-3820.2001.tb00672.x.37) and show that they are equivalent under this regime and do not need to be otherwise. Moreover, under this regime of equivalence, we show that a separation bewteen ecological and evolutionary time scales arises. A fast relaxation toward monomorphism allows us to reduce the complexity of the system of moments, using a slow-fast analysis. This reduction leads us to complete, still in this regime, the analytical description of the bistable asymmetrical equilibria numerically found in Ronce and Kirkpatrick (Evolution 55(8):1520–1531, 2001. https://doi.org/10.1111/j.0014-3820.2001.tb00672.x.37). More globally, we provide explicit modelling hypotheses that allow for such local adaptation patterns to occur.

Appendices

Appendices

A System of moments derived from our model

Here, we derive the system of moments (2) from (1). In the preliminary computations, we will omit the time and deme dependency for the sake of clarity. We will then denote $\varvec{n}$ the trait distribution density, $\varvec{N}$ the size of the population, $\varvec{\overline{z}}$ the mean trait, $\varvec{\sigma }^2$ the mean variance, $\varvec{\psi }^3$ the third central moment and $\varvec{\theta }$ the optimal phenotype.

Preliminary integration of the selection term We have:

$$\begin{aligned} \displaystyle \int _\mathbb {R}(\varvec{z}-\varvec{\theta })^2\varvec{n}(\varvec{z}) d\varvec{z}&= \displaystyle \int _{\mathbb {R}} \left[ (\varvec{z}-\varvec{\overline{z}})^2+(\varvec{\overline{z}}-\varvec{\theta })^2+2(\varvec{z}-\varvec{\overline{z}})(\varvec{\overline{z}}-\varvec{\theta })\right] \varvec{n}(\varvec{z})\,d\varvec{z}\\&= \varvec{\sigma }^2 \varvec{N}+(\varvec{\overline{z}}-\varvec{\theta })^2 \varvec{N}, \end{aligned}$$

and:

$$\begin{aligned}&\displaystyle \int _\mathbb {R}(\varvec{z}-\varvec{\overline{z}})(\varvec{z}-\varvec{\theta })^2\varvec{n} d\varvec{z} \\&\quad = \displaystyle \int _{\mathbb {R}} \left[ (\varvec{z}-\varvec{\overline{z}})^3+(\varvec{z}-\varvec{\overline{z}})(\varvec{\overline{z}}-\varvec{\theta })^2+2(\varvec{z}-\varvec{\overline{z}})^2(\varvec{\overline{z}}-\varvec{\theta })\right] \varvec{n}(\varvec{z})\,d\varvec{z}\\&\quad = 2\varvec{\sigma }^2 (\varvec{\overline{z}}-\varvec{\theta }) \varvec{N}+\varvec{\psi } \varvec{N}. \end{aligned}$$

Size of the subpopulations Recalling that $\varvec{N_i}(\varvec{t}) = \displaystyle \int _\mathbb {R}\varvec{n_i}(\varvec{t},\varvec{z})\,d\varvec{z}$, we get from the preliminary computations by integrating (1):

$$\begin{aligned} \frac{d\varvec{N_i}}{d\varvec{t}}&= \displaystyle \int _\mathbb {R}\frac{\partial \varvec{n_i}}{\partial \varvec{t}}(\varvec{t},\varvec{z}) d\varvec{z}\\&=\displaystyle \int _\mathbb {R}\varvec{r}\varvec{\mathcal {B}}_{\varvec{\sigma }} (\varvec{n_i})(\varvec{t},\varvec{z}) - \varvec{g}(\varvec{z}-\varvec{\theta _i})^2\varvec{n_i}(\varvec{t},\varvec{z}) \\&\qquad - \varvec{\kappa } \varvec{N_i}(\varvec{t})\varvec{n_i}(\varvec{t},\varvec{z})+\varvec{m}\left( \varvec{n}(\varvec{t},\varvec{z})-\varvec{n}(\varvec{t},\varvec{z})\right) \,d\varvec{z} \\&= \left[ \varvec{r}- \varvec{\kappa }\varvec{N_i}(\varvec{t})-\varvec{g}(\varvec{\overline{z}_i}(\varvec{t})-\varvec{\theta _i})^2 - \varvec{g}\varvec{\sigma _i}^2\right] \varvec{N_{i}}(\varvec{t})+\varvec{m}\big (\varvec{N_j}(\varvec{t})-\varvec{N_i}(\varvec{t})\big ). \end{aligned}$$

Local mean trait Recalling that $\varvec{\overline{z}_{i}}(\varvec{t}) =\frac{1}{\varvec{N_{i}}(\varvec{t})}\displaystyle \int _\mathbb {R}\varvec{z}\,\varvec{n_i}(\varvec{t},\varvec{z})\,d\varvec{z}$, we have, thanks to the preliminary computations:

$$\begin{aligned} \frac{d\varvec{z_i}}{d\varvec{t}}&= \frac{1}{\varvec{N_i}}\displaystyle \int _\mathbb {R}\varvec{z}\frac{\partial \varvec{n_i}}{\partial \varvec{t}}(\varvec{t},\varvec{z}) d\varvec{z} - \frac{1}{\varvec{N_i}^2}\frac{d\varvec{N_i}}{d\varvec{t}}\displaystyle \int _\mathbb {R}\varvec{z} \varvec{n_i}(\varvec{t},\varvec{z}) d\varvec{z}\\&=\frac{1}{\varvec{N_i}}\displaystyle \int _\mathbb {R}(\varvec{z}-\varvec{\overline{z}_i})\frac{\partial \varvec{n_i}}{\partial \varvec{t}}(\varvec{t},\varvec{z}) d\varvec{z}\\&=\frac{1}{\varvec{N_i}}\displaystyle \int _\mathbb {R}(\varvec{z}-\varvec{\overline{z}_i}){\left[ - \varvec{g}(\varvec{z}-\varvec{\theta _i})^2\varvec{n_i}(\varvec{t},\varvec{z}) +\varvec{m}\left( \varvec{n_j}(\varvec{t},\varvec{z})-\varvec{n_i}(\varvec{t},\varvec{z})\right) \right] d\varvec{z}} \\&= 2\varvec{g}\varvec{\sigma _i}^2(\varvec{\theta _i}-\varvec{\overline{z}_i})-\varvec{g}\varvec{\psi _i}^3 +\varvec{m}\frac{\varvec{N_j}}{\varvec{N_i}}(\varvec{\overline{z}_j}-\varvec{\overline{z}_i}). \end{aligned}$$

B Equilibria of a dynamical system under the infinitesimal model of reproduction with random mating only

In this subsection, we show that (7) admits any Gaussian of variance $\varepsilon ^2$ as equilibrium. That is equivalent to state that:

Proposition B.1

For $\mu \in \mathbb {R}$, the Gaussian distribution $G_{\mu ,\varepsilon ^2}$ of mean $\mu$ and variance $\varepsilon ^2$ is a fixed point of the operator $\mathcal {B}_\varepsilon$, namely:

$$\begin{aligned} \mathcal {B}_\varepsilon (G_{\mu ,\varepsilon ^2}) = G_{\mu ,\varepsilon ^2}. \end{aligned}$$

Proof

We can first notice that ${\mathcal {B}_\varepsilon }$ can be written using a double convolution product: $\square$

Lemma 10

For $f \in \mathcal {L}_1(\mathbb {R}), \displaystyle \int _\mathbb {R}{f} \ne 0$, we have:

$$\begin{aligned} \mathcal {B}_{\varepsilon }(f) = \frac{4}{\displaystyle \int _\mathbb {R}f(z')\,dz'}G_{0,\frac{\varepsilon ^2}{2}}*F*F, \end{aligned}$$

where $F: z \mapsto f(2z)$.

Proof

(Proof of Lemma 10) For $f \in \mathcal {L}_1(\mathbb {R}), \displaystyle \int _\mathbb {R}f \ne 0$, a straight-forward computation yields:

$$\begin{aligned} \mathcal {B}_{\varepsilon }(f)(z)&= \frac{1}{\sqrt{\pi }\varepsilon }\iint _{\mathbb {R}^2} \exp \left[ \frac{-(z-\frac{z_1+z_2}{2})^2}{\varepsilon ^2}\right] \frac{f(z_1)f(z_2)}{\displaystyle \int _\mathbb {R}f(z')\,dz'}dz_1 dz_2\\&= \frac{1}{\displaystyle \int _\mathbb {R}f(z')\,dz'}\int _\mathbb {R}\int _{\mathbb {R}} G_{0,\frac{\varepsilon ^2}{2}}\big ((z-\frac{z_1}{2})-\frac{z_2}{2}\big ) F(\frac{z_2}{2})\,dz_2{\,}F(\frac{z_1}{2})\,dz_1 \\&= \frac{2}{\displaystyle \int _\mathbb {R}f(z')\,dz'}\int _\mathbb {R}G_{0,\frac{\varepsilon ^2}{2}}*F(z-\frac{z_1}{2}) F(\frac{z_1}{2})\,dz_1\\&= \frac{4}{\displaystyle \int _\mathbb {R}f(z')\,dz'}{\,}G_{0,\frac{\varepsilon ^2}{2}}*F*F(z). \end{aligned}$$

$\square$

If $f = G_{\mu ,\varepsilon ^2}$, then we find $F = \frac{1}{2}\times G_{\frac{\mu }{2},\frac{\varepsilon ^2}{4}}$. Besides, as the convolution product of two Gaussian kernels $G_{\mu _1,\sigma ^2_1}$ and $G_{\mu _2,\sigma ^2_2}$ is the Gaussian kernel $G_{\mu _1+\mu _2,\sigma _1 ^2+\sigma _2 ^2}$, Proposition B.1 is a corollary of the previous lemma.

C Formal expansion within the exponential formalism for $n_{\varepsilon }$

In this subsection, we will remove the deme dependency for the sake of clarity. To formally derive (9), let us consider the following formal expansion of $U_\varepsilon$ with regard to successive orders of $\varepsilon ^2$:

$$\begin{aligned} U_\varepsilon = u_0+\varepsilon ^2u_\varepsilon . \end{aligned}$$

The aim is to characterize $u_0$ thanks to the behaviour of the reproduction term when $\varepsilon \ll 1$, which we expect neither to diverge nor to vanish:

$$\begin{aligned}&\frac{\mathcal {B}_{\varepsilon }(n_{\varepsilon })}{n_{\varepsilon }}(z) \\&\quad = \frac{1}{\sqrt{\pi }\varepsilon }\iint _{\mathbb {R}^2} \frac{\exp \left[ \frac{1}{\varepsilon ^2}\left[ - \left[ z -\frac{z_1+z_2}{2}\right] ^2 +u_{0}(z) - u_{0}(z_1) - u_{0}(z_2)\right] \right] \exp \left[ u_\varepsilon (z) - u_\varepsilon (z_1) - u_\varepsilon (z_2)\right] dz_1 dz_2}{\int _{\mathbb {R}} \exp \left[ -\frac{u_{0}(z')}{\varepsilon ^2}-u(z')\right] dz'} \end{aligned}$$

Then, we have several considerations to make. First, if we assume that $u_{0}$ reaches its minimum at a non degenerate point $z^*$, then the following modified expression of the denominator:

$$\begin{aligned} \frac{1}{\sqrt{\pi }\varepsilon }\int _{\mathbb {R}} \exp \left[ -\frac{1}{\varepsilon ^2}\left[ u_{0}(z') - \min u_{0}\right] -u(z')\right] dz', \end{aligned}$$

will have its integrand concentrate around the minimum of $u_{0}$ and will converge as $\varepsilon \ll 1$. Therefore it is relevant to introduce this minimum both at the numerator and the denominator.

Then, since we expect the numerator not to diverge nor to vanish uniformly as $\varepsilon \ll 1$, we need that:

$$\begin{aligned} \forall z \in \mathbb {R}, \underset{(z_1,z_2)}{\max }\left[ - \left( z -\frac{z_1+z_2}{2}\right) ^2+u_{0}(z) - u_{0}(z_1) - u_{0}(z_2) +\min u_{0} \right] = 0. \nonumber \\ \end{aligned}$$

(27)

As shown in Bouin et al. (2018), thanks to some convexity arguments, this leads necessarily to choose $u_{0}$ as a quadratic function in z, hence its decomposition:

$$\begin{aligned} u_{0}(z) = {u(z^*)} +{\frac{(z - z^*)^2}{2}}, \end{aligned}$$

(28)

where $z^*$ is realizing the minimum of $u_{0}$. Note that $u(z^*) = 0$, due to the Laplace method of integration, since:

$$\begin{aligned} N_{\varepsilon } = \frac{1}{\sqrt{2\pi }\varepsilon } \int _{\mathbb {R}} \exp \left[ -\frac{U_{\varepsilon }(z)}{\varepsilon ^2}\right] dz \underset{\varepsilon \rightarrow 0}{{\approx }} \frac{\exp \left[ -\frac{u(z^*)}{\varepsilon ^2}\right] }{\sqrt{U_\varepsilon ''(z^*)}}. \end{aligned}$$

So either $u(z^*)=0$, either there is extinction or explosion of the population size. That yields (9).

Convexity arguments from Bouin et al. (2018). Let us recall the arguments of convexity involved in Bouin et al. (2018) to show that functional constraint (27) leads in our case to $u_0$ being quadratic:

• 1. First, they show that $u_0$ has some regularities (continuous and has left and right derivative everywhere), for (27) implies that $z\mapsto u_0(z) -z^2$ is concave as minimum of affine functions:

  $$\begin{aligned} \forall z\in \mathbb {R},\quad u_0(z) -z^2 = \underset{(z_1,z_2)}{\min }\left[ -z(z_1+z_2) + \frac{\left( z_1+z_2\right) ^2}{4} + u_{0}(z_1) + u_{0}(z_2) \right] . \end{aligned}$$

• 2. Next, they introduce the Legendre convex conjugate

  $$\begin{aligned}\hat{u_0} : p \mapsto \underset{z\in \mathbb {R}}{\sup }\left[ (z-z^*)p -u(z)\right] ,\end{aligned}$$

  and show that it satisfies the following functional equality, by commuting the different $\sup$ operators while computing $\hat{u}_0(p)$ using (27):

  $$\begin{aligned} \forall p \in \mathbb {R},\quad \hat{u_0}(p) = \frac{p^2}{4}+ 2\,\hat{u_0}\left( \frac{p}{2}\right) . \end{aligned}$$

  (29)

• 3. As $\hat{u_0}$ is convex by definition, it is continuous and admits left and right derivative everywhere. Moreover, $\hat{u_0}$ has a minimum in 0 and $\hat{u_0}(0) = -u(z^*) = 0$. Therefore (29) implies by recursion:

  $$\begin{aligned} \forall p >0\quad (\text {resp.} <0), \quad \hat{u}_0(p) = \frac{p^2}{2}+ \hat{u_0}'(0^+)\,p \quad (\text {resp.} \;\hat{u_0}'(0^-)\,p). \end{aligned}$$

  (30)

  Note that 0 being a minimum of $\hat{u}_0$ implies that: $\hat{u_0}'(0^-)\le 0\le \hat{u_0}'(0^+)$.

• 4. The next step aims at showing that $u_0$ is equal to its convex bi-conjugate

  $$\begin{aligned} \hat{\hat{u_0}} : z \mapsto \underset{p\in \mathbb {R}}{\sup }\left[ p\,(z-z^*) -\hat{u_0}(p)\right] , \end{aligned}$$

  which is computable from (30):

  $$\begin{aligned} \hat{\hat{u_0}} : z\mapsto \begin{aligned} {\left\{ \begin{array}{ll} \frac{(z-z^*-\hat{u}_0(0^-))^2}{2}\text { if }z<z^*+\hat{u}_0(0^-)\\ 0\qquad \quad \text { if }z^*+\hat{u}_0(0^-)\le z \le z^*+\hat{u}_0(0^+)\\ \frac{(z-z^*-\hat{u}_0(0^+))^2}{2}\text { if }z>z^*+\hat{u}_0(0^+). \end{array}\right. } \end{aligned} \end{aligned}$$

  (31)

  Standard convexity analysis shows also that $\hat{\hat{u_0}}$ is the lower convex envelope of $u_0$. The first implication is that $u_0$ and $\hat{\hat{u_0}}$ coincide on $\mathbb {R}\backslash [z^*+\hat{u}_0(0^-),z^*+\hat{u}_0(0^+)]$, because $\hat{\hat{u_0}}$ is strictly convex there. The second implication is that $u_0\left( z^*+\hat{u}_0(0^+)\right) = \hat{\hat{u}}_0\left( z^*+\hat{u}_0(0^+)\right) =0$ (resp. $z^*+\hat{u}_0(0^-)$), since $z^*+\hat{u}_0(0^+)$ (resp. $z^*+\hat{u}_0(0^-)$) is an extremal point of the graph of $\hat{\hat{u}}_0$. One can show using (27) that the midpoint between any zeros of $u_0$ is still a zero of $u_0$ (recall that $u_0\ge 0$). Hence, by density and continuity of $u_0$, $u_0$ vanishes on $[z^*+\hat{u}_0(0^-),z^*+\hat{u}_0(0^+)]$.

• 5. Finally, since $u_0$ satisfies (31) and we need $N_\varepsilon$ not to explode when $\varepsilon$ vanishes, we necessarily obtain that $\hat{u}_0(0^-) = \hat{u}_0(0^+)$. Hence $u_0$ quadratic.

D Formal approximations of the trait distributions moments in the regime of small variance $\varepsilon ^2 \ll 1$

This appendix is dedicated to formally explain (12). We remove the time and the deme dependency for the sake of clarity. We denote $n_\varepsilon$ the trait distribution density, $N_\varepsilon$ the size of the population, $\overline{z}_\varepsilon$ the mean trait, $\sigma _\varepsilon ^2$ the variance and $\psi _\varepsilon$ the third central moment. Let us also recall that the computations are performed using the exponential formalism introduced in (10) while considering the following formal expansion of $u_\varepsilon$ in the regime of small variance:

$$\begin{aligned} u_\varepsilon = u+\varepsilon ^2\, v+\mathcal {O}(\varepsilon ^4). \end{aligned}$$

Size of population We have:

$$\begin{aligned} N_\varepsilon&= \displaystyle \int _\mathbb {R}n_\varepsilon (z)\,dz\\&= \displaystyle \int _\mathbb {R}\frac{1}{\sqrt{2\pi }\varepsilon }e^{-\frac{(z-z^*)^2}{2\varepsilon ^2}}e^{-u(z)-\varepsilon ^2\, v(z)+\mathcal {O}(\varepsilon ^4)}dz\\&= \displaystyle \int _\mathbb {R}\frac{e^{-\frac{y^2}{2}}}{\sqrt{2\pi }}e^{-u(z^*+\varepsilon y)-\varepsilon ^2 v(z^*+\varepsilon y) +\mathcal {O}(\varepsilon ^4)}dy \quad \quad \left( y:=\frac{z-z^*}{\varepsilon }\right) \\&= \displaystyle \int _\mathbb {R}\frac{e^{-\frac{y^2}{2}}}{\sqrt{2\pi }}e^{-[u(z^*)+\varepsilon y u'(z^*) + \frac{\varepsilon ^2 y^2}{2}u''(z^*)+\frac{\varepsilon ^3 y^3}{6}u'''(z^*)+\mathcal {O}(\varepsilon ^4)]-\varepsilon ^2 v(z^*)-\varepsilon ^3yv'(z^*) +\mathcal {O}(\varepsilon ^4)}dy\\&= \displaystyle \int _\mathbb {R}\frac{e^{-\frac{y^2}{2}}}{\sqrt{2\pi }}e^{-u(z^*)}e^{-\left[ \varepsilon y u'(z^*) + \varepsilon ^2 \left[ \frac{y^2u''(z^*)}{2}+v(z^*)\right] +\varepsilon ^3\left[ \frac{y^3}{6}u'''(z^*)+yv'\right] +\mathcal {O}(\varepsilon ^4)\right] }dy\\&= \displaystyle \int _\mathbb {R}\frac{e^{-\frac{y^2}{2}}}{\sqrt{2\pi }}e^{-u(z^*)} \left[ 1-\varepsilon y u'(z^*) - \varepsilon ^2 \left[ \frac{y^2u''(z^*)}{2}+v(z^*)\right] -\varepsilon ^3\left[ \frac{y^3}{6}u'''(z^*)-yv'(z^*)\right] \right. \\&\quad \quad \quad \quad \left. +\frac{1}{2}\left[ \varepsilon ^2y^2u'(z^*)^2 + \varepsilon ^3\left[ y^3u'(z^*)u''(z^*)+2yu'(z^*)v(z^*)\right] \right] - \frac{\varepsilon ^3y^3u'(z^*)^3}{6}+ \mathcal {O}(\varepsilon ^4)\right] \\&=e^{-u(z^*)}\left[ 1+\varepsilon ^2\left[ \frac{u'^2(z^*)}{2}-\frac{u''(z^*)}{2}-v(z^*)\right] \right] + \mathcal {O}(\varepsilon ^4), \end{aligned}$$

from the computations of the moments of a Gaussian.

Mean trait Similarly as above, we have:

$$\begin{aligned} \overline{z}_\varepsilon&= \displaystyle \int _\mathbb {R}z \frac{n_\varepsilon }{N_\varepsilon }dz\\&= \frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}z\frac{1}{\sqrt{2\pi }\varepsilon }e^{-\frac{(z-z^*)^2}{2\varepsilon ^2}}e^{-u(z)-\varepsilon ^2\, v(z)+\mathcal {O}(\varepsilon ^4)}dz\\&= \frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}(z^*+\varepsilon y)\frac{e^{-\frac{y^2}{2}}}{\sqrt{2\pi }}e^{-u(z^*+\varepsilon y)-\varepsilon ^2 v(z^*+\varepsilon y) +\mathcal {O}(\varepsilon ^4)}dy, \quad \quad \left( y:=\frac{z-z^*}{\varepsilon }\right) \\&= \frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}(z^*+\varepsilon y) \frac{e^{-\frac{y^2}{2}}}{\sqrt{2\pi }}e^{-u(z^*)} \left[ 1-\varepsilon y u'(z^*) + \varepsilon ^2 \left[ \frac{y^2u'(z^*)^2}{2}-\frac{y^2u''(z^*)}{2}-v(z^*)\right] \right. \\&\quad \quad \left. +\varepsilon ^3\left[ -\frac{y^3}{6}u'''(z^*)-yv'(z^*)+\frac{y^3u'(z^*)u''(z^*)}{2}+yu'(z^*)v(z^*) - \frac{3y^3u'(z^*)^3}{6}\right] +\mathcal {O}(\varepsilon ^4)\right] \\&= \frac{e^{-u(z^*)}\left[ z^*\left( 1+\varepsilon ^2\left[ \frac{u'^2(z^*)}{2}-\frac{u''(z^*)}{2}-v(z^*)\right] \right) -\varepsilon ^2u'(z^*)\right] +\mathcal {O}(\varepsilon ^4)}{e^{-u(z^*)}\left( 1+\varepsilon ^2\left[ \frac{u'^2(z^*)}{2}-\frac{u''(z^*)}{2}-v(z^*)\right] \right) +\mathcal {O}(\varepsilon ^4)}\\&=z^*-\varepsilon ^2 u'(z^*)+ \mathcal {O}(\varepsilon ^4). \end{aligned}$$

Variance Using the previous formal computations and methodology, we get:

$$\begin{aligned} \sigma ^2_\varepsilon&= \frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}(z-\overline{z}_\varepsilon )^2 n_\varepsilon (z) dz\\&=\frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}\left[ (z-z^*)^2+(z^*-\overline{z}_\varepsilon )^2+2(z-z^*)(z^*-\overline{z}_\varepsilon )\right] n_\varepsilon (z) dz\\&= \frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}\left[ \varepsilon ^2 y^2+2\varepsilon ^3yu'(z^*)+\mathcal {O}(\varepsilon ^4)\right] \left[ 1-\varepsilon yu'+\mathcal {O}(\varepsilon ^2)\right] e^{-u(z^*)} \frac{e^{-\frac{y^2}{2}}}{\sqrt{2\pi }} dy\\&= \frac{\varepsilon ^2e^{-u(z^*)}}{e^{-u(z^*)}\left[ 1+\mathcal {O}(\varepsilon ^2)\right] }\\&= \varepsilon ^2+ \mathcal {O}(\varepsilon ^4). \end{aligned}$$

Third central moment We compute, using the same change in variable $y:= \frac{z-z^*}{\varepsilon }$:

$$\begin{aligned} \psi ^3_\varepsilon&= \frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}(z-\overline{z}_\varepsilon )^3 n_\varepsilon (z) dz\\&=\frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}\left[ (z-z^*)^3+(z^*-\overline{z}_\varepsilon )^3+3(z-z^*)^2(z^*-\overline{z}_\varepsilon )+3(z-z^*)(z^*-\overline{z}_\varepsilon )^2\right] \\&\quad n_\varepsilon (z) dz\\&=\frac{1}{N_\varepsilon }\displaystyle \int _\mathbb {R}\left[ \varepsilon ^3y^3+\mathcal {O}(\varepsilon ^4)\right] \left[ e^{-u(z^*)}+\mathcal {O}(\varepsilon )\right] dz\\&= \mathcal {O}(\varepsilon ^4). \end{aligned}$$

E Fast/slow system: proof of Theorem 3.1

This appendix is dedicated to prove Theorem 3.1.

Let $(z^*_0,\bar{Y}^*_0)\in \mathbb {R}\times \Omega$ (we recall that $\Omega =(\mathbb {R}_+^*)^2\times \mathbb {R}$) be on the slow manifold, ie. such that ${G(z^*_0,\bar{Y}^*_0)} = 0$. From Lemma 6 of fast relaxation towards the slow manifold, the jacobian matrix $J_G(z^*_0,\bar{Y}^*_0)$ is invertible. Consequently, the implicit function theorem gives us U open neighbourhood of $z^*_0$ in $\mathbb {R}$, V open neighbourhood of $(z^*_0,\bar{Y}^*_0)$ in $\mathbb {R}\times \Omega$ and $\phi \in C^\infty (U,V)$ such that :

$$\begin{aligned} \forall (z^*,\bar{Y}^*)\in V, \, {G(z^*,\bar{Y}^*,0)} = 0 \implies \bar{Y}^*=\phi (z^*). \end{aligned}$$

Hence, we can define a notation that we shall use henceforth:

$$\begin{aligned} \forall z\in U, J_{z} := J_G(z,\phi (z)). \end{aligned}$$

If K is a compact subset of U such that $z^*_0 \in \mathring{K}$, we can define the Cauchy problem $(E_0)$ by the following :

$$\begin{aligned} {(E_0)\quad } \begin{aligned} {\left\{ \begin{array}{ll} \frac{dz^*}{dt}=-2gz^*(t)+F\left( \phi (z^*(t))\right) ,\\ z^*(0) = z^*_0, \end{array}\right. } \end{aligned} \end{aligned}$$

(32)

for $t\le t^*$, that we define as the following:

$$\begin{aligned} t^*:=\inf \{t>0, z^*(t) \notin K\}. \end{aligned}$$

It is similar to (20) with the initial conditions $(z^*(0),\bar{Y}^*_0) = (z^*_0,\phi (z^*_0))$. A essential part of the proof relies in the fact that we can define the following uniform positive constant, thanks to Lemma 6 of fast relaxation:

$$\begin{aligned} \lambda _K = -\frac{1}{2}\underset{z\in K}{\max } \{\lambda \in \mathbf {Sp}(J_z)\}>0. \end{aligned}$$

As the first step, we state the following lemma whose proof will be provided at the end of this appendix. It defines a uniform control constant $\gamma >0$:

Lemma 11

There exists $\gamma >0$ such that:

$$\begin{aligned} \underset{z \in K, \, s\ge 0}{\max } {\left| \left| \left| e^{\lambda _Ks}e^{J_z s} \right| \right| \right| }\le \gamma . \end{aligned}$$

$({\left| \left| \left| \cdot \right| \right| \right| }_{\mathcal {M}_3(\mathbb {R})}$ is noted ${\left| \left| \left| \cdot \right| \right| \right| })$.

The next step is to show the convergence of solutions of $(P_\varepsilon )$ (19) towards those of $(P_0)$ (20) on a time interval, yet to be defined, that will be shown to be uniform with regard to $\varepsilon$ and the initial conditions, provided that they are small enough. For that purpose, it is more convenient to consider the system $(R_\varepsilon )$ verified by the residuals $r^\varepsilon _z(t) = z_\varepsilon (t)-z^*(t)$ and $r^\varepsilon _Y(t) = \bar{Y}_\varepsilon (t)-\bar{Y}^*(t)$:

$$\begin{aligned} {(R_\varepsilon )\quad } \begin{aligned} {\left\{ \begin{array}{ll} \varepsilon ^2 \frac{dr^\varepsilon _Y}{dt}= {G(z^*(t)}+r^\varepsilon _z(t),\bar{Y}^*(t)+r^\varepsilon _Y(t))-{G(z^*(t),\bar{Y}^*(t))}-\varepsilon ^2\frac{d\bar{Y}^*}{dt}+\varepsilon ^2\nu _{N,\varepsilon }(t),\\ \\ \frac{dr^\varepsilon _z}{dt} = -2gr^\varepsilon _z(t)+ F(\bar{Y}^*(t)+r^\varepsilon _Y)-F(\bar{Y}^*(t))+\varepsilon ^2\nu _{z,\varepsilon }(t),\\ \\ (r^\varepsilon _z(0),r^\varepsilon _Y(0)) = (z^\varepsilon _0-z^*_0,\bar{Y}^\varepsilon _0-\bar{Y}^*_0), \end{array}\right. } \end{aligned} \end{aligned}$$

(33)

and introduce some further definitions.

Because K is a compact set, there exists $\delta _K>0$ such that the following set is a compact subset of V:

$$\begin{aligned} \bar{K}_{\delta _K} = \{(z,\bar{Y}) \in \mathbb {R}\times \Omega | \exists z^*\in K,\, |(z,\bar{Y}) -(z^*,\phi (z^*))|\le \delta _K\} \subset V. \end{aligned}$$

Let us consider from now $(z^\varepsilon _0,N^\varepsilon _0) \in \bar{K}_{\delta _K}$. Then we define $\Delta = \min \left( \frac{\lambda _K}{4C\gamma },\delta _K\right)$ and $T=\min (t^*,\frac{\lambda _K}{4C'\gamma })$, where:

$$\begin{aligned} C=\max \left( \Vert \partial ^2_{\bar{Y}} G\Vert _{{\infty ,\bar{K}_{\delta _K}}} ,\Vert \partial _z G\Vert _{{\infty ,\bar{K}_{\delta _K}}}, \Vert \partial _{\bar{Y}} F\Vert _ {\infty ,\Pi _\Omega (\bar{K}_{\delta _K})}\right) ) \end{aligned}$$

and :

$$\begin{aligned} C'=\underset{t\le t^*}{\max }{\left| \left| \left| \partial _t J_{z^*(t)} \right| \right| \right| }, \end{aligned}$$

where $\Pi _\Omega$ is the projection from $\mathbb {R}\times \Omega$ on $\Omega$. One can notice from these definitions and from Lemma 11, that $\gamma , \Delta , T, \lambda _K, C, C'$ do not depend on $\varepsilon$ and are uniform on $[0,t^*]$. Specifically taking $\Delta \le \frac{\lambda _K}{4C\gamma }$ and $T\le \frac{\lambda _K}{4C'\gamma }$ will turn out to be important in the proof.

On the time region [0, T], we will show that we can control explicitly the various perturbed terms that arise. We can now state the following proposition, whose proof constitutes the core of the resolution of the problem:

Proposition E.1

As $\max (\varepsilon ,|r^\varepsilon _z(0)|,|r^\varepsilon _Y(0)|) \rightarrow 0$, $(\bar{Y}_\varepsilon ,z_\varepsilon )$ converges toward $(\bar{Y}^*,z^*)$ uniformly on [0, T].

For the final step, we will show that we can reiterate the process on each interval of time $[jT,\min \{(j+1)T,t^*\}]$ with $\forall j \le \lfloor \frac{t^*}{T}\rfloor , jT\le t^*_\varepsilon$. Thus, for sufficiently small $\varepsilon$ and initial conditions, the control remains valid until $t^*$, hence Theorem 3.1.

For convenience, we will denote by $f*g\,(t)$ the convolution product of f and g at time $t>0$ :

$$\begin{aligned} f*g\, (t) = \int _0^t f(\tau )g(t-\tau ) d \tau . \end{aligned}$$

Proof

(Proof of Proposition E.1)

Let $\varepsilon \in ]0,1]$. Let us define an auxiliary time $t^*_\varepsilon$:

$$\begin{aligned} t^*_\varepsilon = \min \left( t^*,\inf \{t>0, |r^\varepsilon _z|+|r^\varepsilon _Y|> \Delta \}\right) . \end{aligned}$$

It ensures that the perturbed trajectory stays inside of $\bar{K}_{\delta _K}$ when $t\le t^*_\varepsilon$.

Let us highlight the main steps of the proof:

• 1. preliminary controls on $r^\varepsilon _Y$ by $|r^\varepsilon _Y(0)|$ and $\frac{1}{\varepsilon ^2}|r^\varepsilon _z|*e^{-\frac{\lambda _K}{2\varepsilon ^2}\cdot }$ thanks to the regularity of G, the fast relaxation properties (Lemma 6 and Lemma 11) and Gronwall’s lemma.

• 2. control $|r^\varepsilon _z|$ by $|r^\varepsilon _z(0)|$ and $|r^\varepsilon _Y|$.

• 3. finish the control on $r^\varepsilon _Y$ by using the latter and Gronwall’s lemma.

• 1. For $t\le \min (T,t^*_\varepsilon )$, we can introduce new terms in the equation from (33) on $r^\varepsilon _Y$ :

  $$\begin{aligned} \frac{dr^\varepsilon _Y}{dt}&=\frac{J_{z^*(0)}}{\varepsilon ^2}r^\varepsilon _Y+\frac{1}{\varepsilon ^2}\left[ {G(z^*(t),\bar{Y}^*(t)+r^\varepsilon _Y(t))-G(z^*(t),\bar{Y}^*(t))}-J_{z^*(0)}r^\varepsilon _Y\right] \\&\quad + \frac{1}{\varepsilon ^2}\left[ {G(z^*(t)+r^\varepsilon _z(t),\bar{Y}^*(t)+r^\varepsilon _Y(t))-G(z^*(t),\bar{Y}^*(t)+r^\varepsilon _Y(t))}\right] \\ {}&\quad -\phi '(z^*(t))(-2gz^*(t)+F(\phi (z^*(t))))+\nu _{N,\varepsilon }(t)\\ \\&= \frac{J_{z^*(0)}}{\varepsilon ^2}r^\varepsilon _Y +A_1(t)+A_2(t)+A_3(t). \end{aligned}$$

  Since $t\le \min (T,t^*_\varepsilon )$ and G is $C^\infty$ on $\bar{K}_{\delta _K}\times [0,1]$, we can control $A_1$:

  $$\begin{aligned} |A_1(t)|\le & {} \frac{1}{\varepsilon ^2}\left[ |{G(z^*(t),\bar{Y}^*(t)+r^\varepsilon _Y(t))-G(z^*(t),\bar{Y}^*(t)}-J_{z^*(t)}r^\varepsilon _Y|\right] \\&+\frac{1}{\varepsilon ^2}\left[ {\left| \left| \left| J_{z^*(t)}-J_{z^*(0)} \right| \right| \right| }\,|r^\varepsilon _Y(t)|\right] \\\le & {} \frac{1}{\varepsilon ^2}\left[ \Vert \partial ^2_{\bar{Y}}G\Vert _ {\infty ,{\bar{K}_{\delta _K}}}\,|r^\varepsilon _Y(t)|^2+T\,\underset{t\le t^*}{\max }{\left| \left| \left| \partial _t J_{z^*(t)} \right| \right| \right| }\,|r^\varepsilon _Y(t)|\right] \\\le & {} \frac{1}{\varepsilon ^2} (C\Delta +C'T) |r^\varepsilon _Y(t)|, \end{aligned}$$

  and $A_2$:

  $$\begin{aligned} |A_2(t)|= & {} \frac{1}{\varepsilon ^2}|{G(z^*(t)+r^\varepsilon _z(t),\bar{Y}^*(t)+r^\varepsilon _Y(t))}-{G(z^*(t),\bar{Y}^*(t)+r^\varepsilon _Y(t))}|\\ {}\le & {} \frac{1}{\varepsilon ^2}{\Vert \partial _z G\Vert _ {\infty ,{\bar{K}_{\delta _K}}}\,|r^\varepsilon _z(t)|\le \frac{C}{\varepsilon ^2}|r^\varepsilon _z(t)|}, \end{aligned}$$

  and $A_3$:

  $$\begin{aligned} |A_3(t)| = |-\phi '(z^*(t))(-2gz^*(t)+F(\phi (z^*(t))))+\nu _{N,\varepsilon }(t)|\le C'', \end{aligned}$$

  for some constant $C''$ independent of $\varepsilon$ and $z^*(0)\in K$. Using Duhamel formulas, we get, for $t\le \min (T,t^*_\varepsilon )$:

  $$\begin{aligned} r^\varepsilon _Y(t) = e^{\frac{J_{z^*(0)}t}{\varepsilon ^2}}r^\varepsilon _Y(0) + \left[ e^{\frac{J_{z^*(0)}\cdot }{\varepsilon ^2}} *(A_1+A_2+A_3) \right] \,(t). \end{aligned}$$

  (34)

  Hence, applying Lemma 11 yields:

  $$\begin{aligned} |r^\varepsilon _Y(t)|&\le \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _Kt}{\varepsilon ^2}}+\frac{\gamma }{\varepsilon ^2}\left[ \left( C|r^\varepsilon _z|+(C\Delta +C'T)|r^\varepsilon _Y|\right) *e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }\right] (t) \\&\quad + \varepsilon ^2\gamma \frac{{C''}}{\lambda _K}\\&\le A^{r^\varepsilon _z}(t) + \frac{\gamma (C\Delta +C'T)}{\varepsilon ^2}\int _0^t|r^\varepsilon _Y(\tau )|\,e^{\frac{\lambda _K}{\varepsilon ^2}(\tau -t)}d\tau , \end{aligned}$$

  where $A^{r^\varepsilon _z} (t) := \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _Kt}{\varepsilon ^2}}+\frac{\gamma C}{\varepsilon ^2}\left( |r^\varepsilon _z|*e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }\right) (t) + \varepsilon ^2\gamma \frac{{C''}}{\lambda _K}$. Applying Gronwall inequality to $r^\varepsilon _Y(t)e^{\frac{\lambda _Kt}{\varepsilon ^2}}$ yields:

  $$\begin{aligned} |r^\varepsilon _Y(t)|\le A^{r^\varepsilon _z}(t) + \frac{\gamma (C\Delta +C'T)}{\varepsilon ^2}\left[ A^{r^\varepsilon _z}*e^{\left( \frac{-\lambda _K}{\varepsilon ^2}+\frac{\gamma (C\Delta +C'T)}{\varepsilon ^2}\right) \cdot }\right] (t). \end{aligned}$$

  (35)

  Having fixed $\Delta \le \frac{\lambda _K}{4C\gamma }$ and $T\le \frac{\lambda _K}{4C'\gamma }$ in the preliminaries ensures that $e^{\left( \frac{-\lambda _K}{\varepsilon ^2}+\frac{\gamma (C\Delta +C'T)}{\varepsilon ^2}\right) \cdot }$ defines a negative exponential term, that we can dominate by $e^{-\frac{\lambda _K}{2\varepsilon ^2} \cdot }$. Hence:

  $$\begin{aligned} |r^\varepsilon _Y(t)|\le A^{r^\varepsilon _z}(t) + \left[ A^{r^\varepsilon _z}*\frac{\lambda _K}{2\varepsilon ^2}e^{-\frac{\lambda _K}{2\varepsilon ^2} \cdot }\right] (t). \end{aligned}$$

  (36)

  Making $A^{r^\varepsilon _z}$ explicit gives:

  $$\begin{aligned}&|r^\varepsilon _Y(t)|\le \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _Kt}{\varepsilon ^2}}+\frac{\gamma C}{\varepsilon ^2}|r^\varepsilon _z|*e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }(t) + \varepsilon ^2\gamma \frac{{C''}}{\lambda _K}\nonumber \\&\qquad + \left[ \left( \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }+\frac{\gamma C}{\varepsilon ^2}\left[ |r^\varepsilon _z|*e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }\right] + \varepsilon ^2\gamma \frac{{C''}}{\lambda _K} \right) *\left( \frac{\lambda _K}{2\varepsilon ^2}e^{-\frac{\lambda _K}{2\varepsilon ^2} \cdot }\right) \right] (t)\nonumber \\&\quad \le \gamma |r^\varepsilon _Y(0)|\left[ e^{-\frac{\lambda _Kt}{\varepsilon ^2}}+ e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }*\left( \frac{\lambda _K}{2\varepsilon ^2}e^{-\frac{\lambda _K}{2\varepsilon ^2} \cdot }\right) (t)\right] +\varepsilon ^2\gamma \frac{{C''}}{\lambda _K}(\left( 1+\int _0^t\frac{\lambda _K}{2\varepsilon ^2}e^{-\frac{\lambda _K}{2\varepsilon ^2}(\tau -t)}dt\right) \nonumber \\&\qquad +\frac{\gamma C}{\varepsilon ^2}|r^\varepsilon _z|*\left( e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }+ e^{-\frac{\lambda _K}{\varepsilon ^2\cdot }}*\frac{\lambda _K}{2\varepsilon ^2}e^{-\frac{\lambda _K}{2\varepsilon ^2} \cdot }\right) (t), \end{aligned}$$

  (37)

  thanks to the associativity of the convolution product. One can compute that, for $t\ge 0$:

  $$\begin{aligned} e^{-\frac{\lambda _Kt}{\varepsilon ^2}}+e^{-\frac{\lambda _K}{\varepsilon ^2}\cdot }*\left( \frac{\lambda _K}{2\varepsilon ^2}e^{-\frac{\lambda _K}{2\varepsilon ^2} \cdot }\right) (t)&= e^{-\frac{\lambda _Kt}{\varepsilon ^2}}+\frac{\lambda _K}{2\varepsilon ^2}\int _0^te^{-\frac{\lambda _K}{\varepsilon ^2}\tau }e^{-\frac{\lambda _K}{2\varepsilon ^2}(t-\tau )}d\tau \\&= e^{-\frac{\lambda _Kt}{\varepsilon ^2}}+\frac{\lambda _K}{2\varepsilon ^2}\int _0^te^{-\frac{\lambda _K}{2\varepsilon ^2}(t+\tau )}d\tau = e^{-\frac{\lambda _K}{2\varepsilon ^2}t}. \end{aligned}$$

  Hence, replacing those terms in (37) yields:

  $$\begin{aligned} |r^\varepsilon _Y(t)|\le \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _K}{2\varepsilon ^2}t}+2\varepsilon ^2\gamma \frac{{C''}}{\lambda _K}+\frac{C\gamma }{\varepsilon ^2}|r^\varepsilon _z|*e^{-\frac{\lambda _K}{2\varepsilon ^2}\cdot }(t). \end{aligned}$$

  (38)

• 2. The next step is to gain similarly some control on $|r^\varepsilon _z|$. Using Duhamel formula on the equation from (33) on $r^\varepsilon _z$ gives, for $t\le \min (T,t^*_\varepsilon )$:

  $$\begin{aligned} r^\varepsilon _z(t) = r^\varepsilon _z(0)e^{-2gt}+ \left( \left[ F(N^*+r^\varepsilon _Y)-F(N^*)+\varepsilon ^2\nu _{z,\varepsilon }\right] *e^{-2g\cdot }\right) (t), \end{aligned}$$

  which yields:

  $$\begin{aligned} |r^\varepsilon _z(t)| \le |r^\varepsilon _z(0)|e^{-2gt}+\varepsilon ^2\frac{\Vert \nu _{z,\varepsilon }\Vert _\infty }{2g}+\Vert \partial _{\bar{Y}} F\Vert _ {\infty ,\Pi _\Omega (\bar{K}_{\delta _K})}\left( |r^\varepsilon _Y|*e^{-2g\cdot }\right) (t). \end{aligned}$$

  Hence:

  $$\begin{aligned} |r^\varepsilon _z(t)| \le |r^\varepsilon _z(0)|e^{-2gt}+\varepsilon ^2\frac{\Vert \nu _{z,\varepsilon }\Vert _\infty }{2g}+C\left( |r^\varepsilon _Y|*e^{-2g\cdot }\right) (t). \end{aligned}$$

  (39)

  At that point, it is clear that it is sufficient to control $|r^\varepsilon _Y|$ and $|r^\varepsilon _z(0)|$ in order to control $|r^\varepsilon _z(t)|$ for sufficiently small $\varepsilon$.

• 3. Plugging the latter in (38) gives:

  $$\begin{aligned}&|r^\varepsilon _Y(t)|\le \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _K}{2\varepsilon ^2}t}+\frac{C\gamma }{\varepsilon ^2}|r^\varepsilon _z(0)|\left( e^{-2g\cdot } *e^{-\frac{\lambda _K}{2\varepsilon ^2}\cdot }\right) (t)+\varepsilon ^2\frac{C\gamma \Vert \nu _{z,\varepsilon }\Vert _\infty }{\lambda _K g}\nonumber \\&\quad +2\varepsilon ^2\gamma \frac{{C''}}{\lambda _K}+\frac{\gamma C^2}{\varepsilon ^2}\left[ |r^\varepsilon _Y|*\left( e^{-2g\cdot }*e^{-\frac{\lambda _K}{2\varepsilon ^2}\cdot } \right) \right] (t). \end{aligned}$$

  (40)

  Similarly as the computation above, we have, for $\varepsilon ^2<\min (\frac{\lambda _K}{8g},1)$ and $t\ge 0$:

  $$\begin{aligned} e^{-2g\cdot }*e^{-\frac{\lambda _K}{2\varepsilon ^2}\cdot }(t) = \frac{1}{\frac{\lambda _K}{2\varepsilon ^2}-2g}\left( e^{-2gt}-e^{-\frac{\lambda _K}{2\varepsilon ^2}t}\right) \le \frac{4\varepsilon ^2}{\lambda _K}e^{-2gt}. \end{aligned}$$

  Hence, for $\varepsilon ^2<\min (\frac{\lambda _K}{8g},1)$, we get from (40):

  $$\begin{aligned} |r^\varepsilon _Y(t)|&\le \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _K}{2\varepsilon ^2}t}+\frac{2\gamma C}{\lambda _K}|r^\varepsilon _z(0)|e^{-2gt}+\varepsilon ^2\frac{C\gamma \Vert \nu _{z,\varepsilon }\Vert _\infty }{\lambda _K g}+2\varepsilon ^2\gamma \frac{{C''}}{\lambda _K}\\ {}&\quad +\frac{2\gamma C^2}{\lambda _K}\left( |r^\varepsilon _Y|*e^{-2g\cdot } \right) (t)\\&\le C^\varepsilon _0(t) +\frac{2\gamma C^2}{\lambda _K}\left( |r^\varepsilon _Y|*e^{-2g\cdot }\right) (t), \end{aligned}$$

  where we define: $C^\varepsilon _0(t):=\gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _K}{2\varepsilon ^2}t}+\frac{2\gamma C}{\lambda _K}|r^\varepsilon _z(0)|e^{-2gt}+\varepsilon ^2\frac{C\gamma \Vert \nu _{z,\varepsilon }\Vert _\infty }{\lambda _K g}+2\varepsilon ^2\gamma \frac{{C''}}{\lambda _K}$. Using once again Gronwall inequality on $|r^\varepsilon _Y|e^{2g\cdot }$ yields:

  $$\begin{aligned} |r^\varepsilon _Y(t)|\le C^\varepsilon _0(t) +\frac{2\gamma C^2}{\lambda _K}\left( C^\varepsilon _0*e^{\left( -2g+\frac{2\gamma C^2}{\lambda _K}\right) \cdot }\right) (t). \end{aligned}$$

  (41)

  Recalling that:

  $$\begin{aligned} C^\varepsilon _0(t) = \gamma |r^\varepsilon _Y(0)|e^{-\frac{\lambda _K}{2\varepsilon ^2}t}+\frac{2\gamma C}{\lambda _K}|r^\varepsilon _z(0)|e^{-2gt}+\varepsilon ^2\frac{C\gamma \Vert \nu _{z,\varepsilon }\Vert _\infty }{\lambda _K g}+2\varepsilon ^2\gamma \frac{{C''}}{\lambda _K}, \end{aligned}$$

  we get that, thanks to (41) and (39), for a given $0<\delta <\Delta$, there exists $\eta _\delta >0$ depending only on $\delta ,g,m,K,t^*,F,G,\Vert \nu _{z,\varepsilon }\Vert _\infty$ such that :

  $$\begin{aligned} \forall (\varepsilon ,|r^\varepsilon _Y(0)|,|r^\varepsilon _z(0)|) \in [0,\eta _\delta ]^3, \underset{t\le \min (T,t^*_\varepsilon )}{\max }|r^\varepsilon _Y(t)|+|r^\varepsilon _z(t)|\le \delta . \end{aligned}$$

  Recalling that $t^*_\varepsilon = \min \left( t^*,\inf \{t>0, |r^\varepsilon _z|+|r^\varepsilon _Y|> \Delta \}\right)$, we get that $T\le t^*_\varepsilon$, for $\delta <\Delta$ and $(\varepsilon ,|r^\varepsilon _Y(0)|,|r^\varepsilon _z(0)|) \in [0,\eta _\delta ]^3$. Consequently, the convergence is uniform on [0, T]. $\square$

Proof

(Proof of Theorem 3.1) One can notice that the control obtained in the proof of Proposition 1 can be applied on any time interval $[a,a+T]$ with $a\in [0,t^*-T]$, provided that $(\varepsilon ,|r^\varepsilon _Y(a)|,|r^\varepsilon _z(a)|)$ are small enough. Therefore, we can reiterate the control a finite number of times on the intervals $[jT,\min \{(j+1)T,t^*\}]$ with $\forall j \le \lfloor \frac{t^*}{T}\rfloor$. Hence, the uniform convergence on $[0,t^*]$. $\square$

Proof

(Proof of Lemma 11) Recall that for all $z \in K$, $J_z$ has real negative eigenvalues, uniformly bounded over K by $-2\lambda _K<-\lambda _K$. Let us define, for $z\in K$:

$$\begin{aligned} f_{\lambda _K,z} : \mathbb {R}_+ \rightarrow \mathbb {R}_+, \quad s \mapsto {\left| \left| \left| e^{J_z s}e^{\lambda _K s} \right| \right| \right| }. \end{aligned}$$

For all $z\in \mathbb {K}$, $f_{\lambda _K,z}$ is continuous. Moreover, Theorem 2.34 of Chicone (1999) ensures that $f_{l,z}$ is bounded for all $l<2\lambda _K$.

We can thus define :

$$\begin{aligned} \Gamma _{\lambda _K} : K \rightarrow \mathbb {R}^*_+,\quad z \mapsto \underset{s\ge 0}{\max }\,f_{\lambda _K,z}(s). \end{aligned}$$

Let us show that $\Gamma _{\lambda _K}$ is a continuous function. Let $z_0\in K$ and $\varepsilon >0$.

One can first notice that, for $s\ge 0$:

$$\begin{aligned} f_{\lambda _K,z}(s) = f_{\frac{3\lambda _K}{2},z}(s)e^{-\frac{\lambda _K}{2}s}<\Gamma _{\frac{3\lambda _K}{2},z}e^{-\frac{\lambda _K}{2}s}. \end{aligned}$$

Thus, $f_{\lambda _K,z}$ vanishes when s goes to infinity. As a consequence, there exists $s_0 \ge 0$ such that:

$$\begin{aligned} \Gamma _{\lambda _K}(z_0)={\left| \left| \left| e^{J_{z_0}s_0}e^{\lambda _K s_0} \right| \right| \right| }. \end{aligned}$$

Furthermore, for $l\in ]\lambda _K,2\lambda _K[$, we have:

$$\begin{aligned} \Gamma _{l}(z_0)={\left| \left| \left| e^{J_{z_0}s_0}e^{l s_0} \right| \right| \right| }=\Gamma _{\lambda _K}(z_0)e^{(l-\lambda _K)s_0}. \end{aligned}$$

We can therefore choose $l\in ]\lambda _K,2\lambda _K[$ such that $\Gamma _{\lambda _K}(z_0) \le \Gamma _{l}(z_0) \le \Gamma _{\lambda _K}(z_0)+\varepsilon$.

As $z\mapsto J_z$ is a continuous function, there exists $\delta >0$ that ensures that for if $z\in K$ and $|z-z_0|\le \delta$, then:

$$\begin{aligned} {\left| \left| \left| J_z-J_{z_0} \right| \right| \right| }<\frac{l-\lambda _K}{2\Gamma _l(z_0)}. \end{aligned}$$

Let us consider such a z.

As $e^{J_z s}$ is solution of the ODE : $y'=J_{z_0}y+(J_{z}-J_{z_0})y$, we obtain, for $s\ge 0$:

$$\begin{aligned} e^{J_z s} = e^{J_{z_0} s} + e^{J_{z_0} \cdot }*(J_z-J_{z_0})e^{J_{z}\cdot } (s). \end{aligned}$$

Hence :

$$\begin{aligned} {\left| \left| \left| e^{J_z t} \right| \right| \right| }\le \Gamma _l(z_0)e^{-ls}+ \frac{l-\lambda _K}{2}{\left| \left| \left| e^{J_z \cdot } \right| \right| \right| }*e^{-l\cdot } \end{aligned}$$

From applying Gronwall’s inequality to $t\mapsto {\left| \left| \left| e^{J_z s} \right| \right| \right| }e^{ls}$, it comes that, for $s\ge 0$:

$$\begin{aligned} {\left| \left| \left| e^{J_z s} \right| \right| \right| }\le & {} \Gamma _l(z_0)e^{-\left( l-\frac{l-\lambda _K}{2}\right) t} \le \Gamma _l(z_0)e^{-\left( \frac{l+\lambda _K}{2}\right) s} \\\le & {} \left[ \Gamma _{\lambda _K}(z_0)+\varepsilon \right] e^{-\lambda _K s}. \end{aligned}$$

Hence:

$$\begin{aligned} \Gamma _{\lambda _K}(z) \le \Gamma _{\lambda _K}(z_0)+\varepsilon . \end{aligned}$$

Moreover, recall that $t_0$ was defined so that :

$$\begin{aligned} \Gamma _{\lambda _K}(z_0)={\left| \left| \left| e^{J_{z_0}s_0}e^{\lambda _K s_0} \right| \right| \right| }. \end{aligned}$$

Then, by continuity of $z\mapsto e^{J_z s_0}$, there exists $\delta '>0$ that ensures that for $|z-z_0|\le \delta '$, we have:

$$\begin{aligned} {\left| \left| \left| e^{J_{z}s_0}e^{\lambda _K s_0} \right| \right| \right| } \ge {\left| \left| \left| e^{J_{z_0}s_0}e^{\lambda _K s_0} \right| \right| \right| }-\varepsilon . \end{aligned}$$

Hence:

$$\begin{aligned} \Gamma _{\lambda _K}(z) \ge \Gamma _{\lambda _K}(z_0) - \varepsilon . \end{aligned}$$

In conclusion, if $|z-z_0| \le \min (\delta ,\delta ')$, then $|\Gamma _{\lambda _K}(z)-\Gamma _{\lambda _K}(z_0)| \le \varepsilon$. Hence $\Gamma _{\lambda _K}$ is continuous over K. Furthermore, as K is a compact set, $\Gamma _{\lambda _K}$ is bounded, by $\gamma$. $\square$

F Proof of Proposition 3.1

This appendix is dedicated to the proof of Proposition 3.1.

Proof

Let $(g,m,z^*)\in \mathbb {R}^*_+\times \mathbb {R}^*_+\times \mathbb {R}_+$ be such that $P_{z^*}$ has a single positive root. From Lemma 1, this root defines a fast equilibrium if it is greater than $f_1(z^*)$. From Lemma 2, that is the case if and only if $f_1(z^*)$ is negative or $P_{z^*}(f_1(z^*))$ is negative.

First, regarding the sign of $f_1(z^*)$, we have:

$$\begin{aligned} f_1(z^*)<0 \iff (z^*+1)^2 < \frac{1-m}{g}, \end{aligned}$$

which requires that $m<1$. If $m<1$ then:

$$\begin{aligned} f_1(z^*)<0 \iff 0\le z^* < \sqrt{\frac{1-m}{g}}-1, \end{aligned}$$

which requires that $m+g< 1$. Hence:

$$\begin{aligned} f_1(z^*)<0 \iff [m+g<1]\wedge [z^* < \sqrt{\frac{1-m}{g}}-1]. \end{aligned}$$

Next, regarding the sign of $P_{z^*}(f_1(z^*))$, we compute:

$$\begin{aligned} P_{z^*}(f_1(z^*))&= f_1(z^*)f_2(z^*)-1\\&= \left( 1+\frac{g}{m}(z^*+1)^2-\frac{1}{m}\right) \left( 1+\frac{g}{m}(z^*-1)^2-\frac{1}{m}\right) -1\\&= \frac{g^2}{m^2} \left[ {z^*}^4+{z^*}^2\,\frac{2(m-g-1)}{g}+\frac{(g-1)(2m+g-1)}{g^2}\right] \end{aligned}$$

Let us define:

$$\begin{aligned} Q(X) = X^2+X\,\frac{2(m-g-1)}{g}+\frac{(g-1)(2m+g-1)}{g^2}, \end{aligned}$$

$z_1,z_2$ its two roots and $\Delta = \frac{4}{g^2}\left[ m^2-4g\,(m-1)\right]$ its discriminant. From the computation above,

$$\begin{aligned} P_{z^*}(f_1(z^*)) <0 \iff [\;\Delta >0\;]\wedge \left[ \;{z^*}^2 \in ]z_1,z_2[\;\right] . \end{aligned}$$

We have:

$$\begin{aligned} \Delta>0=&\iff m^2-4\,g\,m+4\,g>0 \\&\iff [g<1]\vee \left[ [g\ge 1] \wedge \left[ \left[ 0<m<2g\left( 1-\sqrt{1-\frac{1}{g}}\right) \right] \right. \right. \\&\left. \left. \vee \left[ m>2g\left( 1+\sqrt{1-\frac{1}{g}}\right) \right] \right] \right] \end{aligned}$$

and:

$$\begin{aligned} z_1z_2 = \frac{(g-1)(2m+g-1)}{g^2},\quad z_1+z_2 =\frac{2(g+1-m)}{g}. \end{aligned}$$

Consequently:

$\diamond$:

if $g\ge 1$, then $2m+g-1>0$ and then $z_1z_2\ge 0$. If, additionally, $m<2g\left( 1-\sqrt{1-\frac{1}{g}}\right)$, then $m< 2 \le g+1$ ($g\mapsto 2g-2\sqrt{g^2-g}$ is decreasing on $[1,+\infty [$). Therefore, we get: $z_1+z_2> 0$ and thus, $z_2>0$ and $z_1\ge 0$. At last, if $m>2g\left( 1+\sqrt{1-\frac{1}{g}}\right)$, then $m>2\,g\ge g+1$, which implies $z_1+z_2<0$ and thus $z_1< 0, z_2\le 0$.

$\diamond$:

if $g< 1$, then $z_1+z_2\ge 0$ if and only if $m\le g+1$ and $z_1z_2\ge 0$ if and only if $m\le \frac{1-g}{2}$ (which is lower than $g+1$).

Hence the result. $\square$

G Proof of Lemma 9

This section is dedicated to proving Lemma 9, which concludes the proof of Proposition 4.2.

Proof

(Proof of Lemma 9)

Let $(m,g)\in {\mathbb {R}_+^*}^2$ verify (25). Then, from the first part of the proof of Proposition 4.2, there exists a unique $\rho ^*>0$ that is solution of the equation in (23). Let us define $N_1^*$ and $N_2^*$ such as in (26). Then we have: $0<\rho ^* = \frac{N_2^*}{N_1^*}$. Thus:

$$\begin{aligned} N_1^*>0 \iff N_2^*>0 \iff \frac{1}{m}(N_1^*+N_2^*) > 0. \end{aligned}$$

Borrowing once again the notations: $a = \frac{4g}{m}$, $b=\frac{1}{m}$ and $y^*=\rho ^*+\frac{1}{\rho ^*}$ (unique root of S larger than 2), (26) leads to:

$$\begin{aligned} \frac{1}{m}(N_1^*+N_2^*)= & {} 2\left( \frac{1}{m}-1\right) + {y^*} -\frac{4g}{m}\frac{{y^*}^2-2}{{y^*}^2}\\= & {} \frac{1}{{y^*}^2}\left[ {y^*}^3+\left[ \frac{1-2m}{m}+\frac{1}{m}-\frac{4g}{m}\right] {y^*}^2+2\times \frac{4g}{m}\right] \\= & {} \frac{1}{{{y^*}^2}}\left[ S({y^*}) + (\frac{1-2m}{m}){y^*}^2+\frac{4g}{m}{y^*}+\frac{4g}{m}\right] . \end{aligned}$$

As $S({y^*}) = 0$, we get:

$$\begin{aligned} N_1^*>0\iff N_2^*>0\iff (1-2m){y^*}^2+4g{y^*}+4g>0. \end{aligned}$$

This is always true whenever $m \le \frac{1}{2}$. Otherwise, let us suppose henceforth that $2m>1$. The condition above is equivalent to:

$$\begin{aligned} {y^*} < c+\sqrt{c^2+2c},\quad \text { where: }c = \frac{2g}{2m-1} >0. \end{aligned}$$

Let us show that: $c+\sqrt{c^2+2c}\ge 2$. It is sufficient to show that: $c\ge \frac{2}{3}$, which is equivalent to having: $3g+1\ge 2m$. In this proof, we are considering $(m,g) \in {\mathbb {R}_+^*}^2$ such that $1+2m<5g$ and $4g\,(m-1)<m^2$. Let us show that such pairs verify $3g+1\ge 2m$:

$\diamond$:

if $g\le 1$, then $m < \frac{5g-1}{2} \le \frac{3g+1}{2}$.

$\diamond$:

if $g\ge 1$, then $m<2g-2\sqrt{g^2-g}$ which is a decreasing function on $[1,+\infty [$, which takes the value 2 when $g=1$. Hence it is always dominated by $g\mapsto \frac{3g+1}{2}$ on this interval.

Hence $c+\sqrt{c^2+2c}\ge \frac{2}{3}+\sqrt{\frac{4}{9}+\frac{4}{3}}=2$. Therefore, as $y^*$ is the only root of S greater than 2, we get the following equivalence:

$$\begin{aligned} y^*<c+\sqrt{c^2+c}\iff S\left( c+\sqrt{c^2+2c}\right) > 0. \end{aligned}$$

The rest of the proof is dedicated to examine the conditions on (m, g) under which:

$$\begin{aligned} S\left( c+\sqrt{c^2+2c}\right) > 0. \end{aligned}$$

Let us set $Q:=\sqrt{c^2+2c} = \sqrt{4g\,\frac{g+2m-1}{(2m-1)^2}}$. Tedious computations done with the help of Mathematica show that: $S(c) = Q^2\left[ \frac{g(4-6m)+(2m-1)^2}{{m}\,(2m-1)}\right]$, and we next compute:

$$\begin{aligned} S(c+Q)= & {} S(c)+Q^2\left[ 3c+\frac{1-4g}{m}\right] +Q\left[ Q^2+3c^2+2c\,\frac{(1-4g)}{m}-\frac{4g}{m}\right] \\ {}= & {} Q^2\left[ \frac{g(4-6m)+(2m-1)^2}{{m}\,(2m-1)}+\frac{6g}{2m-1}+\frac{1-4g}{m}\right] \\&+Q\left[ 4c^2+2c\,\frac{(m+1-4g)}{m}-\frac{4g}{m}\right] \\= & {} Q\left[ 2Q\frac{\left( 2m^2-m-4g\,(m-1)\right) }{m (2 m-1)}-\frac{4 g \left( 4 g\,(m-1)+2 m^2-5 m+2\right) }{m(2m-1)^2}\right] . \end{aligned}$$

Hence:

$$\begin{aligned}&S(c+Q)>0\\&\iff Q\left( 2m^2-m-4g\,(m-1)\right)> 2g\frac{\left( 4 g\,(m-1)+2 m^2-5 m+2\right) }{(2m-1)}\\&\iff \sqrt{g+2m-1}\left( 2m^2-m-4g\,(m-1)\right) \\&\quad >2\sqrt{g}\left( 4 g\,(m-1)+2 m^2-5 m+2\right) . \end{aligned}$$

Let us study different cases corresponding to different ranges of value of $m>\frac{1}{2}$.

If $\underline{m=1}$, then the last line is equivalent to :

$$\begin{aligned} \sqrt{1+g}>-2\sqrt{g}, \end{aligned}$$

which is true for all $g>0$.

If $\underline{\frac{1}{2}<m<1}$, then:

$$\begin{aligned} 4g\,(m-1) + 2m^2-5m+2 = 4g\,(m-1) +2(m-2)\left( m-\frac{1}{2}\right) <0, \end{aligned}$$

and:

$$\begin{aligned} 2m^2-m-4g\,(m-1) = 2m\left( m-\frac{1}{2}\right) + 4g(1-m)>0. \end{aligned}$$

Hence, for all g such that $1+2m<5g$ and $m^2>4g\,(m-1)$:

$$\begin{aligned} \sqrt{g+2m-1}\left( 2m^2-m-4g\,(m-1)\right) >2\sqrt{g}\left( 4 g\,(m-1)+2 m^2-5 m+2\right) . \end{aligned}$$

If $\underline{m>1}$, then:

$$\begin{aligned} 2m^2-m>m^2>4g\,(m-1). \end{aligned}$$

Hence, if: $4g\,(m-1)+2 m^2-5 m+2<0$, then, for all g such that $1+2m<5g$ and $m^2>4g\,(m-1)$:

$$\begin{aligned} \sqrt{g+2m-1}\left( 2m^2-m-4g\,(m-1)\right) >2\sqrt{g}\left( 4 g\,(m-1)+2 m^2-5 m+2\right) . \end{aligned}$$

Otherwise, if $4g\,(m-1)+2 m^2-5 m+2\ge 0$, then:

$$\begin{aligned}&S(c+Q)>0\\&\iff \sqrt{g+2m-1}\left( 2m^2-m-4g\,(m-1)\right)>2\sqrt{g}\left( 4 g\,(m-1)+2 m^2-5 m+2\right) \\&\iff \left( 1+\frac{2m-1}{g}\right) \left( 2m^2-2-4g\,(m-1)\right) ^2 > 4\left( 4g\,(m-1)+2m^2-5m+2\right) ^2. \end{aligned}$$

Let us note $x:=\frac{2m-1}{g}$. Then, the latter is equivalent to:

$$\begin{aligned}&(1+x)\left[ (m-1)x + \left( x-4(m-1)\right) \right] ^2 - \left[ (m-1)x - \left( x-4(m-1)\right) \right] ^2>0\\&\iff 4(m-1)x(x-4(m-1)) + x\left( mx-4(m-1)\right) ^2>0\\&\iff 4(m-1)x-16(m-1)^2+m^2x^2-8mx(m-1)+16(m-1)^2>0\\&\iff m^2x^2+4x(m-1)(1-2m)>0\\&\iff m^2x^2-4x^2g\,(m-1)>0\\&\iff m^2>4g\,(m-1). \end{aligned}$$

$\square$

H Details of the numerical analysis carried out in Sect. 2

Domains We consider a bounded trait domain $[-\varvec{z_{\max }},\varvec{z_{\max }}]$, discretised in a mesh $\left( \varvec{z}_k\right) _{0\le k < K}$ (K odd) with regard to the step length $\delta z>0$, and a time domain $[0,\varvec{T_{\max }}]$, discretised in a mesh $\left( \varvec{t}^l\right) _{0\le l< L}$ with regard to the step length $\delta t>0$. In the the simulations involved in Fig. 2, we use the following values for the parameters:

$$\begin{aligned} \varvec{z_{\max }} = 7,\quad \varvec{T_{\max }} = 1000, \quad \delta z = 1.6\times 10^{-2}, \quad \delta t = 5\times 10^{-3}. \end{aligned}$$

Scheme For $i\in \{1,2\}, 0\le l <L$, we approximate the trait distributions $\varvec{n_i}(\varvec{t}^l,\cdot )$ by $\left( \varvec{n}_{i,k}^l\right) _{0\le k<K}$ with the following semi-implicit scheme:

$$\begin{aligned} \frac{\varvec{\sigma }^2}{\delta t}\left( \varvec{n}_{i,k}^{l+1} -\varvec{n}_{i,k}^{l} \right) = \varvec{r}\,\varvec{B}_{i,k}^l -\left( \varvec{g}\,\left( \varvec{z_k}-\varvec{\theta _i}\right) ^2 +\varvec{\kappa }\,\varvec{N}_i^l +\varvec{m}\right) \,\varvec{n}_{i,k}^{l+1} + \varvec{m}\,\varvec{n}_{j,k}^{l+1}, \end{aligned}$$

where $\varvec{N}_i^{l} = \sum _{k=0}^{K-2} \varvec{n}_{i,k}^l\,\delta z$ and $\varvec{B}_{i,k}^l$ is a discretisation of the reproduction operator $\varvec{\mathcal {B}_\sigma }(\varvec{n_i}(t^l,z_k)$. In the next paragraph, we detail how we compute $\left( \varvec{B}_{i,k}^l\right) _{0\le k<K}$.

We approximate the system of moments of Ronce and Kirkpatrick (2001) following a similar semi-implicit scheme.

Discretization of the reproduction operator The discretization of the reproduction operator is in accordance with the double convolution form shown in Lemma 10, as it increases greatly the computational speed in comparison to a double loop. However, the half-arguments involved in Lemma 10 calls for a special attention to the meshes involved.

Let us define two auxiliary trait meshes

• 1. $\left( \varvec{\tilde{z}}_{k'}\right) _{0\le k' < 2K-1}$ on $[-\varvec{z_{\max }},\varvec{z_{\max }}]$, with step length $\frac{\delta z}{2}$,

• 2. $\left( \varvec{\hat{z}}_{k''}\right) _{0\le k'' < 4K-3}$ on $[-2\varvec{z_{\max }},2\varvec{z_{\max }}]$, with step length $\frac{\delta z}{2}$.

We define the vector $\left( G_{k'}\right) _{0\le k' < 2K-1}$ discretising the Gaussian kernel involved in our reproduction operator on the trait grid $\left( \varvec{\tilde{z}}_{k'}\right) _{0\le k' < 2K-1}$:

$$\begin{aligned} G_{k'} = \frac{1}{\sqrt{\pi }\varvec{\sigma }}\exp \left[ -\frac{\varvec{\tilde{z}}_{k'}^2}{\varvec{\sigma ^2}}\right] . \end{aligned}$$

We next define the vector $\left( \hat{B}_{i,k''}^l\right) _{0\le k''<4K-3}$ resulting from the following double discrete convolution (denoted $*$):

$$\begin{aligned} \left( \hat{B}_{i,k''}^l\right) _{0\le k''<4K-3}&=\frac{1}{\varvec{N}_i^l} \left( \varvec{n}_{i,k}^l\right) _{0\le k<K} *\left( \varvec{n}_{i,k}^l\right) _{0\le k<K} *\left( G_{k'}\right) _{0\le k' < 2K-1} \end{aligned}$$

We use a convolution algorithm with default settings: the size of the output is the sum of entry vector sizes minus one, and out of bounds index entries are extrapolated as 0. A straight-forward computation shows that $\left( \hat{B}_{i,k''}^l\right) _{0\le k''<4K-3}$ is the approximation of the reproduction operator on the mesh $\left( \varvec{\hat{z}}_{k''}\right) _{0\le k'' < 4K-3}$:

$$\begin{aligned} \hat{B}_{i,k''}^l&=\frac{\delta z^2}{\varvec{N}_i^l}\sum _{k_1=0}^{4K-4}\varvec{n}_{i,k_1}^l\sum _{k_2=0}^{3K-2}\varvec{n}_{i,k_2}^l\,G_{k''-k_1-k_2}\\&= \frac{\delta z^2}{\varvec{N}_i^l} \sum _{k_1=0}^{4K-4}\varvec{n}_{i,k_1}^l\sum _{k_2=0}^{3K-2}\varvec{n}_{i,k_2}^l\,\frac{1}{\sqrt{\pi }\varvec{\sigma }}\exp \left[ -\frac{\varvec{\tilde{z}}_{k''-k_1-k_2}^2}{\varvec{\sigma ^2}}\right] \\&= \frac{\delta z^2}{\varvec{N}_i^l}\sum _{k_1=0}^{4K-4}\varvec{n}_{i,k_1}^l\sum _{k_2=0}^{3K-2}\varvec{n}_{i,k_2}^l\,\frac{1}{\sqrt{\pi }\varvec{\sigma }}\exp \left[ -\frac{\left( -\varvec{z}_{\max }+\frac{\delta z}{2} (k''-k_1-k_2)\right) ^2}{\varvec{\sigma ^2}}\right] \\&= \frac{\delta z^2}{\varvec{N}_i^l}\sum _{k_1=0}^{4K-4}\varvec{n}_{i,k_1}^l\sum _{k_2=0}^{3K-2}\varvec{n}_{i,k_2}^l\,\frac{1}{\sqrt{\pi }\varvec{\sigma }}\\&\quad \exp \left[ -\frac{\left( -2\varvec{z}_{\max }+\frac{k''\delta z}{2} -\frac{(-\varvec{z}_{\max }+k_1\delta z)+(-\varvec{z}_{\max }+k_2\delta z)}{2}\right) ^2}{\varvec{\sigma ^2}}\right] \\&= \frac{\delta z^2}{\varvec{N}_i^l}\sum _{k_1=0}^{4K-4}\varvec{n}_{i,k_1}^l\sum _{k_2=0}^{3K-2}\varvec{n}_{i,k_2}^l\,\frac{1}{\sqrt{\pi }\varvec{\sigma }}\exp \left[ -\frac{\left( \varvec{\hat{z}}_{k''} -\frac{\varvec{z}_{k_1}+\varvec{z}_{k_2}}{2}\right) ^2}{\varvec{\sigma ^2}}\right] . \end{aligned}$$

Thus, we interpolate $\left( \hat{B}^l_{i,k''}\right) _{0\le k''<4K-3}$ at the entries corresponding to $\left( \varvec{z}_k\right) _{0\le k < K}$ to obtain $\left( \varvec{B}_{i,k}^l\right) _{0\le k<K}$.

I Numerical outcomes details: Figs. 6 and 7

Numerical setting The lower panel of Fig. 6 has been produced by running 3600 simulations, one for each couple of migration rate $m \in [0.01,\,3]$ and intensity of selection $g\in [0.01,\,3]$, for $t\le T_{max} \in \left[ \frac{300}{\varepsilon ^2},\frac{600}{\varepsilon ^2}\right]$, with a criteria to cut the simulation short at a time greater than $\frac{300}{\varepsilon ^2}$ if the difference between two consecutive steps is small enough. The value of the other parameters are the same for each simulation: $r = 1, \theta = 1, \kappa = 1, \varepsilon = 0.05$, as well as the initial state:

$$\begin{aligned}\begin{aligned} {\left\{ \begin{array}{ll} n_1^0(z) = 0.99\,\times \,\frac{e^{-\frac{(z-0.2)^2}{2\varepsilon ^2}}}{\sqrt{2\pi }\varepsilon },\\ n_2^0(z) = \frac{e^{-\frac{(z-0.2)^2}{2\varepsilon ^2}}}{\sqrt{2\pi }\varepsilon }. \end{array}\right. } \end{aligned} \end{aligned}$$

The initial state is taken as monomorphic, as the aim of this figure is to be compared to the theoretical outcomes that are predicted within the scope of the slow-fast analysis as stated in Theorem 3.1 (so when the initial state is close enough from the slow manifold).

Scoring Each simulation final state $(n_1^f,n_2^f)$ is attributed a score between 0 and 1 according to the following scheme:

• 1. If $\max \left( N_1^f,N_2^f\right) <0.01$, then the score is 0 (for extinction) and is corresponding to the deep purple color. Else, the score is a positive number (lower than 1) according to what follows.

• 2. If the variance in trait of the metapopulation is greater than $2\,\varepsilon ^2$, the score is 1 (corresponding to the color yellow). This would be the case if the final state is dimorphic, but more generally, this is to highlight the simulations whose final state does not fall in the small segregational variance regime analysis prediction (which in particular predicts that the distribution of trait in the metapopulation is monomorphic (see Sect. 3), with a variance of order $\varepsilon ^2$ (see (12)).

• 3. If both conditions above are not met, then the score S is given according to the following formula:

  $$\begin{aligned} S = \frac{5}{6} - \frac{1}{3}\frac{\left| N_2^f-N_1^f\right| }{N_1^f+N_2^f}. \end{aligned}$$

  This formula discriminates between symmetrical equilibria (which are characterized by equal population sizes, see Proposition 4.1), which typically have a score of $\frac{5}{6}$ (corresponding to the color light green), and asymmetrical equilibria, which have a discrepancy in local population sizes and therefore have a typically much lower score (in the blue tones).

Adjustments for Fig. 7 The methodology is the same for the lower panel of Fig. 6 and both panels of Fig. 7, at the exception of the initial state, set as:

$$\begin{aligned}\begin{aligned} {\left\{ \begin{array}{ll} n_1^0(z) = 0.9\,\times \,\frac{e^{-\frac{(z+1)^2}{2\varepsilon ^2}}}{\sqrt{2\pi }\varepsilon },\\ n_2^0(z) = \frac{e^{-\frac{(z-1)^2}{2\varepsilon ^2}}}{\sqrt{2\pi }\varepsilon }. \end{array}\right. } \end{aligned}. \end{aligned}$$

and of the time step for the lower panel of Fig. 7, which is refined to keep up with the smaller value of $\varepsilon ^2$.