# Impacts of Tempo and Mode of Environmental Fluctuations on Population Growth: Slow- and Fast-Limit Approximations of Lyapunov Exponents for Periodic and Random Environments

Pierre Monmarché, Sebastian J. Schreiber, Édouard Strickler

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

## To cite this version:

Pierre Monmarché, Sebastian J. Schreiber, Édouard Strickler. Impacts of Tempo and Mode of Environmental Fluctuations on Population Growth: Slow- and Fast-Limit Approximations of Lyapunov Exponents for Periodic and Random Environments. 2024. hal-04769603

HAL Id: hal-04769603

https://hal.science/hal-04769603v1

Preprint submitted on 6 Nov 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Impacts of Tempo and Mode of Environmental Fluctuations on Population Growth: Slow- and Fast-Limit Approximations of Lyapunov Exponents for Periodic and Random Environments

Pierre Monmarché¹,², Sebastian J. Schreiber³, and Édouard Strickler⁴

¹Sorbonne Université, 4 place Jussieu 75005 Paris (France)

²Institut universitaire de France (IUF)

³Department of Evolution and Ecology and Center for Population Biology, University of California, Davis 95616 USA

⁴Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy (France)

August 22, 2024

## Abstract

Populations consist of individuals living in different states and experiencing temporally varying environmental conditions. Individuals may differ in their geographic location, stage of development (e.g. juvenile versus adult), or physiological state (infected or susceptible). Environmental conditions may vary due to abiotic (e.g. temperature) or biotic (e.g. resource availability) factors. As survival, growth, and reproduction of individuals depend on their state and the environmental conditions, environmental fluctuations often impact population growth. Here, we examine to what extent the tempo and mode (i.e. periodic versus random) of these fluctuations matter for population growth. To this end, we model population growth for a population with *d* individual states and experiencing *N* different environmental states. The models are switching, linear ordinary differential equations $x'(t) = A(\sigma(\omega t))x(t)$ where $x(t) = (x_1(t),..., x_d(t))$ corresponds to the population densities in the *d* individual states, $\sigma(t)$ is a piece-wise constant function representing the fluctuations in the environmental states $1,..., N$, $\omega$ is the frequency of the environmental fluctuations, and $A(1),..., A(n)$ are Metzler matrices representing the population dynamics in the environmental states $1,..., N$. $\sigma(t)$ can either be a periodic function or correspond to a continuous-time Markov chain. Under suitable conditions, there exists a Lyapunov exponent $\Lambda(\omega)$ such that $\lim_{t\to\infty} \frac{1}{t} \log \sum_i x_i(t) = \Lambda(\omega)$ for all non-negative, non-zero initial conditions $x(0)$ (with probability one in the random case). For both random and periodic switching, we derive analytical first-order and second-order approximations of $\Lambda(\omega)$ in the limits of slow ($\omega \to 0$) and fast ($\omega \to \infty$) environmental fluctuations. When the order of switching and the average switching times are equal, we show that the first-order approximations of $\Lambda(\omega)$ are equivalent in the slow-switching limit, but not in the fast-switching limit. Hence, the mode (random versus periodic) of switching matters for population growth. We illustrate our results with applications to a simple stage-structured model and a general spatially structured model. When dispersal rates are symmetric, the first order approximations suggest that population growth rates increase with the frequency of switching – consistent with earlier work on periodic switching. In the absence of dispersal symmetry, we demonstrate that $\Lambda(\omega)$ can be non-monotonic in $\omega$. In conclusion, our results show that population growth rates often depends both on the tempo ($\omega$) and mode (random versus deterministic) of the environmental fluctuations.
# 1 Introduction

Central to many questions in ecology, evolution, and epidemiology is identifying how environmental conditions determine population growth. This growth depends on the survival, growth, and reproduction of individuals, which in turn depends on the state of an individual (i-state). Individuals may differ in their i-state based on their geographic location, their stage of development, their genotype, or their behavior [DeAngelis, 2018]. One way to account for these difference in i-state on the population state (p-state) is using Matrix Population Models (MPM) [Caswell, 2001, DeAngelis, 2018]. In these models, there are a finite number d of i-states and the p-state corresponds to the vector, $x = (x_1, x_2, ..., x_d)$, of population densities. For continuous-time MPMs, a matrix A with non-negative off diagonal elements, $a_{ij}$ with $i \neq j$, define the rates at which individuals in state j contribute to individuals in different state j. These contributions may correspond to individuals transition to the other state (e.g. dispersing from one geographic location to another, developing from one stage to another, changing their behavior) or producing individuals of the other state. Diagonal elements, $a_{ii}$, of the matrix A can be positive, zero, or non-negative and correspond to the net rate at which individuals leave their current state (including death) and produce more individual of the same state. Under constant environmental conditions, the continuous-time MPM is a linear system of differential equations

$$x'(t) = Ax(t) \quad (1)$$

where A, with its non-negative off-diagonal entries, is called a Metzler matrix [Mitkowski, 2008]. These continuous-time MPM commonly arise when linearizing a disease-free equilibrium for an epidemiological model, the extinction equilibrium of a demographic model, or a fixation equilibrium of a population genetics model [Bürger, 2000, Hethcote, 2000, Kon et al., 2004]. When the matrix A is irreducible, the Perron-Frobenius Theorem implies that for any non-negative, non-zero initial condition $x(0)$

$$\lim_{t \to \infty} \frac{1}{t} \ln \sum_{i} x_{i}(t) = \lambda(A) \quad (2)$$

where $\lambda(A)$ is the spectral abscissa of A i.e. the maximum of the real parts of the eigenvalues of A. Hence, $\lambda(A)$ determine the long-term growth rate of the population. In particular, when $\lambda(A)$ is positive, the population persists and grows. When $\lambda(A)$ is negative, the population declines exponentially quickly to the extinction state.

Most populations experience fluctuations in environmental conditions such as temperature, precipitation, and nutrient availability. As these environmental conditions determine survival, growth, and reproduction of individuals, environmental fluctuations lead to fluctuations in demographic rates. Understanding the implication of these fluctuations on population growth is an active topic [Tuljapurkar, 1990, Boyce et al., 2006, Fay et al., 2020, Hilde et al., 2020, Kortessis et al., 2023b, Paniw et al., 2018, Schreiber, 2023, Tuljapurkar and Haridas, 2006]. Nearly all of this earlier work, see however [Benaïm and Strickler, 2019, Benaïm et al., 2023, Benaim et al., 2024], is for discrete-time MPMs. One of the simplest approaches to modeling this time-dependency in continuous-time MPMs is to allow A in (1) to switch between a finite number of Metzler matrices $A_1, A_2, ..., A_N$ corresponding to N environmental states e.g. summer versus winter, wet versus dry period. If $\sigma(t)$ denotes the environmental state at time t, then we get a switched, continuous-time MPM

$$x'(t) = A_{\sigma(t)}x(t).$$

If $\sigma(t)$ is a periodic piecewise constant function or given by an irreducible Markov chain, then there exists the dominant Lyapunov exponent $\Lambda$ such that

$$\Lambda = \lim_{t \to \infty} \frac{1}{t} \ln \sum_{i} x_{i}(t) \quad (3)$$
holds (with probability one in the random case) for any non-negative, non-zero initial condition $x(0)$ (see Lemmas 1 and 2 below). As in the constant environment case (1), the Lyapunov exponent measures the long-term population growth rate of the system. In particular, by a linearization procedure, it allows to determine the persistence or extinction for switched non-linear flows [Benaïm and Strickler, 2019].

Environmental fluctuations may occur at daily, yearly, and multi-decadal time scales. This raises the question to what extent does the rate of switching between environmental states influence the population growth rate $\Lambda$? To tackle this question, consider a fixed choice of the Meltzer matrices $A_1, \dots, A_N$ and environmental trajectory $\sigma(t)$. To manipulate the rate of switching, let $\omega > 0$ be the frequency of environmental switching in which case our MPM becomes

$$x'(t) = A_{\sigma(\omega t)}x(t). \qquad (4)$$

Let $\Lambda(\omega)$ denote the dominant Lyapunov exponent (3) associated with the MPM with frequency $\omega$. When switching is fast ($\omega \to \infty$), solutions of the switching model (4) converges to the solution of an ODE $x' = \bar{A}x$ where the averaged matrix $\bar{A}$ corresponds to the appropriate convex combination of the matrices $A_i$. In this limit, $\Lambda(\omega)$ converges to $\lambda(\bar{A})$ [Benaïm et al., 2014, Benaïm et al., 2023, Benaim et al., 2024, Benaïm and Strickler, 2019, Chitour et al., 2021, Du et al., 2021]. On the contrary, when the switching is slow ($\omega \to 0$), the Lyapunov exponent converge to a weighted average of the Lyapunov exponents $\lambda(A_i)$ of the matrices $A_i$ [Benaim et al., 2023, Benaim et al., 2024]. As the Lyapunov exponents of these limits, in general, are different, these earlier results imply that whether populations grow or decline can depend on the frequency of environmental switching.

These earlier results raise two interesting questions: Do the population growth rates increase or decrease as one approaches these limits? To what extent does the answer depend on whether the environmental switching is periodic or random? Whether periodic switching or random switching between some matrices yields the same behavior is a recurrent question in the study of switched linear systems. For instance, Chitour et al. [2021] provides conditions on the matrices $A_1, \dots, A_N$ under which the maximal Lyapunov exponent obtained by considering the worst deterministic signal $\sigma$ is strictly greater than the maximal probabilistic Lyapunov exponent obtained by considering the worst Markov chain. In particular, this shows that there exists matrices $A_1, \dots, A_N$ such that a suitable deterministic switching signal makes the system (4) unstable (the population persists), while it is stable for all Markovian signals (the population goes extinct). For applications in population dynamics, one is often interested either in the sign of the Lyapunov exponent (i.e. growth or extinction) or in its monotonicity with respect to some key parameters. For specific models, these questions can be answered in the fast or slow switching regimes with the explicit formula obtained in the present work. In particular, the monotonicity of the population growth in the switching rate is explicitly given in these two regimes by the sign of the first order term in its expansion.

To address these questions, let $\Lambda_p(\omega)$ and $\Lambda_M(\omega)$ denote the Lyapunov exponents in the periodic case and the Markovian case, respectively. For the fast Markovian (random) switching, Monmarché and Strickler [2023] answered one of these questions by deriving an asymptotic expansion of $\Lambda_M(\omega)$ in terms of $1/\omega$. Here, we address the remaining questions by deriving an asymptotic expansion of $\Lambda_p(\omega)$ in terms of $1/\omega$, and deriving first order expansions in terms of $\omega$ of $\Lambda_p(\omega)$ and $\Lambda_M(\omega)$ is the slow-switching limit. For slow-switching, we show that the periodic and random approximations of $\Lambda(\omega)$ are equivalent. However, in the limit of fast-switching, these approximation are not equivalent. Indeed, when there is fast switching between two environments, we show the first-order correction in terms of $1/\omega$ always vanishes and derive a second-order correction in terms of $1/\omega^2$. In contrast, the first order term for Markovian switching is non-zero in general. We present these results in Section 3 after providing a more detailed description of the periodic and Markovian models in Section 2. To illustrate the applicability of our approximations, we apply them to a simple stage-structure model and explore their implications for spatially structure populations in Section 4. Consistent with earlier work of Katriel [2022] for periodic environments, our approximations suggest that population growth rates increase with switching
rates for spatially structured populations with symmetric dispersal. However, we use these approximations to illustrate how asymmetric dispersal can result in a non-monotonic relationship between switching frequency and population growth rates. Proofs of the main results are in Sections 6 and 8.

## 2 Lyapunov exponents of periodic and random switching models

Here, we provide more details about the switching continuous-time matrix population model in (4) and their Lyapunov exponents $\Lambda(\omega)$ in (3). For these models, $N \ge 1$ corresponds to the number of environmental states, $\llbracket 1, N \rrbracket = \{1, 2, \dots, N\}$ is the set of environmental states, and the demographic matrices $A_1, A_2, \dots, A_N$ are $d \times d$ Metzler matrices i.e. have non-negative off-diagonal entries. We consider two approaches for specifying the piece-wise constant function $\sigma(t)$ that takes values $\llbracket 1, N \rrbracket$. One approach is to assume that $\sigma(t)$ is a periodic, piecewise constant, right-continuous function. The other approach is to assume that $\sigma(t)$ corresponds to a continuous-time Markov chain on $\llbracket 1, N \rrbracket$. These approaches are described in detail below.

The state space for the matrix population model is the non-negative cone $\mathbb{R}_+^d = [0, \infty)^d$ of $\mathbb{R}^d$. If $A_i$ is Metzler for all $i$ in the model (4), then $x(0) \in \mathbb{R}_+^d$ implies that $x(t) \in \mathbb{R}_+^d$ for all $t \ge 0$ i.e. $\mathbb{R}_+^d$ is forward invariant for the dynamics.

For each form of the model, we also present results ensuring the existence of dominant Lyapunov exponent $\Lambda(\omega)$. To do this for the periodic case, let $r(A)$ denote the spectral radius of $A$ (i.e. the eigenvalue of $A$ with the largest absolute value) and $\lambda(A)$ denote the spectral abscissa of $A$ (i.e. the eigenvalue of $A$ with largest real part). For a Metzler and irreducible matrix $A$, we can apply Perron-Frobenius Theorem to $A + rI$ for some large enough $r$, to conclude that $\lambda(A)$ is a simple eigenvalue of $A$. Associated with $A$ is a unique positive right eigenvector $x$ with $\sum_i x_i = 1$ and a unique left eigenvector $y$ with $\sum_i x_i y_i = 1$.

### 2.1 The periodic case

For the periodic version of the model (3), we assume that the function $\sigma: \mathbb{R} \to \llbracket 1, N \rrbracket$ is piece-wise constant, right-continuous, has period 1 and successively takes the values $1, \dots, N$ as follows

$$ \sigma(t) = i \text{ for all } t \in [\tau_i, \tau_{i+1}), $$

where $\tau_0 = 0$ and $\tau_{i+1} = \tau_i + \alpha_i$ for some $\alpha_i \ge 0$ satisfying $\sum_{i=1}^N \alpha_i = 1$. $\alpha_i$ represents the proportion of time spent in environmental state $i$ during one period of length $T := \omega^{-1}$.

Let $M(\omega)$ be the *monodromy matrix* associated with model over the time interval $[0, 1]$:

$$ M(\omega) = e^{\omega^{-1} \alpha_N A_N} \dots e^{\omega^{-1} \alpha_1 A_1} $$

where $e^M$ denotes the matrix exponential of the matrix $M$. For all initial condition $x(0)$ of (4), one can verify that $x(\omega^{-1}) = M(\omega)x(0)$. $M(\omega)$ is a matrix with nonnegative entries, so that its principal eigenvalue $\lambda(M(\omega))$ is well-defined. We have the following (see e.g. [Benaïm et al., 2024, Thm. 2])

**Lemma 1.** If the matrices $A_1, \dots, A_N$ are Metzler and $\sum_{i=1}^N A_i$ is irreducible, then (3) holds for all non-zero initial conditions $x(0) \in \mathbb{R}_+^d$ and equals

$$ \Lambda_p(\omega) = \omega \log(\lambda(M(\omega))). $$
## 2.2 The Markovian case

For random environmental switching, let $\sigma(t)$ be a continuous-time Markov chain with transition rates $\xi_{ij} \ge 0$ from state $i$ to state $j \ne i$. Let $\Xi$ be the transition matrix with entries $\xi_{ij}$ for $i \ne j$ and $\xi_{ii} = -\sum_{j \ne i} \xi_{ij}$. Assume $\Xi$ is irreducible in which case there is a unique invariant probability measure $\alpha = (\alpha_1, \dots, \alpha_N)$ on the environmental state space $\llbracket 1, N \rrbracket$. Consider $(x(t), \sigma(\omega t))_{t \ge 0}$ where $x$ solves (4). $(x(t), \sigma(\omega t))_{t \ge 0}$ is known as a Piecewise Deterministic Markov Process (PDMP).

For a non-negative, non-zero initial condition $x(0) \in \mathbb{R}_+^d$, decompose $x(t)$ as $x(t) = \rho(t)\theta(t)$ where $\rho(t) = 1 \cdot x(t) > 0$ with $1 = (1, \dots, 1) \in \mathbb{R}^d$ and $\theta(t) = x(t)/\rho(t) \in \Delta = \{x \in \mathbb{R}_+^d, x_1 + \dots + x_d = 1\}$. In this coordinate system, the switching model (4) is

$$ \rho'(t) = (1 \cdot A_{\sigma(\omega t)}\theta(t))\rho(t), \quad \theta'(t) = F_{\sigma(\omega t)}(\theta(t)) \quad (5) $$

with

$$ F_i(\theta) = A_i \theta - (1 \cdot A_i \theta) \theta. \quad (6) $$

The following is proven in Benaim et al. [2023].

**Lemma 2.** Assuming that the matrices $A_1, \dots, A_N$ are Metzler and $\sum_{i=1}^N A_i$ is irreducible, then the Markov process $(\theta(t), \sigma(\omega t))_{t \ge 0}$ admits a unique invariant measure $\mu_\omega$ on $\Delta \times \llbracket 1, N \rrbracket$ and, moreover, the limit (3) exists almost surely, is deterministic, independent from the initial condition $x(0) \in \mathbb{R}_+^d \setminus \{0\}$ and is given by

$$ \Lambda_M(\omega) = \int_{\Delta \times \llbracket 1, N \rrbracket} 1 \cdot A_i \theta \, \mu_\omega(d\theta d\mu). \quad (7) $$

# 3 General results

Our most general results are first-order expansions of $\Lambda_p(\omega)$, $\Lambda_M(\omega)$ in terms of $\omega$ in the slow switching limit $(\omega \to 0)$ and in terms of $\frac{1}{\omega}$ in the fast switching limit $(\omega \to \infty)$. In the case of only two environments ($N=2$), we show that the first-order correction term of $\Lambda_p(\omega)$ is zero in the fast switching limit, and derive a second order approximation in terms of $\frac{1}{\omega^2}$. We also provide a simpler representation of the first-order approximation of $\Lambda_M(\omega)$ in the random case with two environments.

## 3.1 First order expansions of $\Lambda_p(\omega)$, $\Lambda_M(\omega)$ for fast and slow switching.

For the slow limit $(\omega \to 0)$, we need the right eigenvectors $x_i$ and left eigenvectors $y_i$ of $A_i$ associated with the eigenvalue $\lambda(A_i)$. We assume that these eigenvectors are normalized such that $x_i^\top y_i = 1$ and $1_d^\top x_i = 1$ where $1_d$ is the column vectors with all entries equal to 1. For the fast limit $(\omega \to \infty)$, define the averaged matrix

$$ \bar{A} = \sum_{i=1}^{N} \alpha_i A_i $$

with weights $\alpha = (\alpha_1, \dots, \alpha_N)$ given either by $\alpha_i = \tau_{i+1} - \tau_i$ in the periodic settings of Section 2.1, or by the invariant measure of $\sigma$ in the Markovian settings of Section 2.2. Let $\bar{x}$ and $\bar{y}$ the right and left eigenvectors of $\bar{A}$ associated with eigenvalue $\lambda(\bar{A})$ such that $\bar{x}^\top \bar{y} = 1$ and $1_d^\top \bar{x} = 1$.

To state our main results, recall that the matrix commutation of two $d \times d$ matrix A and B is

$$ [A, B] = AB - BA. $$

Furthermore, we need a generalized inverse for non-invertible matrices.
**Definition 1.** Let M be a matrix of index 1, meaning that M and M² have the same rank. Then, the group inverse of M is the unique solution X to the matrix equation

$$MXM = M, \quad XMX = X, \quad XM = MX.$$

When it is well-defined, we denote the group inverse of M by $M^{-1}$.

Importantly for us, if A is an irreducible, Metzler matrix, then $M = A - \lambda(A)I$ is an index 1, singular matrix. One way of computing $M^{-1}$ is to use its rank factorization. Specifically if M has rank $k < d$, then there is a $d \times k$ matrix C with rank $k$ and a $k \times d$ matrix F with rank $k$ such that $M = CF$. Then $M^{-1} = C(FC)^{-2}F$. There are multiple ways to get the rank factorization. For example, if $M = U \begin{pmatrix} \Sigma_k & 0 \\ 0 & 0 \end{pmatrix} V$ is a singular value decomposition of M where $\Sigma_k$ is a $k \times k$ diagonal matrix with the non-zero singular values of M, then one can choose C to be the first k columns of U and F to be the matrix product of $\Sigma_k$ and the first k rows of V.

The following theorem provides first order expansions of the Lyapunov exponents $\Lambda_p, \Lambda_M$ in the slow and fast limits.

**Theorem 1.** Let $A_1, \dots, A_N$ be Metzler matrices, such that $\bar{A}$ is irreducible.

1. (Slow periodic switching) If $A_1, \dots, A_N$ are irreducible, then

$$\Lambda_p(\omega) = \sum_{i=1}^{N} \alpha_i \lambda(A_i) + c_{sp}\omega + \underset{\omega \to 0}{o}(\omega) \quad (8)$$

with

$$c_{sp} = \sum_{i=1}^{N} \ln(x_i^\top y_{i+1}).$$

2. (Fast periodic switching)

$$\Lambda_p(\omega) = \lambda(\bar{A}) + c_{fp}\omega^{-1} + \underset{\omega \to \infty}{o}(\omega^{-1}) \quad (9)$$

with

$$c_{fp} = \bar{y}^\top \left( \frac{1}{2} \sum_{1 \le i < j \le N} \alpha_j \alpha_i [A_j, A_i] \right) \bar{x}.$$

For the random switching case, assume, without loss of generality that, $\max_j \sum_{i \ne j} \xi_{i,j} = 1$, and define $Q$ to be the jump matrix given by $q_{i,j} = \xi_{i,j}$ for $i \ne j$ and $q_{i,i} = 1 - \sum_{j \ne i} q_{i,j}$. Then

3. (Slow Markovian switching) If $A_1, \dots, A_N$ are irreducible, then

$$\Lambda_M(\omega) = \sum_{i=1}^{N} \alpha_i \lambda(A_i) + c_{sM}\omega + \underset{\omega \to 0}{o}(\omega) \quad (10)$$

with

$$c_{sM} = \sum_{i,j=1}^{N} \alpha_i q_{ij} \ln(y_j^\top x_i).$$
### 4. (Fast Markovian switching)

$$ \Lambda_M(\omega) = \lambda(\bar{A}) + c_{fM}\omega^{-1} + \underset{\omega \to \infty}{o}(\omega^{-1}) \quad (11) $$

with

$$ c_{fM} = \bar{y}^{\top} \left( \sum_{i,j} \alpha_i (Q - I)_{i,j}^{-1} A_j (\bar{x}\bar{y}^{\top} - I) A_i \right) \bar{x}. $$

The fast Markovian case (11) follows from [Monmarché and Strickler, 2023, Proposition 2]. The proofs for the three other cases are given in Section 6.

The signs of the first order terms determine whether the Lyapunov exponents are increasing or decreasing in the slow and fast switching limits. For instance, in the slow switching limit, a positive value of $c_{sp} > 0$ (resp. $c_{Mp}$) implies that there exists a positive frequency $\omega_0 > 0$ such that $\Lambda_p$ (resp. $\Lambda_M$) is increasing over the interval $[0, \omega_0]$. Alternatively, in the fast switching limit, a positive value of $c_{fp} > 0$ (resp. $c_{fp}$) implies that there exists a positive frequency $\omega_\infty > 0$ such that $\Lambda_p$ (resp. $\Lambda_M$) is decreasing over the interval $[\omega_\infty, \infty)$.

**Correspondence between periodic and random cases.** For the slow switching case, there is a strong correspondence between the periodic and random approximations, as we now show. Let $\alpha_i$ be given for the periodic version of the model (see section 2.1). The natural, random counterpart of this deterministic model follows the environmental states in the same order and, on average, remains in each environmental state for the same amount of time. To this end, let $\beta = \min_i \alpha_i$ and define the transition matrix $\Xi$ by $\xi_{i,i+1} = \beta/\alpha_i$ for $1 \le i \le N-1$, $\xi_{N,1} = \beta/\alpha_N$ and all other off-diagonal coefficients are 0. This transition matrix is such that the Markov process $\sigma$ visits successively the states 1, ..., N, its unique invariant probability measure is $\alpha = (\alpha_1, ..., \alpha_N)$, and $\max_i \sum_{j \ne i} \xi_{i,j} = 1$. The average time for $\sigma$ to visit the N successive states (corresponding to one period of the deterministic signal) is $1/\beta$. Therefore, $\Lambda_M(\omega)$ should be compared with $\Lambda_p(\beta\omega)$. The asymptotic expansion (10) for the slow switching rate reads

$$ \begin{aligned} \Lambda_M(\omega) &= \sum_{i=1}^{N} \alpha_i \lambda(A_i) + \omega \sum_{i,j} \alpha_i q_{i,j} \ln(y_j^{\top} x_i) + \underset{\omega \to 0}{o}(\omega) \\ &= \sum_{i=1}^{N} \alpha_i \lambda(A_i) + \omega \beta \sum_{i=1}^{N} \ln(y_j^{\top} x_i) + \underset{\omega \to 0}{o}(\omega), \end{aligned} \quad (12) $$

which is exactly the formula of $\Lambda_p(\beta\omega)$ in the periodic case (8). This correspondence doesn't occur in the fast switching regime. Indeed, even when $N = 2$, Propositions 1 and 2 in the next section provide examples where the periodic and Markovian cases differ at first order.

**The problem with reducible $A_i$ in the slow limit.** In Theorem 1, the slow switching formulae are proven under the additional assumption that each $A_i$ is irreducible. To see why this assumption is necessary, we provide an example involving three, reducible $2 \times 2$ matrices such that $\lim_{\omega \to 0} \Lambda_p(\omega) \ne \sum_i \alpha_i \lambda(A_i)$. Namely, the expansion (8) can fail at zero$^{th}$ order for reducible matrices. Consider the three, reducible matrices

$$ A_1 = \begin{pmatrix} 2 & 1 \\ 0 & 1 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix}, \quad \text{and} \quad A_3 = \begin{pmatrix} 2 & 0 \\ 1 & 1 \end{pmatrix} $$
with $\lambda(A_1) = \lambda(A_2) = \lambda(A_3) = 2$. Let $\alpha_1 = \alpha_2 = \alpha_3 = 1/3$. The matrix $\bar{A} = \frac{1}{3} \sum_i A_i$ is irreducible. Writing $T = \omega^{-1}$, one can compute

$$ \prod_{i=1}^{3} e^{A_i T/3} = e^T \begin{pmatrix} e^{2T/3} & (e^{T/3} - 1)e^{T/3} \\ (e^{T/3} - 1)e^{T/3} & (e^{T/3} - 1)^2 + e^{T/3} \end{pmatrix} $$

whose dominant eigenvalue is

$$ \frac{\sqrt{\left(2e^{\frac{5T}{3}} - e^{\frac{4T}{3}} + e^T\right)^2 - 4e^{3T}}}{2} + e^{\frac{5T}{3}} - \frac{e^{\frac{4T}{3}}}{2} + \frac{e^T}{2} = 2e^{5T/3} + O(e^{-T/3}). $$

It follows that

$$ \lim_{\omega \to 0} \Lambda_p(\omega) = \frac{5}{3} < 2 = \frac{1}{3} \sum_i \lambda(A_i). \quad (13) $$

## 3.2 Fast switching with two matrices

In this section, we investigate the case when switching only between two matrices. In the periodic case, we show that the first order term of the expansion is always 0, and give a second order expansion.

**Proposition 1.** In the fast periodic case of Theorem 1, consider the case $N = 2$. Then, $c_{fp} = 0$. Moreover, for any Metzler matrices $A_1, A_2$ of size $d \times d$,

$$ \Lambda_p(\omega) = \lambda(\bar{A}) + c_{fp,2}\omega^{-2} + \mathcal{O}\left(\frac{\omega^{-2}}{\omega}\right) \quad (14) $$

with

$$ c_{fp,2} = \bar{y}^{\top} \left( \frac{1}{12} [\alpha_2 A_2 - \alpha_1 A_1, [\alpha_2 A_2, \alpha_1 A_1]] \right) + \frac{1}{4} \alpha_1^2 \alpha_2^2 [A_2, A_1] (\bar{A} - \lambda(\bar{A})I)^{-1} [A_2, A_1] \bar{x}. $$

*Proof.* First, we show that the first order term is zero. By (9), this term is given by $c_{fp} = \frac{1}{2}\bar{y}^{\top}[\alpha_2 A_2, \alpha_1 A_1]\bar{x}$. Set $\bar{\lambda} = \lambda(\bar{A})$. By definition of $\bar{x}$ and $\bar{y}$, one has

$$ \alpha_2 A_2 \bar{x} = (\bar{\lambda} - \alpha_1 A_1) \bar{x} \quad \text{and} \quad \bar{y}^{\top} \alpha_2 A_2 = \bar{y}^{\top} (\bar{\lambda} - \alpha_1 A_1). $$

Thus,

$$ \begin{aligned} \bar{y}^{\top}[\alpha_2 A_2, \alpha_1 A_1]\bar{x} &= \bar{y}^{\top}\alpha_2 A_2 \alpha_1 A_1 \bar{x} - \bar{y}^{\top}\alpha_1 A_1 \alpha_2 A_2 \bar{x} \\ &= \bar{y}^{\top}(\bar{\lambda} - \alpha_1 A_1)\alpha_1 A_1 \bar{x} - \bar{y}^{\top}\alpha_1 A_1(\bar{\lambda} - \alpha_1 A_1)\bar{x} = 0, \end{aligned} $$

and $c_{fp} = 0$. Using once again the Baker-Campbell-Hausdorff formula, we have $M(T) = \exp(TB(T))$, with

$$ B(T) = \bar{A} + \frac{T}{2}[\alpha_2 A_2, \alpha_1 A_1] + \frac{T^2}{12}[\alpha_2 A_2 - \alpha_1 A_1, [\alpha_2 A_2, \alpha_1 A_1]] + o(T^2). $$

By Theorem 4.1 in Haviv et al. [1992], we have for a matrix $F$ and $T$ small enough (recall that $\bar{\lambda}I - \bar{A}$ admits a group inverse)

$$ \lambda(\bar{A} + TF) = \lambda(\bar{A}) + T\bar{y}^{\top}F\bar{x} + T^2\bar{y}^{\top}F(\bar{\lambda}I - \bar{A})^{-1}F\bar{x} + o(T^2). $$
Thus,

$$ \Lambda_p(T^{-1}) = \bar{\lambda} + T\bar{y}^{\top} \left( \frac{1}{2}[\alpha_2 A_2, \alpha_1 A_1] + \frac{T}{12}[\alpha_2 A_2 - \alpha_1 A_1, [\alpha_2 A_2, \alpha_1 A_1]] \right) \bar{x} + T^2 \bar{y}^{\top} \frac{1}{4} \alpha_1^2 \alpha_2^2 [A_2, A_1] (\bar{\lambda}I - \bar{A})^{-1} [A_2, A_1] \bar{x} + o(T^2), $$

which gives the announced result.

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

We now provide a simpler formula for the first order term in (11) in the Markov switching case. In this case, the transition matrix is of the form

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

where $p = \xi_{1,2} \in (0, 1)$.

**Proposition 2.** In the fast Markov case of Theorem 1, consider the case $N = 2$, with a transition matrix given by (15). Then, for any Metzler matrices $A_1, A_2$ of size $d \times d$, (11) holds with

$$ c_{fM} = p(1-p) \left[ \bar{y}^{\top}(A_1 - A_2)^2 \bar{x} - (\bar{y}^{\top}(A_1 - A_2)\bar{x})^2 \right]. \qquad (16) $$

In particular, if $A_1 - A_2$ is a diagonal matrix with entries $(\gamma_i)_{1 \le i \le d}$ on the diagonal, then

$$ c_{fM} = (1-p)p \left[ \sum_{i} \gamma_{i}^{2} \bar{x}_{i} \bar{y}_{i} - \left( \sum_{i} \gamma_{i} \bar{x}_{i} \bar{y}_{i} \right)^{2} \right] \ge 0, $$

with strict inequality as soon as the $\gamma_i$ are not all the same.

In particular, when switching randomly between two matrices, the first order term in the fast regime expansion is, in general, not zero, contrary to the periodic switching case.

*Proof.* The invariant measure of $Q$ is $\alpha = (1-p, p)$, and it is easily checked that $(Q-I)^{-1} = (Q-I)$. Hence, (11) reads

$$ c_{fM} = (1-p)p\bar{y}^{\top}(A_1PA_2 + A_2PA_1 - A_1PA_1 - A_2PA_2)\bar{x}, $$

with $P = I - \bar{x}\bar{y}^{\top}$. Straightforward computations lead to Equation (16). Finally, since $\sum_i \bar{x}_i \bar{y}_i = 1$, the fact that $c_1 \ge 0$ in the case where $A_1 - A_2$ is diagonal follows from Jensen's inequality. $\square$

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

# 4 Applications

## 4.1 Simple stage-structured models

To illustrate the use of these approximations, we examine a stage-structured model with juveniles ($x_1$) and adults ($x_2$). Juveniles mature to adults at rate 1. We consider two versions of the model: one with fluctuating birth rates and one with fluctuating mortality rates. For both versions, there are two environments and the environment transitions between its states at rate 1 when $\omega = 1$. Hence, the transition probabilities of the jump process is $Q = \begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix}$.
For the model with fluctuating birth rates, adults die at rate 1, reproduce at rate $a > 0$ in environment 1, and reproduce at rate $b > 0$ in environment 2. Assume that $a > b$ i.e. the birth rate is higher in environment 1. Hence,

$$A_1 = \begin{pmatrix} -1 & a \\ 1 & -1 \end{pmatrix} \quad \text{and} \quad A_2 = \begin{pmatrix} -1 & b \\ 1 & -1 \end{pmatrix}. \qquad (17)$$

In the slow switching limit ($\omega \to 0$), we get the approximation

$$\Lambda(\omega) = \bar{\lambda}_s + c_s \omega + o(\omega^2)$$

where

$$\bar{\lambda}_s = \frac{\sqrt{a}}{2} + \frac{\sqrt{b}}{2} - 1 \quad \text{and} \quad c_s = \log \left( \frac{\left(\sqrt{a} + \sqrt{b}\right)^2}{4\sqrt{ab}} \right).$$

The zero-th order term $\bar{\lambda}_s$ corresponds to the average $\frac{1}{2} \sum_i \lambda(A_i)$ and is positive if $\sqrt{a} + \sqrt{b} > 2$. In particular, one needs $a > 1$ to ensure that individuals have a chance of replacing themselves in their lifetime. The second order term equals zero if $a = b$ and is positive otherwise. Hence, provided there is variation in birth rates, the population growth rate increases with frequency $\omega$ (at low frequencies).

In fast randomly switching environments ($\omega \to \infty$),

$$\Lambda_M(\omega) = \lambda(\bar{A}) + c_{fM} \frac{1}{\omega} + o\left(\frac{1}{\omega^2}\right)$$

The zero-th order term equals

$$\lambda(\bar{A}) = \frac{\sqrt{2(a+b)}}{2} - 1$$

and is positive only if the average birth rate $\frac{a+b}{2}$ is greater than one. Moreover, as we have $\lambda(\bar{A}) > \bar{\lambda}_s$, the long-term population growth rate is higher in the fast limit $\omega \to \infty$ than in the slow limit $\omega \to 0$. The first order correction term is

$$c_{fM} = -\frac{(a-b)^2}{8a + 8b}.$$

This term is negative whenever there is variation in the birth rates. In the fast periodically switching environments ($\omega \to \infty$),

$$\Lambda_p(\omega) = \lambda(\bar{A}) + c_{fp,2} \frac{1}{\omega^2} + o\left(\frac{1}{\omega^3}\right)$$

where the second-order term equals

$$c_{fp,2} = -\frac{7\sqrt{2}(a-b)^2}{384\sqrt{a+b}}.$$

This term is also negative whenever there is variation in the birth rates.

Taken together these approximations suggest that (i) the long-term population growth rate $\Lambda(\omega)$ increases with frequency $\omega$ and (ii) the long-term population growth rate is higher in the periodic environment than the random environment. To explore what happens at intermediate frequencies, we numerically calculated $\Lambda_p(\omega)$ for $a = 2.5$ (good times) and $b = 0.01$ (bad times) (Fig. 1A). For these birth rates, $\Lambda(0) < 0$ and $\lim_{\omega \to \infty} \Lambda(\omega) > 0$. As predicted by the approximations, the long-term population growth rate is increasing with $\omega$ and is always lower for the random environment. Interestingly, at intermediate frequencies (around $\omega \approx 0.5$), the long-term
<figure><img src="image_5.png" /><figcaption>Figure 1: In A and B, the Lyapunov exponents for the juvenile-adult models (17) and (18), respectively, as functions of the switching frequency $\omega$. Solid lines correspond to numerically approximated Lyapunov exponents for periodic (blue) and random (red) switching. The dark grey dashed lines correspond to the analytical approximations. Parameters: In A, per-capita birth rates are $a = 2.5$ in environment 1 and $b = 0.01$ in environment 2. In B, the per-capita birth rate is $a = 100$ in both environments, the per-capita death rates for juveniles and adults are $b = 40$ and 1, respectively, in environment 1, and the per-capita death rates are 1 and $b = 40$, respectively, in environment 2.</figcaption></figure>

population growth rate is positive for a periodic environment but negative for the random environment. Hence at these intermediate frequencies, populations persist in periodically switching environments but not randomly switching environments.

For the model with fluctuating death rates, adults have a per-capita birth rate of $a$. In environment 1, adults have a per-capita death rate of $b$ and juveniles have a per-capita death rate of 1. In environment 2, juveniles have a per-capita death rate $b$, while adults have a per-capita death rate 1. Hence,

$$A_1 = \begin{pmatrix} -b & a \\ 1 & -1 \end{pmatrix} \quad \text{and} \quad A_2 = \begin{pmatrix} -1 & a \\ 1 & -b \end{pmatrix} \qquad (18)$$

We assume $b > 1$ i.e. environment 1 is worse for adults, while environment 2 is worse for juveniles. In the slow switching limit ($\omega \to 0$), we get the approximation

$$\Lambda(\omega) = \bar{\lambda}_s + c_s \omega + o(\omega^2)$$

where

$$\bar{\lambda}_s = -\frac{b}{2} + \frac{\sqrt{4a + b^2 - 2b + 1}}{2} - \frac{1}{2} \quad \text{and} \quad c_s = \log \left( \frac{4a}{4a + (b-1)^2} \right).$$

The zero-th order term $\bar{\lambda}_s$ corresponds to the average $\frac{1}{2} \sum_i \lambda(A_i)$ and is positive if the birth rate $a$ is greater than the maximal per-capita death rate $b > 1$. Under our assumption $b > 1$, the first-order term is negative. Hence, at lower frequencies, the population growth rate decreases with frequency.
Alternatively, in fast randomly switching environments ($\omega \to \infty$),

$$ \Lambda_M(\omega) = \lambda(\bar{A}) + c_{fM} \frac{1}{\omega} + o\left(\frac{1}{\omega^2}\right) $$

The zero-th order term

$$ \lambda(\bar{A}) = \sqrt{a} - \frac{b}{2} - \frac{1}{2} $$

is positive only if the square root of the average birth rate $\sqrt{a}$ is greater than the average death rate $\frac{b+1}{2}$. Moreover, as we have $\lambda(\bar{A}) > \bar{\lambda}_s$, the long-term population growth rate is higher in the fast limit $\omega \to \infty$ than the slow limit $\omega \to 0$. The first order correction term for the fast, randomly switching environment

$$ c_{fM} = \frac{(b-1)^2}{4} $$

is always positive. In fast, periodically switching environments ($\omega \to \infty$),

$$ \Lambda_p(\omega) = \lambda(\bar{A}) + c_{fp,2} \frac{1}{\omega^2} + o\left(\frac{1}{\omega^3}\right) \quad (19) $$

where the second-order term

$$ c_{fp,2} = \frac{7\sqrt{a}(b-1)^2}{96} $$

is also positive.

Taken together these approximations suggest that (i) the long-term population growth rate $\Lambda(\omega)$ decreases with frequency $\omega$ and (ii) the long-term population growth rate is lower in the periodic environment than the random environment. To explore what happens at intermediate frequencies, we numerically calculated $\Lambda_p(\omega)$ for $a = 100$ and $b = 40$ (Fig. 1B). For these birth and death rates, $\Lambda(0) > 0$ and $\lim_{\omega \to \infty} \Lambda(\omega) < 0$. As predicted by the approximations, the long-term population growth rate is decreasing with $\omega$ and is always higher for the random environment. Hence, at intermediate frequencies (around $\omega \approx 2$), the long-term population growth rate is positive for a random environment but negative for the periodic environment.

## 4.2 Metapopulation models

To model the effects of spatial heterogeneity and environmental fluctuations on population growth, we consider a metapopulation occupying $d$ distinct patches i.e. a population of $d$ populations coupled by dispersal [Hanski, 1999]. The population density in patch $i$ is $x_i$ with a per-capita growth rate of $r_{\sigma(\omega t)}^i$ at time $t$. Individuals from patch $j$ disperse to patch $i$ at rate $l_{i,j}$. The population dynamics in patch $i$ are

$$ x'_{i}(t) = r_{\sigma(\omega t)}^{i} + \sum_{j \neq i} l_{i,j}(x_{j} - x_{i}). \quad (20) $$

In matrix form, the metapopulation dynamics are

$$ x'(t) = \underbrace{(R_{\sigma(\omega t)} + L)}_{=:A_{\sigma(\omega t)}} x(t). \quad (21) $$

where $R_\sigma$ is a diagonal matrix with entries $r_\sigma^1, \dots, r_\sigma^d$ and $L$ is the dispersal matrix with entries $l_{i,j}$.

For a symmetric dispersal matrix and continuous (rather than piece-wise continuous), periodic $r^i: \mathbb{R} \to \mathbb{R}$, Katriel [2022] proved the following result about monotonicity of the $\Lambda_p(\omega)$.
<figure><img src="image_6.png" /><figcaption>Figure 2: The Lyapunov exponents and density-dependent dynamics for a two patch metapopulation model (21). In A, the Lyapunov exponent as a function of the switching frequency $\omega$. Solid lines correspond to numerically approximated Lyapunov exponents for periodic (blue) and random (red) switching. The dark gray, dashed lines correspond to the analytical approximations. Random (in C) and periodic (in C) simulations at frequency</figcaption></figure>

$\omega = 0.1$ for a nonlinear model whose linearization at the origin is given by (21). Parameters: $L = \begin{pmatrix} -0.1 & 0.1 \\ 0.1 & -0.1 \end{pmatrix}$, $R_1 = \begin{pmatrix} 0.5 & 0 \\ 0 & -1.5 \end{pmatrix}$, and $R_2 = \begin{pmatrix} -1.5 & 0 \\ 0 & 0.5 \end{pmatrix}$. For B and C, the dynamics are $x'(t) = A_{\sigma(\omega t)}x(t) - 0.01x(t) \circ x(t)$ where $\circ$ denotes the Hadamard product.
**Theorem 2.** [Katriel, 2022] Assume L is irreducible and symmetric, $r^i: \mathbb{R} \to \mathbb{R}$ are continuous, period 1 functions that are not all equal. Then the dominant Floquet multiplier $\Lambda_p(\omega)$ of

$$x'_{i}(t) = r^{i}(\omega t) + \sum_{j \neq i} l_{i,j}(x_{j} - x_{i}).$$

is a strictly decreasing function of $\omega$.

Theorem 2 suggests that slower switching promotes metapopulation persistence i.e. a higher growth rate $\Lambda_p(\omega)$. In particular, if $\sum_{i=1}^N \alpha_i \lambda(A_i) > 0 > \lambda(\sum_{i=1}^N \alpha_i A_i)$, then there is a critical frequency $\omega^*$ such that the metapopulation persists if and only if the frequency $\omega$ of environmental switching lies below $\omega^*$ (see black curve in Figure 2).

Using Theorem 1, we show that Theorem 2 may extend to random and periodic, piecewise continuous switching and may partially extend to asymmetric dispersal matrices L.

**Theorem 3.** Assume L is irreducible and the $R_i$ in (21) are all not equal.

1. (Fast switching) $c_{fM} > c_{fp} = 0$.

2. (Slow switching) If L is symmetric or $d = 2$, then $c_{sp} \le 0$ and $c_{sM} \le 0$. Moreover, equality holds if and only if the right eigenvectors $x_i$ of $A_i = R_i + L$ are all equal.

In the case of random switching, Theorem 3 implies that sufficiently fast switching always leads to lower population growth rates. Specifically there exists $\omega^* \ge 0$ such that $\omega \mapsto \Lambda_M(\omega)$ decreases with the frequency $\omega$ of environmental switching for $\omega \ge \omega^*$. This occurs for all forms of dispersal - symmetric or asymmetric. If dispersal is symmetric or there are only two patches, Theorem 3 ensures that $\Lambda(\omega)$ decreases for sufficiently slow frequencies. Hence, we make the following conjecture for the random case:

**Conjecture 1.** Assume L is irreducible and the eigenvectors $x_i$ of $A_i = R_i + L$ in (21) are all not equal. If L is symmetric or $d = 2$, then $\omega \mapsto \Lambda_M(\omega)$ is strictly decreasing.

Figure 2 illustrates Theorems 2 and 3 and Conjecture 1 numerically for a two-patch, two environment model ($d = N = 2$). Both for random and periodic switching, the population growth rate $\Lambda(\omega)$ decreases with the switching frequency. Hence, the critical frequency $\omega^*$ below which the metapopulation persists is higher for the randomly switching environment than the periodically switching environment.

What happens when there are $d \ge 3$ patches and the dispersal matrix L is asymmetric? It turns out that $c_{sM}$ and $c_{sp}$ can be positive and, consequently, $\omega \to \Lambda(\omega)$ need not be monotonic. To illustrate this possibility, consider patches lying along a circle and individuals moving clockwise along this circle. Specifically, let L be the circular permutation matrix:

$$L = \begin{pmatrix} -1 & 0 & \dots & 0 & 1 \\ 1 & -1 & \dots & 0 & 0 \\ 0 & 1 & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & -1 \end{pmatrix}$$

Namely $L_{i,i} = 1$, $L_{1,d} = 1$, $L_{i+1,i} = 1$ for $1 \le i \le d-1$, and $L_{i,j} = 0$ otherwise. We assume there are $N = d$ environmental states and the per-capita growth of patch i in environmental state i equals $\rho+\beta \in \mathbb{R}$ and otherwise is $\rho \in \mathbb{R}$. Namely, $r_i^i = \beta + \rho$ and $r_i^j = \rho$ for $i \ne j$. Finally, we consider the Markovian signal $\sigma$ with transition
<figure><img src="image_7.png" /><figcaption>Figure 3: The plain arrows represent the population migration, the dotted arrows the random transitions of the signal, i.e. of the state where there is a non-zero growth rate. Left: $\Xi = L^T$, the signal moves in the same direction as the population. Right: $\Xi = L$, the signal moves in the reverse direction.</figcaption></figure>

matrix $\Xi$ being either $L$ or $L^T$. When $\Xi = L^T$, the signal moves in the same direction as the population (we call this the synchronized case), when $\Xi = L$ it goes backward (we call this the asynchronous case), see Figure 3. As discussed at the end of Section 3.1, the first order term expansion of $\Lambda_p(\omega)$ when $\sigma$ deterministically goes from $i$ to $i+1$ (synchronized case) or from $i$ to $i-1$ (asynchronous) is given by $c_{sp} = dc_{sM}$.

**Proposition 3.** Let $\eta$ be the unique positive solution of $(\eta - \beta)\eta^{d-1} = 1$.

1. (synchronized case) If $\Xi = L^T$, then

$$c_{sM} = c_{sM}^{sync} := \ln \left( \frac{(d-2)(\eta-\beta) + 2\eta}{(d-1)(\eta-\beta)\eta + \eta^2} \right).$$

In the case $d=3$, the sign of $c_{sM}$ is the opposite of the sign of $\beta$. For $d > 3$, its sign is the sign of $\beta_d - \beta$ where $\beta_d = \eta_* - \eta_*^{1-d} > 0$ with $\eta_* > 2(d-1)/(d+1) > 1$ the unique root of $X^{d+1} - 2X^d + (d-1)X - d + 2$ over $(1, \infty)$.

2. (asynchronous case) If $\Xi = L$, then

$$c_{sM} = c_{sM}^{async} := \ln \left( \frac{d\eta}{d-1 + \eta^d} \right),$$

which is negative for all $\beta \neq 0$, $d \ge 3$.

For the periodic case where $\sigma$ deterministically goes from $i$ to $i+1$ (synchronized case) or from $i$ to $i-1$ (asynchronous case), the first order term expansion of $\Lambda_p(\omega)$ is given by $c_{sp} = dc_{sM}$.

The Proof of Proposition 3 is in Section 8. Figure 4 illustrates the main conclusions of this proposition with $d=3$ patches, $\beta = -1$, and $\rho = 0.3$.
<figure><img src="image_8.png" /><figcaption>Figure 4: The Lyapunov exponent $\Lambda_p(\omega)$ from the circular metapopulation model with $d = 3$ patches. In panel A, the asynchronous case where $\Xi = L$. In panel B, the synchronous case where $\Xi = L^r$. The per-capita growth $r_i^j$ in patch $i$ equals $-0.7$ in environment $j = i$ and $0.3$ otherwise.</figcaption></figure>

# 5 Discussion

Populations experience environmental fluctuations at multiple time scales from diurnal cycles to multi-decadal climate cycles [Hasegawa et al., 2022, Gorenstein et al., 2023]. Here, we derived analytical approximations for how the frequency of these environmental fluctuations influence the long-term growth rate of structured populations. In the limit of low-frequency fluctuations, we derived new analytic approximations of the form $\Lambda(\omega) = \sum_i \alpha_i \lambda(A_i) + c\omega + o(\omega)$ with explicit formulas for $c$ in the periodic and random cases. In this low frequency limit, we showed that the mode of fluctuations, random versus periodic, has no effect on the population growth rate i.e. the correction terms $c_p$ and $c_M$ are equivalent when the random and periodic signal are comparable. This differs sharply from the high-frequency limit. In the limit of high-frequency fluctuations, we derived an analytical approximation for periodic environments complementing the work of Monmarché and Strickler [2023] for random environments. For both approximations, the population growth rate in the infinitely fast switching limit are equal (i.e. $\lim_{\omega \to \infty} \Lambda_p(\omega) = \lim_{\omega \to \infty} \Lambda_M(\omega) = \lambda(\sum_i \alpha_i A_i)$), but the first-order correction terms with respect to the period $1/\omega$ of fluctuations, in general, differ (i.e. $c_{fp} \neq c_{fM}$). Indeed, for switching between two environments, the first-order correction term is zero for the periodic case, but (in general) non-zero for the random case. This differences implies that whether populations persist or not may depend on the mode of environmental switching.

## Biological interpretations of the slow and fast limits

The slow limit approximation has a clear biological interpretation in terms of the stable state distributions $x_i$ and the vectors of reproductive values $y_i$ – two fundamental quantities of population modeling [Caswell, 2001]. When environmental shifts are rare, the distribution of the population approaches the stable state distribution of the current environment, say $x_i$ in environment $i$. Reproductive values for this environment, the components of $y_i$, correspond to the contributions of an individual in each state to the long-term population [Caswell, 2001]. This interpretation of the absolute, rather than relative, contributions to the growth rate relies on the vector
of reproductive values being normalized, so that, the expected value of the reproductive value of a randomly
chosen individual from the stable state distribution is one i.e. $y_i^T x_i = 1$ [Caswell, 2001]. When the environment
shifts to a new environment, say $j$, the expected reproductive value of a randomly chosen individual shifts to
$y_j^T x_i$ i.e. the expected reproductive value in the new environment of a randomly chosen individual from the old
environment. When this new expectation is greater than one, individuals, on average, have greater reproductive
value immediately after the environmental shift. This increase in reproductive value boosts the population
growth rate by $\log(y_j^T x_i)$. If environmental shifts, on average, boost reproductive values, then the long-term
population growth rate $\Lambda(\omega)$ increases with frequency $\omega$– at least for $\omega$ sufficiently low. This phenomena is
illustrated with our model (17) of a population of juveniles and adults where birth rates fluctuate between low
and high values. In low birth rate environments, the stable state distribution mostly consists of adults and
adults have slightly higher reproductive values. In the high birth rate environments there are mostly juveniles
in the stable state distribution but adults have much high reproductive value. Hence, rare shifts from the low
birth rate to high birth rate environments substantially increase the average reproductive value of an individual.
In contrast, rare shifts from the high birth rate to low birth rate environments slightly decrease the average
reproductive value of an individual. Hence, the net effect of environmental shifts is positive and the long-term
population growth rate $\Lambda(\omega)$ increases with frequency. In contrast, when the juvenile and adult death rates
are fluctuating i.e. equation (18), we get the opposite trend. In the environment with higher adult mortality,
there are more juveniles at the stable state distribution and juveniles have higher reproductive value. In the
environment with higher juvenile adult mortality, adults are more common at the stable state distribution and
adults have higher reproductive value. Hence, switching between environments always results in a decrease in
the average reproductive value of an individual and a corresponding decrease in the population growth rate.
Similar reasoning applies to the two patch metapopulation model in Figure 2: switching always results in more
individuals in the lower quality patch and a reduction in the population growth rate.

The interpretation of the fast limit approximations are less clear in general. However, we can draw some
conclusions in the case of random switching between two environments. When only the contribution of one
population state to another population state fluctuates (i.e. a non-diagonal entry of $A$), Proposition 2 implies
that higher frequency fluctuations always increase the long-term population growth rate. In contrast, when
only the contribution of population states to themselves rapidly fluctuate (i.e. only diagonal terms fluctuate),
higher frequencies decrease the long-term population growth rate. In particular, this conclusion applies to
metapopulations with constant dispersal rates whose per-capita growth rates fluctuate between two values.

Implications for metapopulations

For metapopulations with symmetric dispersal matrices, Katriel [2022] showed that the long-term population
growth rate $\Lambda_p(\omega)$ always decreases with frequency when per-capita growth rates exhibit periodic, continuous
fluctuations. Hence, metapopulations exhibiting diffusive-like movement grow faster in low frequency, periodic
environments. Here, we showed the same conclusion holds for randomly fluctuating environments in the limits
of fast and slow switching (see Theorem 3). Moreover, in the fast switching limit, long-term population growth
rates are higher in randomly switching environments than periodically switching environments. An open problem
(see Conjecture 1) is whether these conclusions hold at intermediate frequencies of random switching. When
metapopulations exhibit asymmetric dispersal, however, we show that the long-term population growth rate
may not vary monotonically with the environmental frequency $\omega$. For example, we constructed an example (see
Figure 4) in which $\Lambda(\omega)$ is maximized at intermediate environmental frequencies. A likely hard open problem is
characterizing which asymmetric dispersal matrices result in similar behaviors.
<figure><img src="image_9.png" /><figcaption>Figure 5: Demographic variances (22) and (23) for the stage structured models (17)-(18) with fluctuating birth rates (A) and death rates (B). Parameters are as in Figure 1</figcaption></figure>

## Relationship to stochastic demography

The effect of temporal variation in environmental conditions on the growth rate of structured population is a primary focus of the field of stochastic demography [Tuljapurkar, 1990, Caswell, 2001, Boyce et al., 2006, Hastings and Gross, 2012]. This extensive body of work, however, focuses exclusively on discrete-time models. Hence, to make comparisons to this literature requires defining an appropriate metric of variation of instantaneous demographic rates, say an entry $a(\omega t)$ of the matrix $A_{\sigma(\omega t)}$. One natural option is to consider the variation in the average demographic rate at stationarity over a fixed time interval, say of length one. In a random environment with stationary distribution $\pi$, this metric of variation equals

$$ 
\begin{aligned}
\mathrm{Var}_M(\omega) &= \mathbb{E}_\pi \left[ \left( \int_0^1 a(\omega t) dt \right)^2 \right] - \mathbb{E}_\pi \left[ \int_0^1 a(\omega t) dt \right]^2 \\
&= \mathbb{E}_\pi \left[ \left( \frac{1}{\omega} \int_0^\omega a(t) dt \right)^2 \right] - \mathbb{E}_\pi \left[ \frac{1}{\omega} \int_0^\omega a(t) dt \right]^2
\end{aligned}
\quad (22) $$

where $\mathbb{E}_\pi$ corresponds to the expectation with $P[\sigma(0) = i] = \pi_i$. $\mathrm{Var}_M(\omega)$ is a smooth function of $\omega$ and converges to zero as $\omega \to 0$ (see, e.g. [Monmarché and Strickler, 2023]). Similarly, one can define variation of $a$ in the periodic environment as

$$ \mathrm{Var}_p(\omega) = \int_0^1 \left[ \left( \frac{1}{\omega} \int_0^\omega a(s+t)dt \right)^2 \right] ds - \left[ \int_0^1 \frac{1}{\omega} \int_0^\omega a(t+s)dt ds \right]^2 \quad (23) $$

where $a(t)$ is started randomly in the interval $[0, 1]$. $\mathrm{Var}_p(\omega)$ also converges to zero as $\omega \to \infty$. However, unlike the random case, this convergence isn't monotone as $\mathrm{Var}_p(\omega)$ equals zero at integer values. None the less, as
illustrated in Figure 5, average birth and death rates for the juvenile-adult models tend to decrease with frequency for the periodic case and continually decrease for random fluctuations. These trends facilitates comparisons to the stochastic demography literature.

A central dogma of stochastic demography is [Hastings and Gross, 2012, page 96] “temporal variation in the vital rates will typically lead to a reduction in population growth rate.” This genesis of this dogma were discrete-time geometric random walks and the observation that geometric mean is less than the arithmetic mean [Lewontin and Cohen, 1969]. For random, discrete-time matrix models $x(t + 1) = A(t + 1)x(t)$, Tuljapurkar [1990] extended this dogma provided the non-negative matrices $A(1), A(2),...$ are serially uncorrelated: the dominant Lyapunov exponent of this random product decreases with variation in the matrix entries. The continuous-time models considered here, however, are inherently temporally autocorrelated over sufficiently short time scales. Consistent with our results, the central dogma of stochastic demography doesn’t apply, in general, to discrete-time models with temporal auto-correlations: temporally auto-correlated fluctuations can increase population growth [Tuljapurkar, 1990]. This “inflationary effect” of auto-correlated fluctuations have been observed in stage-structured [Tuljapurkar and Haridas, 2006] and spatially structured models [Roy et al., 2005, Tuljapurkar and Haridas, 2006, Schreiber, 2010, Kortessis et al., 2023a, Schreiber, 2023]. Consistent with these earlier studies, we found fluctuations in mortality for stage-structured models and per-capita growth rates in spatially-structured models can generate this inflationary effect: lower frequencies of fluctuations generate greater variation in demographic rates averaged over a time step and higher population growth rates. Notably, our approximations suggest that this inflationary effect in spatially-structured population models occurs whenever there is symmetric dispersal.

## Concluding remarks

Our approach of modeling environmental fluctuations as a piece-wise constant, continuous process offers an analytically more tractable method for assessing the impact of environmental variability on the long-term growth rates of structured populations. The analytical approximations we’ve developed for growth rates in the limits of slow and fast environmental switching may serve as valuable new tools for theoretical biologists and applied mathematicians. These tools enable a deeper examination of how both the tempo (slow versus fast) and mode (periodic versus random) of environmental fluctuations influence critical population dynamics, including growth, persistence, and extinction risk. This refined approach not only enhances our understanding of population responses to environmental variability but also opens new avenues for research in population biology and applied mathematics with potential applications to epidemiology, ecology, population genetics, and conservation biology.

# 6 Proof of Theorem 1

## 6.1 Slow periodic switching

We prove formula (8). Recall that $T = \omega^{-1}$. We write $r_i = \mu(A_i)$. By Perron-Frobenius Theorem, for all $i$, for all $T > 0$,

$$e^{-r_i \alpha_i T} e^{\alpha_i T A_i} = p_i + o(e^{-\gamma T}),$$

where $\gamma > 0$ and $p_i = x_i y_i^T$ is the projection matrix on $x_i$ in the direction of $y_i$. Thus,

$$\prod_{i=1}^{N} e^{-r_i \alpha_i T} e^{T \alpha_i A_i} = \prod_{i=1}^{N} p_i + o(e^{-\gamma T}),$$
which leads to

$$ \lambda(M(T)) = \prod_{i=1}^{N} e^{r_i \alpha_i T} \lambda \left( \prod_{i=1}^{N} p_i + o(e^{-\gamma T}) \right) $$

Now, $\lambda(\prod_{i=1}^{N} p_i + o(e^{-\gamma T})) = \lambda(\prod_{i=1}^{N} p_i) + o(\frac{1}{T})$ and

$$ \prod_{i=1}^{N} p_i = x_N y_N^T \dots x_1 y_1^T = \left( \prod_{i=1}^{N-1} x_i^T y_{i+1} \right) x_N y_1^T. $$

In particular, $\lambda(\prod_{i=1}^{N} p_i) = \prod_{i=1}^{N-1} x_i^T y_{i+1} \lambda(x_N y_1^T)$, where $x_N y_1^T$ is the projection matrix on $x_N$ in the direction of $y_1^T$. Hence, $\lambda(x_N y_1^T) = x_N^T y_1$. This concludes the proof by Lemma 1.

## 6.2 Fast periodic switching

This section is devoted to the proof of Formula (9). On the one hand, for an irreducible Metzler matrix $B$, we have

$$ \lambda(e^B) = e^{\lambda(B)}. $$

Indeed, since $B$ is Metzler and irreducible, $e^B$ is a nonnegative irreducible matrix, hence $\lambda(e^B) = r(e^B)$. Now it is easily seen that $r(e^B) = e^{\lambda(B)}$ because the spectrum of $e^B$ is the exponential of the spectrum of $B$. On the other hand, by an iteration of the Baker - Campbell - Hausdorff formula, one has

$$ M(T) = \exp \left( T \bar{A} + \frac{T^2}{2} \sum_{1 \le i < j \le N} [\alpha_j A_j, \alpha_i A_i] + o(T^2) \right) $$

Thus,

$$ \frac{1}{T} \log \lambda(M(T)) = \lambda(B(T)), \qquad (24) $$

where $B(T) = \bar{A} + \frac{T}{2} \sum_{1 \le i < j \le N} \alpha_j \alpha_i [A_j, A_i] + o(T)$. Now, since $\bar{A}$ has a unique dominant eigenvalue, the same holds for $B(T)$, for $T$ small enough. Moreover, one has $B'(0) = \frac{1}{2} \sum_{1 \le i < j \le N} \alpha_j \alpha_i [A_j, A_i]$. Thus, by [Horn and Johnson, 2013, Theorem 6.3.12] or [Haviv et al., 1992, Theorem 4.1], the function $T \mapsto \lambda(B(T))$ is differentiable at 0 and its derivative at 0 is given by $\bar{y}^T B'(0) \bar{x} = \bar{y}^T (\frac{1}{2} \sum_{1 \le i < j \le N} \alpha_j \alpha_i [A_j, A_i]) \bar{x}$. This, together with Lemma 1 and Equation (24), proves (9).

## 6.3 Slow Markovian switching

In this section, devoted to the proof of (10), $\sigma$ is a Markov chain as in Section 2.2. We use the notations of Section 2.2. In view of (7), the goal is now to give an expansion of $\mu_\omega$ when $\omega$ vanishes.

In the following, for $i \in [1, N]$, we denote by $(\varphi_t^i)_{t \ge 0}$ the flow on $\Delta$ associated to $F_i$; namely, for $z \in \Delta$, $\theta_i(t) := \varphi_t^i(z)$ solves $\theta_i'(t) = F_i(\theta_i(t))$ for all $t \ge 0$ with $\theta_i(0) = z$. Recall that $x_i$ stands for a unit right eigenvector of $A_i$. The following follows e.g. from [Benaïm and Strickler, 2019, Lemma 6].

**Lemma 3.** Assuming that $A_1, \dots, A_N$ are Metzler and irreducible, there exist $C, a > 0$ such that for all $z \in \Delta$, $i \in [1, N]$ and $t \ge 0$,

$$ |\varphi_t^i(x) - x_i| \le Ce^{-at}. \qquad (25) $$
In the slow switching regime, since the signal $\sigma(\omega t)$ stays constant for very long times, Lemma 3 shows that the process $\theta$ in (5) spends most of its time close to the points $(x_i)_{i \in [1, N]}$, and $\mu_\omega$ is then expected to be close to

$$ \tilde{\mu}_0 = \sum_{i=1}^{N} \alpha_i \delta_{x_i, i}. $$

The next statement provides the first order expansion of this convergence. Notice that its proof doesn't rely on the specific form of the vector fields (6), but only on the contraction established in Lemma 3 (so that the next result can be straightforwardly extended to slow switching between any contracting flows – i.e. assuming (25) – over a compact manifold)

**Proposition 4.** *Assuming that $A_1, \dots, A_N$ are Metzler and irreducible, for all $f \in C^1(\Delta \times [1, N])$,*

$$ \mu_\omega f = \tilde{\mu}_0 f + \omega \tilde{c}_1 + \mathcal{O}_{\omega \to 0}(\omega) $$

with

$$ \tilde{c}_1 = \sum_{i,j=1}^{N} \alpha_i q_{ij} \int_{0}^{\infty} (f(\varphi_s^j(x_i), j) - f(x_j, j)) ds. $$

*Proof.* Fix $f \in C^1(\Delta \times [1, N])$. Denote by $(T_k)_{k \in \mathbb{N}}$ the jump times of the Markov chain $(\sigma(\omega t))_{t \ge 0}$ and by $\nu_\omega$ the invariant measure of the discrete-time chain $(\theta(T_k), \sigma(\omega T_k))_{k \in \mathbb{N}}$, called the skeleton chain associated to $(\theta(t), \sigma(\omega t))_{t \ge 0}$. Recall the following relation between $\nu_\omega$ and $\mu_\omega$ (see [Davis, 1993, Theorem 34.31 p.123]):

$$ \mu_\omega f = \int_{\Delta \times [1, N]} \int_0^\infty f(\varphi_s^i(z), i) \omega e^{-\omega s} ds \nu_\omega(dz di). \quad (26) $$

Thanks to Lemma 3, for any $(z, i) \in \Delta \times [1, N]$ we can write

$$ \begin{aligned} \int_0^\infty f(\varphi_s^i(z), i) \omega e^{-\omega s} ds &= f(x_i, i) + \omega \int_0^\infty (f(\varphi_s^i(z), i) - f(x_i, i)) e^{-\omega s} ds \\ &= f(x_i, i) + \omega \int_0^\infty (f(\varphi_s^i(z), i) - f(x_i, i)) ds + \mathcal{O}_{\omega \to 0}(\omega), \end{aligned} \quad (27) $$

by dominated convergence. On the other hand, since $\nu_\omega$ is invariant for the Markov kernel $\mathcal{T}$ given by

$$ \mathcal{T} f(z, i) = \sum_{j=1}^{N} q_{ij} \int_{0}^{\infty} f(\varphi_{s}^{i}(z), j) \omega e^{-\omega s} ds, $$

using that $\nu_\omega = \nu_\omega \mathcal{T}$, that the marginal on $[1, N]$ of $\nu_\omega$ is $\alpha$ and that, thanks to Lemma 3,

$$ \left| \mathcal{T} f(x, i) - \sum_{j=1}^{N} q_{ij} f(x_i, j) \right| \le C \| \nabla_x f \|_{\infty} \omega $$

for some $C > 0$ independent from $\omega$, we immediately get

$$ \nu_\omega f \underset{\omega \to 0}{\longrightarrow} \sum_{i,j=1}^{N} \alpha_i q_{ij} f(x_i, j). $$
Using this after (27) (and, again, that the marginal of $\nu_\omega$ on $[1, N]$ is $\alpha$) in (26) we end up with

$$ 
\begin{aligned}
\mu_\omega f &= \tilde{\mu}_0 f + \omega \int_{\Delta \times [1,N]} \int_0^\infty (f(\varphi_s^j(z), i) - f(x_j, j)) \, ds \, \nu_\omega(dz) \, di + \mathcal{o}_{\omega \to 0}(\omega) \\
&= \tilde{\mu}_0 f + \omega \sum_{i,j=1}^N \alpha_i q_{ij} \int_0^\infty (f(\varphi_s^j(x_i), j) - f(x_j, j)) \, ds + \mathcal{o}_{\omega \to 0}(\omega).
\end{aligned}
$$

This concludes the proof.

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

In view of (7) and Proposition 4, the proof of the expansion (10) in the slow Markovian switching case is concluded by the following.

**Lemma 4.** In Proposition 4, for $f(z, i) = 1 \cdot A_i z$, we have

$$ \tilde{c}_1 = \sum_{i,j=1}^{N} \alpha_i q_{ij} \ln(y_j^T x_i). $$

*Proof.* Notice that, for all $z \in \Delta_d$,

$$ f(\varphi_s^j(z), j) = \frac{\mathbf{1}_d^T A_j e^{sA_j} z}{\mathbf{1}_d^T e^{sA_j} z} = \frac{d}{ds} \ln(\mathbf{1}_d^T e^{sA_j} z). $$

Theorefore, using that $x_j = \varphi_s^j(x_j)$, one has

$$ 
\begin{aligned}
\int_0^\infty (f(\varphi_s^j(x_i), j) - f(x_j, j)) \, ds &= \int_0^\infty \frac{d}{ds} \left( \ln \left( \frac{\mathbf{1}_d^T e^{sA_j} x_i}{\mathbf{1}_d^T e^{sA_j} x_j} \right) \right) \, ds \\
&= \lim_{s \to \infty} \ln \left( \frac{\mathbf{1}_d^T e^{sA_j} x_i}{\mathbf{1}_d^T e^{sA_j} x_j} \right) \\
&= \ln(y_j^T x_i),
\end{aligned}
$$

where we have used (once again) Perron-Frobenius Theorem.

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

# 7 Proof of Theorem 3

We first prove a more general result about random switching of Metzler matrices that are symmetric. This result applies to the special case of (21) when the dispersal matrix $L$ is symmetric.

**Proposition 5.** In the settings of Theorem 1, assume that all the matrices $A_1, \dots, A_N$ are symmetric. Then, in (8) and (10), $c_{sp} \le 0$ and $c_{sM} \le 0$. Moreover, either all the $x_i$ are equal, in which case $\omega \mapsto \Lambda_p(\omega)$ and $\omega \mapsto \Lambda_M(\omega)$ are constant, or otherwise $c_{sp} < 0$ and $c_{sM} < 0$.

*Proof.* We only give the proof for the Markov case, since in the slow switching regime, $c_{sp}$ is positively proportional to $c_{sM}$ for a particular Markov chain, as described at the end of Section 3.1. Since $A_i$ is symmetric, its left
eigenvector $y_i$ is proportional the right eigenvector $x_i$. Moreover, due to the normalization condition $y_i^\top x_i = 1$, we have $y_i = \|x_i\|^{-2}x_i$. Therefore,

$$ \prod_{i,j} (y_j^\top x_i)^{q_{i,j} \alpha_i} = \prod_{i,j} (\|x_j\|^{-2} x_j^\top x_i)^{q_{i,j} \alpha_i} \\ \le \prod_{i,j} \left( \frac{\|x_i\|}{\|x_j\|} \right)^{q_{i,j} \alpha_i}, $$

where the last line follows from Cauchy-Schwarz inequality. Now, in this last product, the exponent of each $\|x_i\|$ is $\sum_j \alpha_i q_{i,j} - \sum_j \alpha_j q_{j,i}$, which is zero since $\alpha$ is invariant for $Q$. This entails that $c_{sM} \le 0$. Moreover, $c_{sM} = 0$ if and only if there is equality in all the Cauchy-Schwarz inequalities we used, namely, if and only if all the $x_i$ are positively correlated, and thus, equal. In that case, one easily checks that $\mu_\omega$, the invariant probability measure of the PDMP on the sphere considered in Section 6.3 for the proof of (10), is $\delta_{x_1} \otimes \alpha$ for all values of $\omega$. In view of (7), $\omega \mapsto \Lambda_M(\omega)$ is thus constant. In the periodic case, $\Lambda_p(\omega)$ can also be written as the integral of a function of the unique globally asymptotically periodic solution of the system on the sphere (see Formula 29 in Theorem 8 in Benaim et al. [2024]). Here again, it is easy to check that the constant solution equal to $x_1$ is the unique periodic solution, and thus $\Lambda_p(\omega)$ is constant. $\square$

The next two propositions provide show that $c_{fM} > c_{fp} = 0$ in the fast-switching limit.

**Proposition 6.** In the case (21), $c_{fp} = 0$.

*Proof.* Recall that

$$ c_{fp} = \frac{1}{2} \sum_{i<j} \alpha_i \alpha_j \bar{y}^\top [A_j, A_i] \bar{x}. $$

Since $A_i = R_i + L$ and $R_i$ and $R_j$ commute, we end up with $[A_j, A_i] = [L, R_i] + [R_j, L]$. Now,

$$ \begin{aligned} \bar{y}^\top [L, R_i] \bar{x} &= \bar{y}^\top L R_i \bar{x} - \bar{y}^\top R_i L \bar{x} \\ &= \bar{y}^\top (\lambda(\bar{A}) - \bar{R}) R_i \bar{x} - \bar{y}^\top R_i (\lambda(\bar{A}) - \bar{R}) \bar{x} \\ &= 0, \end{aligned} $$

where we have used that $\bar{R}$ and $R_i$ commute and that by definition,

$$ \bar{y}^\top (\bar{R} + L) = \lambda(\bar{A}) \bar{y}^\top, \quad (\bar{R} + L) \bar{x} = \lambda(\bar{A}) \bar{x}. $$

**Proposition 7.** In the case (21), $c_{fM} > 0$.

*Proof.* First, let us prove that, in the case (21),

$$ c_{fM} = \sum_{i,j} \alpha_i (Q - I)_{i,j}^{-1} (\bar{y}^\top R_j \bar{x} \bar{y}^\top R_i \bar{x} - \bar{y}^\top R_j R_i \bar{x}) . \quad (28) $$

Indeed,

$$ \bar{y}^\top A_j \bar{x} \bar{y}^\top A_i \bar{x} - \bar{y}^\top A_j A_i \bar{x} = \bar{y}^\top R_j \bar{x} \bar{y}^\top R_i \bar{x} - \bar{y}^\top R_j R_i \bar{x} + B_i + C_j + D, $$
where

$$B_i = \bar{y}^\top L \bar{x} \bar{y}^\top R_i \bar{x} - \bar{y}^\top L R_i \bar{x}$$

$$C_j = \bar{y}^\top R_j \bar{x} \bar{y}^\top L \bar{x} - \bar{y}^\top R_j L \bar{x}$$

$$D = \bar{y}^\top L \bar{x} \bar{y}^\top L \bar{x} - \bar{y}^\top L^2 \bar{x}.$$

Using that $\sum_i \alpha_i (Q - I)_{ij}^{-1} = 0 = \sum_j (Q - I)_{i,j}^{-1}$, we end up with $\sum_{i,j} \alpha_i (Q - I)_{i,j}^{-1} (B_i + C_j + L) = 0$, leading to (28).

For $k \in [1, N]$, write $a_k = \bar{x}_k \bar{y}_k$ and $(R_i)_{k,k} = r^{i,k}$. Since the $R_i$'s are diagonal, we get that $\bar{y}^\top R_i \bar{x} = \sum_k a_k r^{i,k}$. Now, let $X = (X_1, \dots, X_N)$ be a random variable taking the value $(r^{1,k}, \dots, r^{N,k})$ with probability $a_k$ for all $k \in [1, N]$ and define $\varphi: \mathbb{R}^N \to \mathbb{R}$ by

$$\varphi : (x_1, \dots, x_N) \mapsto -\sum_{i,j} \alpha_i (Q - I)_{i,j}^{-1} x_i x_j.$$

It is readily checked that

$$\sum_{i,j} \alpha_i (Q - I)_{i,j}^{-1} (\bar{y}^\top R_j \bar{x} \bar{y}^\top R_i \bar{x} - \bar{y}^\top R_j R_i \bar{x}) = \mathbb{E}(\varphi(X)) - \varphi(\mathbb{E}(X)).$$

It remains to prove that $\varphi$ is convex, which will imply by Jensen inequality that $c_{fM} > 0$. Since $\varphi$ is a quadratic form, it is convex if and only if it is positive. Yet, $\varphi(x) = -\int f(i)(Q - I)^{-1} f(i) \alpha(di)$, where $f(i) = x_i$, which can be interpreted as an asymptotic variance, and is thus positive (see the discussion after Proposition 4 in Monmarché and Strickler [2023]). $\square$

Finally, we consider the case of $d = 2$ patches. The matrices $A_i$ in case (21) can be written as

$$A_i = \begin{pmatrix} a_i & \alpha \\ \beta & b_i \end{pmatrix}, \quad i = 1, \dots, N.$$

for some $a_i, b_i \in \mathbb{R}$ and $\alpha, \beta > 0$ (recall that the matrices are Metzler and irreducible). The right and left eigenvectors can be computed explicitely as follows

$$x_i = \frac{1}{2\alpha + g_i} [2\alpha, g_i]^\top, \quad y_i = \frac{2\alpha + g_i}{4\alpha\beta + g_i^2} [2\beta, g_i]^\top, \qquad (29)$$

with

$$g_i = \sqrt{(b_i - a_i)^2 + 4\alpha\beta + b_i - a_i > 0}.$$

Note that for $i \neq j$, $g_i = g_j$ if and only if $x_i = x_j$, if and only if $A_i = A_j + \delta_i I$ for some $\delta_i \in \mathbb{R}$.

**Proposition 8.** In case (21), when $d = 2$,

$$c_{sM} = \sum_{i,j} \alpha_i q_{ij} \ln \left( \frac{4\alpha\beta + g_i g_j}{4\alpha\beta + g_j^2} \right) \le 0.$$

Moreover, either for all $j > 1$, there is $\delta_j$ such that $A_j = A_1 + \delta_j I$ and $\omega \mapsto \Lambda_M(\omega)$ is constant and in particular, $c_{sM} = 0$; or otherwise $c_{sM} < 0$.
*Proof.* By Equations (10) and (29), we get that

$$c_{sM} = \sum_{i,j} \alpha_i q_{ij} \ln \left( \frac{4\alpha\beta + g_i g_j}{4\alpha\beta + g_j^2} \frac{2\alpha + g_j}{2\alpha + g_i} \right).$$

Note that by invariance of $\alpha$ with respect to $Q$, we have for all $(\eta_i)_{1\le i\le N}$

$$\sum_{i,j} \alpha_i q_{i,j} \eta_i = \sum_{i,j} \alpha_i q_{i,j} \eta_j = \sum_i \alpha_i \eta_i \qquad (30)$$

This yields

$$c_{sM} = \sum_{i,j} \alpha_i q_{ij} \ln \left( \frac{4\alpha\beta + g_i g_j}{4\alpha\beta + g_j^2} \right).$$

Now, applying for all $i$ and $j$ Cauchy-Schwarz inequality to the two-dimensional vectors $(\sqrt{4\alpha\beta}, g_i)$ and $(\sqrt{4\alpha\beta}, g_j)$, we have

$$(4\alpha\beta + g_i g_j) \le (4\alpha\beta + g_i^2)^{1/2} (4\alpha\beta + g_j^2)^{1/2},$$

with equality if and only if $g_i = g_j$. Therefore,

$$\begin{align*}
\sum_{i,j} \alpha_i q_{ij} \ln (4\alpha\beta + g_i g_j) &\le \sum_{i,j} \alpha_i q_{ij} \ln ((4\alpha\beta + g_i^2)^{1/2} (4\alpha\beta + g_j^2)^{1/2}) \\
&= \frac{1}{2} \sum_{i,j} \alpha_i q_{ij} \ln (4\alpha\beta + g_i^2) + \frac{1}{2} \sum_{i,j} \alpha_i q_{ij} \ln (4\alpha\beta + g_j^2) \\
&= \sum_{i,j} \alpha_i q_{ij} \ln (4\alpha\beta + g_i^2),
\end{align*}$$

where the last equality comes from (30), and the first inequality is an equality if and only if $g_i = g_j$ for all $i$ and $j$. Hence, we get that $c_{sM} \le 0$, and $c_{sM} < 0$ except if all the $g_i$ are equal, that is, if and only if all the $x_i$ are equal. But, as explained in the proof of Proposition 5, this implies that the map $\omega \mapsto \Lambda_M(\omega)$ is constant. $\square$

# 8 Proof of Proposition 3

Because of the symmetries of the problem, the principal right eigenvector $x_i$ of $A_i = L + R_i$ is simply $x_1$ up to a circular relabelling of the states, i.e. $x_{i,k} = x_{1,k-i+1}$ for all $k \in [1, d]$ (recall that indexes are understood modulo $d$). Similarly, denoting by $\hat{y}_i$ the principal left eigenvector of $A_i$ normalized so that $\mathbb{1}_d^\top \hat{y}_i = 1$ (so that $y_i = (x_i^\top \hat{y}_i)^{-1} \hat{y}_i$), which is the right eigenvector of $R_i + L^T$ (similar to $R_i + L$ except that migration turns in the other direction), it is obtained from $x_i$ by $\hat{y}_{i,k} = x_{i,2i-k}$. Then, in (10), in the synchronized case,

$$c_{sM} = c_{sM}^{sync} := \frac{1}{d} \sum_{i=1}^{d} \ln(y_{i+1}^{\top} x_i) = \ln(y_2^{\top} x_1) = \ln(\hat{y}_2^{\top} x_1) - \ln(\hat{y}_1^{\top} x_1),$$

while, in the a-synchronized case,

$$c_{sM} = c_{sM}^{async} := \frac{1}{d} \sum_{i=1}^{d} \ln(y_{i-1}^{\top} x_i) = \ln(\hat{y}_d^{\top} x_1) - \ln(\hat{y}_1^{\top} x_1),$$
Denoting for brevity $\lambda = \lambda(A_1)$, $x_1$ solves

$$ (\beta - 1)x_{1,1} + x_{1,d} = \lambda x_{1,1}, \quad -x_{1,i} + x_{1,i-1} = \lambda x_{1,i} \quad \forall i \in [1, d], $$

which, writing $r = \lambda + 1$, is solved as

$$ x_{1,i} = (r - \beta)r^{d-i}x_{1,1} \quad \forall i \in [2, d] $$

with $x_{1,1}$ being fixed by the normalization. Moreover, the equation $x_{1,1} = rx_{1,2}$ shows that $r$ is the unique positive solution of

$$ (r - \beta)r^{d-1} = 1 \qquad (31) $$

In particular, $r > \beta$ and $\lambda = r - 1$ is always larger than $-1$, and has the same sign as $\beta$ (since the solution for $\beta = 0$ is $r = 1$, and the solution of this equation is an increasing function of $\beta$).

Using the symmetries previously mentioned,

$$ \begin{aligned} \hat{y}_2^\top x_1 &= (r - \beta)^2 x_{1,1}^2 \begin{pmatrix} r^{d-2} \\ r^{d-1} \\ 1 \\ \vdots \\ r^{d-4} \\ r^{d-3} \end{pmatrix} \cdot \begin{pmatrix} r^{d-1} \\ r^{d-2} \\ r^{d-3} \\ \vdots \\ r \\ 1 \end{pmatrix} &= (r - \beta)^2 x_{1,1}^2 [(d - 2)r^{d-3} + 2r^{2d-3}] \\ &= x_{1,1}^2 [(d - 2)(r - \beta)r^{-2} + 2r^{-1}], \end{aligned} $$

while

$$ \begin{aligned} \hat{y}_1^\top x_1 &= (r - \beta)^2 x_{1,1}^2 \begin{pmatrix} r^{d-1} \\ 1 \\ r \\ \vdots \\ r^{d-3} \\ r^{d-2} \end{pmatrix} \cdot \begin{pmatrix} r^{d-1} \\ r^{d-2} \\ r^{d-3} \\ \vdots \\ r \\ 1 \end{pmatrix} &= (r - \beta)^2 x_{1,1}^2 [(d - 1)r^{d-2} + r^{2d-2}] \\ &= x_{1,1}^2 [(d - 1)(r - \beta)r^{-1} + 1], \end{aligned} $$

and

$$ \begin{aligned} \hat{y}_d^\top x_1 &= (r - \beta)^2 x_{1,1}^2 \begin{pmatrix} 1 \\ r \\ r^2 \\ \vdots \\ r^{d-2} \\ r^{d-1} \end{pmatrix} \cdot \begin{pmatrix} r^{d-1} \\ r^{d-2} \\ r^{d-3} \\ \vdots \\ r \\ 1 \end{pmatrix} &= (r - \beta)^2 x_{1,1}^2 dr^{d-1} \\ &= x_{1,1}^2 (r - \beta) d. \end{aligned} $$
We have thus obtained that

$$
\begin{align*}
c_{sM}^{\text{sync}} = \ln \left( \frac{\hat{y}_2^\top x_1}{\hat{y}_1^\top x_1} \right) &= \ln \left( \frac{(d-2)(r-\beta) + 2r}{(d-1)(r-\beta)r + r^2} \right) \\
&= \ln \left( \frac{(d-2)r^{1-d} + 2r}{(d-1)r^{2-d} + r^2} \right) \\
&= \ln \left( \frac{(d-2) + 2r^d}{(d-1)r + r^{d+1}} \right).
\end{align*}
$$

This has the same sign as $f(r)$ where

$$
f(s) = (d-2) + 2s^d - (d-1)s - s^{d+1}.
$$

When $d=3$, we see that

$$
f(s) = -s^4 + 2s^3 - 2s + 1 = (1-s)^3(1+s),
$$

which for $s > 0$ has the same sign as $1-s$. As we saw, $1-r = -\lambda$ has the same sign as $-\beta$, and as a conclusion in the case $d=3$, the sign of $c_{sM}$ is the opposite of the sign of $\beta$.

Consider the case $d > 3$ (for $d = 2$, $L$ is symmetric). Computing

$$
\begin{align*}
f'(s) &= 2ds^{d-1} - (d-1) - (d+1)s^d \\
f''(s) &= 2d(d-1)s^{d-2} - d(d+1)s^{d-1},
\end{align*}
$$

we see that $f(1) = f'(1) = 0$, and $f''$ cancels exactly twice, at $0$ and $s_* = \frac{2(d-1)}{d+1} \in (1,2)$, and is positive on $(0, s^*)$ and negative on $(s^*, +\infty)$. In particular, $f$ is convex decreasing positive over $(0,1)$, increasing convex over $(1, s_*)$ and then it becomes strictly concave, being eventually decreasing and then going to $-\infty$. In particular, it admits exactly one root $r^*$ over $(1, \infty)$, with $r_* > s_*$, $f(s) > 0$ if $s \in (1, r_*)$ and $f(s) < 0$ if $s > r_*$. Since the solution $r_\beta$ of (31) is an increasing function of $\beta$ and goes to infinity as $\beta$ goes to infinity, we obtain that, for each $d$, there is a unique $\beta_d > 0$ such that $r_\beta = r_*$ (which is simply $\beta = r_* - r_*^{1-d}$). For $\beta < \beta_d$ (resp. >), $r_\beta < r_*$ (resp. >), hence $f(r_\beta) > 0$ (resp. <), so that $c_{sM} > 0$ (resp. < 0).

In the asynchronised case,

$$
\begin{align*}
c_{sM}^{\text{async}} = \ln \left( \frac{\hat{y}_2^\top x_1}{\hat{y}_1^\top x_1} \right) &= \ln \left( \frac{(r-\beta)dr}{(d-1)(r-\beta) + r} \right) \\
&= \ln \left( \frac{dr^{2-d}}{(d-1)r^{1-d} + r} \right) = \ln \left( \frac{dr}{d-1 + r^d} \right).
\end{align*}
$$

This has the same sign as $f(r)$ with $f(s) = -s^d + ds - (d-1)$. Since $f'(s) = d(1-s^{d-1})$, $f$ is increasing over $[0,1]$ and decreasing over $[1, +\infty)$, reaching its global maximum at $1$, with $f(1) = 0$. In particular, for $\beta \neq 0$, $r \neq 1$ and thus $f(r) < 0$ from which $c_{sM}^{\text{async}} < 0$.

The periodic case follows from the discussion at the end of Section 3.1

Acknowledgments

The research of P. Monmarché is supported by the projects SWIDIMS (ANR-20-CE40-0022) and CONVIVIAL-
ITY (ANR-23-CE40-0003) of the French National Research Agency. SJS was supported by the U.S. National
Science Foundation Grant DEB-2243076.
References

M. Benaïm and E. Strickler. Random switching between vector fields having a common zero. *Ann. Appl. Probab.*, 29(1):326-375, 02 2019. doi: 10.1214/18-AAP1418. URL https://doi.org/10.1214/18-AAP1418.

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. ISSN 1050-5164. doi: 10.1214/13-AAP924. URL https://doi.org/10.1214/13-AAP924.

M. Benaïm, C. Lobry, T. Sari, and E. Strickler. Untangling the role of temporal and spatial variations in persistence of populations. *Theoretical Population Biology*, 154:1-26, 2023.

M. Benaïm, C. Lobry, T. Sari, and É. Strickler. A note on the top Lyapunov exponent of linear cooperative systems. *arXiv e-prints*, art. arXiv:2302.05874, Feb. 2023. doi: 10.48550/arXiv.2302.05874.

M. Benaim, C. Lobry, T. Sari, and E. Strickler. When can a population spreading across sink habitats persist? *Journal of Mathematical Biology*, 88(2):1-56, 2024.

M. Benaïm, C. Lobry, T. Sari, and E. Strickler. Dispersal-induced growth or decay in a time-periodic environment. *arXiv preprint arXiv:2407.07553*, 2024.

M. S. Boyce, C. V. Haridas, C. T. Lee, and the NCEAS Stochastic Demography Working Group. Demography in an increasingly variable world. *Trends in Ecology and Evolution*, 21:141-148, 2006.

R. Bürger. *The mathematical theory of selection, recombination, and mutation*. Chichester: Wiley, 2000.

H. Caswell. *Matrix Population Models*. Sinauer, Sunderland, Massachusetts, 2001.

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*, art. arXiv:2112.07005, Dec. 2021.

M. Davis. *Markov Models and Optimization*. Monographs on Statistics and Applied Probability. Chapman and Hall, 1993.

D. L. DeAngelis. *Individual-based models and approaches in ecology: populations, communities and ecosystems*. CRC Press, 2018.

N. H. Du, A. Hening, D. H. Nguyen, and G. Yin. Dynamical systems under random perturbations with fast switching and slow diffusion: Hyperbolic equilibria and stable limit cycles. *J. Differential Equations*, 293:313-358, 2021. ISSN 0022-0396. doi: https://doi.org/10.1016/j.jde.2021.05.032. URL https://www.sciencedirect.com/science/article/pii/S0022039621003223.

R. Fay, S. Michler, J. Laesser, J. Jeanmonod, and M. Schaub. Can temporal covariation and autocorrelation in demographic rates affect population dynamics in a raptor species? *Ecology and Evolution*, 10(4):1959-1970, 2020.

I. Gorenstein, I. Wainer, F. S. Pausata, L. F. Prado, M. Khodri, and P. L. S. Dias. A 50-year cycle of sea surface temperature regulates decadal precipitation in the tropical and south atlantic region. *Communications Earth & Environment*, 4(1):427, 2023.

I. Hanski. *Metapopulation Ecology*. Oxford Series in Ecology and Evolution. Oxford University Press, 1999.
H. Hasegawa, N. Katsuta, Y. Muraki, U. Heimhofer, N. Ichinnorov, H. Asahi, H. Ando, K. Yamamoto, M. Murayama, T. Ohta, et al. Decadal–centennial-scale solar-linked climate variations and millennial-scale internal oscillations during the early cretaceous. *Scientific reports*, 12(1):21894, 2022.

A. Hastings and L. J. Gross. *Encyclopedia of theoretical ecology*. Number 4. Univ of California Press, 2012.

M. Haviv, Y. Ritov, and U. G. Rothblum. Taylor expansions of eigenvalues of perturbed matrices with applications to spectral radii of nonnegative matrices. *Linear Algebra Appl.*, 168:159–188, 1992. ISSN 0024-3795. doi: 10.1016/0024-3795(92)90293-J. URL https://doi.org/10.1016/0024-3795(92)90293-J.

H. W. Hethcote. The mathematics of infectious diseases. *SIAM Review*, 42(4):599–653, 2000. doi: 10.1137/S0036144500371907. URL http://link.aip.org/link/?SIR/42/599/1.

C. H. Hilde, M. Gamelon, B.-E. Sæther, J.-M. Gaillard, N. G. Yoccoz, and C. Pélabon. The demographic buffering hypothesis: evidence and challenges. *Trends in Ecology & Evolution*, 35(6):523–538, 2020.

R. A. Horn and C. R. Johnson. *Matrix analysis*. Cambridge University Press, Cambridge, second edition, 2013. ISBN 978-0-521-54823-6.

G. Katriel. Dispersal-induced growth in a time-periodic environment. *Journal of Mathematical Biology*, 85(3): 24, 2022.

R. Kon, Y. Saito, and Y. Takeuchi. Permanence of single-species stage-structured models. *Journal of Mathematical Biology*, 48:515–528, 2004.

N. Kortessis, G. Glass, A. Gonzalez, N. Ruktanonchai, M. Simon, B. Singer, and R. Holt. Neglected consequences of spatio-temporal heterogeneity and dispersal: Metapopulations, the inflationary effect, and real-world consequences for public health. 2023a.

N. Kortessis, G. Glass, A. Gonzalez, N. W. Ruktanonchai, M. W. Simon, B. Singer, and R. D. Holt. Neglected consequences of spatio-temporal heterogeneity and dispersal: Metapopulations, the inflationary effect, and real-world consequences for public health. *bioRxiv*, pages 2023-10, 2023b.

R. C. Lewontin and D. Cohen. On population growth in a randomly varying environment. *Proceedings of the National Academy of Sciences USA*, 62:1056–1060, 1969.

W. Mitkowski. Dynamical properties of metzler systems. *Bulletin of the Polish Academy of Sciences: Technical Sciences*, 56(4):309–312, 2008.

P. Monmarché and E. Strickler. Asymptotic expansion of the invariant measure for Markov-modulated ODEs at high frequency. *arXiv e-prints*, art. arXiv:2309.16464, Sept. 2023. doi: 10.48550/arXiv.2309.16464.

M. Paniw, A. Ozgul, and R. Salguero-Gómez. Interactive life-history traits predict sensitivity of plants and animals to temporal autocorrelation. *Ecology letters*, 21(2):275–286, 2018.

M. Roy, R. Holt, and M. Barfield. Temporal autocorrelation can enhance the persistence and abundance of metapopulations comprised of coupled sinks. *American Naturalist*, 166:246–261, 2005.

S. Schreiber. Interactive effects of temporal correlations, spatial heterogeneity, and dispersal on population persistence. *Proceedings of the Royal Society: Biological Sciences*, 277:1907–1914, 2010.
S. Schreiber. Partitioning the impacts of spatial-temporal variation in demography and dispersal on metapopulation growth rates. *bioRxiv*, pages 2023-11, 2023.

S. Tuljapurkar. *Population Dynamics in Variable Environments*. Springer-Verlag, New York, 1990.

S. Tuljapurkar and C. Haridas. Temporal autocorrelation and stochastic population growth. *Ecology Letters*, 9 (3):327-337, 2006.