Vincent Bansaye
email: vincent.bansaye@polytechnique.edu

Amaury Lambert
email: amaury.lambert@upmc.fr

New approaches of source-sink metapopulations decoupling the roles of demography and dispersal

Keywords:

Source-sink systems are metapopulations of habitat patches with different, and possibly temporally varying, habitat qualities, which are commonly used in ecology to study the fate of spatially extended natural populations.

We propose new techniques that allow to disentangle the respective contributions of demography and dispersal to the dynamics and fate of a single species in a source-sink metapopulation. Our approach is valid for a general class of stochastic, individual-based, stepping-stone models, with density-independent demography and dispersal, provided the metapopulation is finite or else enjoys some transitivity property.

We provide 1) a simple criterion of persistence, by studying the motion of a single random disperser until it returns to its initial position; 2) a joint characterization of the long-term growth rate and of the asymptotic occupancy frequencies of the ancestral lineage of a random survivor, by using large deviations theory. Both techniques yield formulae decoupling demography and dispersal, and can be adapted to the case of periodic or random environments, where habitat qualities are autocorrelated in space and possibly in time.

 metapopulations. 1 Key words.

Introduction 1.Ecological background

Stochastic models of population dynamics play a prominent role in epidemiology and in ecology [START_REF] Kot | Elements of mathematical ecology[END_REF], in predicting the fate of natural populations (persistence vs extinction, disease outbreak, species invasion, competitive coexistence,...), and in computing some related quantities of interest (extinction probability, long-term growth rate, mid-term equilibrium distribution, stable age distribution, parameter elasticities).

Detailing these models is indispensable to understand the effect on these predictions of internal or external characteristic features, like spatial structure, age structure, intra-and interspecific interactions, or environmental change. In particular, metapopulation models [START_REF] Hanski | Metapopulation biology: Ecology, genetics, and evolution[END_REF], where the spatial structure is explicit, are used to infer the processes which have shaped contemporary range distributions, to predict migration trends or invasion fronts in response to biotic or abiotic changes, to understand the evolution of dispersal, to design protected areas and natural reserves, etc.

When the landscape is heterogeneous in terms of habitat suitability, even the mere question of predicting persistence can be a complicated task, since persistence is the result of the intricate interplay between population growth in suitable habitats, population depletion in unsuitable habitats and of how dispersal connects different habitat patches. In ecology, metapopulation models where habitat suitability is spatially heterogeneous are commonly referred to as sourcesink systems [START_REF] Dias | Sources and sinks in population biology[END_REF][START_REF] Holt | Population dynamics in two-patch environments: some anomalous consequences of an optimal habitat distribution[END_REF][START_REF] Pulliam | Sources, sinks, and population regulation[END_REF]. Roughly speaking, even if the definition of sources and sinks have been subject to debate [START_REF] Pulliam | Sources, sinks, and population regulation[END_REF][START_REF] Runge | The role of local populations within a landscape context: defining and classifying sources and sinks[END_REF], sources designate habitat patches where the habitat is suitable enough for the population to persist in the absence of dispersal (fundamental niche), and sinks are habitat patches where the population would become extinct in the absence of dispersal, or from which mortality during dispersal is too high to compensate growth. Spatial heterogeneity can be due to biotic environmental variables (predation risk, resource availability) or to abiotic environmental variables, which can either be constant through time (altitude or depth, latitude) or variable through time (precipitation, moisture, irradiance, pH, salinity).

To study the persistence of a single species in a metapopulation, it is common to further assume that population dynamics are density-independent. This assumption does certainly not hold for all natural populations, but can at least be used for populations whose persistence is guaranteed whenever their abundance is large enough to make this approximation unrealistic. It is also particularly relevant when asking about the establishment success of a new variant arising in few copies (immigrants, genetic mutants, infectives).

The assumption of density-independence allows theoretical ecologists to make use of linear models: matrix population models [START_REF] Caswell | Matrix population models: Construction, analysis, and interpretation[END_REF] for deterministic dynamics, multitype branching processes [START_REF] Asmussen | Branching processes, volume 3 of Progress in Probability and Statistics[END_REF][START_REF] Athreya | Branching processes[END_REF][START_REF] Haccou | Branching processes: Variation, growth, and extinction of populations[END_REF][START_REF] Jagers | Branching processes with biological applications[END_REF] for stochastic dynamics. These models are parsimonious in the number of parameters, and the associated mathematical theory is extremely well developed. The extinction probability has a very simple power dependence upon initial population size and composition, and under suitable assumptions, conditional on long-term survival, the geographic distribution of the population stabilizes over time, whereas its overall abundance grows exponentially with an exponent called the Malthusian growth rate, or long-term growth rate, or simply growth rate. In addition, the stable geographic distribution and the long-term growth rate are solutions to a well-known spectral problem. Namely, the growth rate is the maximal eigenvalue of the mean offspring matrix (encompassing both demography and dispersal), and the stable distribution is an associated eigenvector [START_REF] Li | Applications of Perron-Frobenius theory to population dynamics[END_REF][START_REF] Seneta | Non-negative matrices and Markov chains[END_REF].

A lot of work has been dedicated to extend these results to more complicated situations, like infinite metapopulations [START_REF] Metz | How should we define fitness in structured metapopulation models? Including an application to the calculation of evolutionarily stable dispersal strategies[END_REF], or, as earlier stressed, because spatial heterogeneity can itself be time-variable, to multitype branching processes in random environment [START_REF] Athreya | On branching processes with random environments, I: Extinction probabilities[END_REF][START_REF] Athreya | Branching processes with random environments, II: Limit theorems[END_REF][START_REF] Benaïm | Persistence of structured populations in random environments[END_REF][START_REF] Haccou | Establishment probability in fluctuating environments: a branching process model[END_REF][START_REF] Haccou | Establishment success and extinction risk in autocorrelated environments[END_REF][START_REF] Kaplan | Some results about multidimensional branching processes with random environments[END_REF][START_REF] Tanny | Limit theorems for branching processes in a random environment[END_REF]. More ecologically-related work has investigated which dispersal strategies are more likely to persist in metapopulations with random environment [START_REF] Gonzalez | The inflationary effects of environmental fluctuations in sourcesink systems[END_REF][START_REF] Schreiber | Interactive effects of temporal correlations, spatial heterogeneity and dispersal on population persistence[END_REF][START_REF] Schreiber | Invasion dynamics in spatially heterogeneous environments[END_REF], which such metapopulations are more prone to persistence [START_REF] Bascompte | Patchy populations in stochastic environments: Critical number of patches for persistence[END_REF], and which introduction strategies are more successful (single large vs several small) [START_REF] Haccou | Establishment success and extinction risk in autocorrelated environments[END_REF][START_REF] Wilcox | Conservation strategy: The effects of fragmentation on extinction[END_REF]. Specific attention has been given to coupled sinks, that is, metapopulations where each habitat patch is a (time-averaged) sink, but where populations might still persist thanks to dispersal in sparse favourable periods [START_REF] Evans | Stochastic population growth in spatially heterogeneous environments[END_REF][START_REF] Jansen | Populations can persist in an environment consisting of sink habitats only[END_REF][START_REF] Roy | Temporal autocorrelation can enhance the persistence and abundance of metapopulations comprised of coupled sinks[END_REF][START_REF] Schreiber | Interactive effects of temporal correlations, spatial heterogeneity and dispersal on population persistence[END_REF].

Goals and outline of the paper

In the present paper, our aim is to develop new methods in order to disentangle the contributions of demography and dispersal to the dynamics and outcome of source-sink systems with possibly varying environment. We will be interested in criteria for global persistence and in the computation of the long-term growth rate, and of the occupation frequencies of long-lived lineages.

One of the main problems of the spectral approach to the study of metapopulations is that the computation of eigenvalues and eigenvectors is totally opaque to biological interpretation. In particular, the respective contributions of dispersal and demography to the value of the long-term growth rate are very hard, if not impossible in general, to disentangle. As regards the question of persistence, we could ask for an alternative criterion, equivalent to, but simpler than, the positivity of this growth rate, which would avoid computing directly this eigenvalue. Similarly as in [START_REF] Evans | Stochastic population growth in spatially heterogeneous environments[END_REF][START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF][START_REF] Krkošek | An R 0 theory for source-sink dynamics with application to Dreissena competition[END_REF][START_REF] Runge | The role of local populations within a landscape context: defining and classifying sources and sinks[END_REF][START_REF] Schreiber | Invasion dynamics in spatially heterogeneous environments[END_REF], we will first seek to provide such an alternative criterion.

For example, in (st)age-structured models, it is easy to compute the net reproductive number R 0 , which is the expected total progeny produced in the lifetime of a single individual. Then thanks to a simple renewal argument, the condition R 0 > 1 is seen to be equivalent to possible survival. More rigorously, the set of juvenile offspring of a focal juvenile ancestor forms what is called a stopping line, for which it is known that an extended branching property holds [START_REF] Chauvin | Sur la propriété de branchement[END_REF]. This idea of the next generation-stopping line has been adapted to the spatial context in [START_REF] Krkošek | An R 0 theory for source-sink dynamics with application to Dreissena competition[END_REF][START_REF] Runge | The role of local populations within a landscape context: defining and classifying sources and sinks[END_REF], but remains of limited applicability. In the first part of this work, the key idea is to use as an alternative stopping line the set of descendants of a focal ancestor who are the first to return to the ancestor patch. Then by the extended branching property, the population will persist with positive probability iff the expected number, say R, of individuals on the stopping line is larger than 1.

If, as we first assume, the dispersal scheme does not depend on the state of the environment, then R can be expressed separately in terms of the mean offspring numbers in each patch (and in each environmental state) and of the motion of a single random disperser. More specifically, a random disperser is a single walker on the metapopulation which follows the dispersal stochastic scheme. We denote by X n its position at time n, so that (X n ) is a Markov chain with transitions given by the dispersal matrix which will be denoted by D. In the case when the environment is constant, we let m i be the mean number of offspring begot in patch i, and we prove that

R = m 1 E T -1 n=1 m Xn ,
where T is the first time the random disperser returns to patch 1 (assumed to be the initial patch). The population persists with positive probability iff R > 1. This way, our formulae are seen to disentangle the effects of demography and dispersal. If all other habitat patches than patch 1 have the same mean offspring m, then the last equality specializes into

R = m 1 E m T -1 ,
where the expectation in the last display can now be seen as the probability generating function of the random variable T -1 evaluated at m. We will also compute this expectation in some special cases of interest. It is interesting to note that the formulae obtained in [START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF] by a totally different method (expanding principle minors of the mean offspring matrix minus the identity matrix) feature numerous multiplicative terms also evoking closed reproductive paths.

In a second part, we will use large deviations techniques to prove that the logarithm of the long-term growth rate ρ and the asymptotic fraction (ϕ i ) of time spent in each patch of the ancestral lineage of a random survivor, are given respectively by the maximum and the unique argmax of a functional R-I defined on the set F of frequencies indexed by the metapopulation, where R only depends on the reproduction/survival scheme and I only depends on the dispersal scheme. Our formulae are then seen to decouple once again demography and dispersal. Namely,

log(ρ) = sup{R(f ) -I(f ) : f ∈ F} = R(ϕ) -I(ϕ),
where R is a linear functional of frequencies only depending on the mean offspring numbers in each patch

R(f ) := i f i log(m i ),
and I is a (more complicated) functional which only depends on the dispersal matrix D

I(f ) := sup i f i log(v i /(vD) i ) : v ≫ 0 ,
where v ≫ 0 denotes a positive row vector, that is, v i > 0 for each i. We find that (ϕ i ) never equals the stationary distribution of the single disperser, except in the case when all habitat qualities are identical (i.e., m i = m for all i). We compute ρ and ϕ in the case of a fully mixing metapopulation, i.e., when the probability for an individual to migrate from patch i to patch j does not depend on i, a case also referred to as parent-independent migration.

Addressing those questions is much more difficult when the model is enriched with a variable environment affecting simultaneously all habitat qualities. We can nevertheless adapt our arguments to the case when the environment is periodic or given by an ergodic sequence of random variables. We illustrate our speculations with two-patch metapopulations and a twostate environment. We make computations for fully mixing metapopulations.

Finally, we prove that for both periodic and ergodic environments, we can find parameters for which there is possible survival in coupled sinks, a result which was previously known in the absence of demographic stochasticity for fully mixing [START_REF] Jansen | Populations can persist in an environment consisting of sink habitats only[END_REF] and even partially mixing [START_REF] Schreiber | Interactive effects of temporal correlations, spatial heterogeneity and dispersal on population persistence[END_REF][START_REF] Evans | Stochastic population growth in spatially heterogeneous environments[END_REF] metapopulations.

Finally, we extend naturally our approach to a wide class of infinite metapopulations, called finitely transitive, in the sense that they can be naturally built by connecting copies of a finite subgraph called motif.

The paper ends with a short discussion on the uses of and possible extensions to our method (other models, local vs global persistence).

Preliminaries

Model

We consider a stochastic, individual-based model of spatially structured population dynamics. The spatial structure is a metapopulation of patches that can be of different habitat qualities. We label by i = 1, . . . , K the patches so that the model can be described by a labeled finite graph with weighted oriented edges. Vertices represent the patches, an oriented edge from vertex i to vertex j bears a weight d ij equal to the probability of dispersal from patch i to patch j. We let D be the square matrix with generic element d ij and we call it dispersal matrix. When D has identical columns, we will speak of parent-independent migration, or of a fully mixing metapopulation.

We assume a simple asexual life cycle with discrete non-overlapping generations and no density-dependence. At each generation, as a net result of reproduction and survival (including survival to possible migration), all individuals of patch i, independently from one another, leave to the next generation a random number of individuals, called offspring, all distributed as some random variable N i . The mean per capita number of offspring in patch i will be denoted by m i = E(N i ). In the second part of the paper, we will also enrich the model with a variable environment affecting simultaneously all patches. When the environment is in state w, we will denote by m i (w) the mean offspring in a patch of type i.

Immediately after local growth, each individual from the new generation migrates independently, from patch i to patch j with probability d ij . Since we assumed that mortality during dispersal is encompassed in the growth phase (see Discussion), we have K j=1 d ij = 1 for all i, i.e., D is a stochastic matrix.

Reproduction, survival and dispersal probabilities are assumed not to depend on local densities. Thanks to this assumption of density-independence, and because the mean offspring numbers encompass migration-induced mortality, we have the following classification. If m i > 1, we say that patch i is a source, and if m i ≤ 1, we say that patch i is a sink.

It will be convenient to assume that m 1 ≥ m 2 ≥ • • • ≥ m K . In addition, the problem of persistence is more interesting in the case when m 1 > 1 ≥ m K . Indeed, even in the presence of sinks, the metapopulation might persist thanks to local growth on sources replenishing sinks by dispersal. The case m K ≤ . . . ≤ m 1 ≤ 1 (resp. 1 < m K ≤ . . . ≤ m 1 ) leads trivially to extinction (resp. to persistence with positive probability).

Two natural examples with one source type and one sink type

Let us describe two examples with two possible habitat qualities, one source type and one sink type, that will be treated as a special case throughout the paper:

M := m 1 > 1, m := m 2 = • • • = m K ≤ 1.
First, we will be interested in the simple case with two patches, patch 1 with mean offspring M , and patch 2 with mean offspring m. In this case, we will always use the simplified notation p = d 12 and q = d 21 (see Figure 1). Second, we will consider the case when each source is only connected to sinks and two adjacent sources are separated by an array of n identical sinks. An example of such graph is the cyclic finite graph with one source and n sinks, or two sources connected by n sinks, or an infinite array with period n (see Figure 2)...

Method

The number of individuals located in patch i in generation n is denoted by

Z (i) n . The process Z = (Z (i) n , i = 1, • • • , K, n ≥ 0) is a multitype Galton-Watson process.
It is known from the mathematical literature [START_REF] Asmussen | Branching processes, volume 3 of Progress in Probability and Statistics[END_REF][START_REF] Athreya | Branching processes[END_REF] that either the population becomes extinct or it grows exponentially (under Assumption (A1) below). More specifically, we see that m i d ij is equal to the mean number of offspring of an individual living in patch i which will land into patch j in one time step, and therefore we call mean offspring matrix the matrix A defined as

A := (m i d ij : i, j = 1, • • • , K).
The maximal eigenvalue (see e.g. [START_REF] Seneta | Non-negative matrices and Markov chains[END_REF]) of A is the long-term growth rate, or simply growth rate of the metapopulation. Indeed [START_REF] Asmussen | Branching processes, volume 3 of Progress in Probability and Statistics[END_REF][START_REF] Athreya | Branching processes[END_REF], if ρ ≤ 1, the metapopulation dies out with probability 1, Notation Interpretation vertex i patch oriented edge probability of dispersal from patch i to patch j with weight d ij m i mean number of offspring in a patch with habitat type i M = m 1 mean growth rate in a source habitat (case

m 1 > 1 ≥ m 2 = • • • = m K ) m = m 2
mean growth rate in a sink habitat (case

m 1 > 1 ≥ m 2 = • • • = m K ) p = d 12
probability of dispersal from the source to the sink (case of 2 patches) q = d 21 probability of dispersal from the sink to the source (case of 2 patches) X random walk on the graph following the dispersal probabilities

Z n ρ n n→∞ -→ W,
where W is a component-wise non-negative and finite random vector (under Assumption (A2') below).

In the case of a fully mixing metapopulation, all column vectors of D are identical to some vector δ, say. Then A is the rank 1 matrix A = δµ, where µ is the row vector µ := (m 1 , . . . , m K ). In this case, the spectral approach is straightforward, since A n = (µδ) n-1 A, so that ρ = µδ = K j=1 m j δ j , and δ and µ are respectively right and left eigenvectors of A associated with ρ. We call random disperser a single individual who moves on the graph at discrete time steps following the dispersal probabilities. In other words, if X n denotes the position of such a random disperser after n time steps, then (X n ) is the Markov chain with transition matrix D

P(X n+1 = j | X n = i) = d ij .
The goal of this paper is to display new persistence criteria, along with results regarding the asymptotic growth rate and the asymptotic fraction of time spent in each patch (by an individual taken at random in the surviving population). In contrast with the method involving the maximal eigenvalue of the mean offspring matrix, this one can yield quite simple, interpretable and partially explicit criteria. In addition, these criteria decouple the contributions of dispersal and demography on population survival. In a number of remarks, we will also provide sufficient conditions for survival which are explicit, in particular in the case of a fully mixing metapopulation.

This approach is still valid when the graph is an infinite graph which can be reduced to a finite graph by transitivity. It is then called finitely transitive. The associated finite graph is called a motif, which is repeated to obtain the whole graph in such a way that the graph seen from any motif looks the same. A practical example is given by sources with the same quality connected by corridors of identical sinks and of the same length (see Figure 2 for an example). A finite-transitive graph could also be an (infinite) chessboard where whites are sinks and blacks are sources, the square lattice Z 2 where sources have coordinates of type (n, n) (diagonal) or of type (n, 0) (horizontal array), and so on (see Figures 3 and4

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q a) b)

Figure 3: Two examples of finite-transitive graphs ; a) the chessboard; b) a square grid where four-degree vertices are sources separated by n sinks (here n = 1).

Assumptions

Here, we list the assumptions we will use throughout the paper.

(A1) For convenience, we do not consider the degenerate case when any individual in patch 1 leaves one single offspring a.s. Thus, we assume that P(N 1 = 1) < 1.

(A2) We assume that in each habitat i the offspring number N i has a finite first moment, i.e., E(N i ) < ∞. It ensures that the total size of the population has a finite first moment in every generation.

(A2') We assume that the offspring number N i of any individual living in patch i satisfies

E(N i log + N i ) < ∞ (finite N log N moment).
(A3) To get the growth rate of the population, we will need to assume that m i = 0 for all i = 1 . . . K.

(A4) Irreducibility. For any ordered pair (i, j) of habitat patches, there is an integer n, such that P(X n = j | X 0 = i) = 0. That is, the random disperser can go in finite time from any patch to any other patch by using edges with positive weight.

(A5) Aperiodicity. For any patch i, the greatest common divisor (GCD) of the set of times n such that P(X n = i | X 0 = i) = 0, is equal to one. As a simple example, we mention the case when the graph is irreducible and at least one loop-edge has a positive weight, that is, there is at least one patch in which the probability of staying put is nonzero.

When the graph is both irreducible and aperiodic, we say that it is strongly irreducible or primitive. Then the Markov chain (X n ; n ≥ 0) is strongly irreducible and both the matrices D (and A if (A3) holds) are strongly irreducible, which is equivalent to the existence of n 0 > 0 such that all the coefficient values of M n 0 are positive [START_REF] Seneta | Non-negative matrices and Markov chains[END_REF].

Assumption (A2') ensures the convergence of Z n /ρ n to a non degenerate r.v. W which is non-negative on the survival event. This r.v. has only positive components if (A3) and (A4) ? ? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? ?

? --

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q a) b)

Figure 4: Two examples of graphs embedded in Z 2 which have an infinite motif (they still enjoy transitivity, but not finite) ; a) a diagonal of sources; b) a horizontal array of sources.

are also in force.

A first result on global persistence 3.1 General case

We now give a criterion for metapopulation persistence in terms of the random disperser X.

For that purpose, we assume from now on that the random disperser starts in patch 1 (X 0 = 1) and we denote by T the first return time of the random disperser into patch 1,

T := min{n ≥ 1 : X n = 1}.
Theorem 3.1 We assume (A1, A2, A4). Then the population persists with positive probability iff

m 1 E T -1 n=1 m Xn > 1.
In the case of a fully mixing metapopulation, D has all its columns equal to some column vector δ and it is known that ρ = K j=1 δ j m j (see previous section). It is easy to see that in this case, (X n ) is a sequence of i.i.d. random variables whose common distribution is given by δ, so that T is geometrically distributed with success parameter δ 1 , and

m 1 E T -1 n=1 m Xn = m 1 n≥0 δ 1 (1 -δ 1 ) n K j=2 δ j m j 1 -δ 1 n = δ 1 m 1 1 -K j=2 δ j m j
which is larger than 1 iff K j=1 δ j m j is larger than 1. Thus, we recover the criterion obtained with the spectral approach.

Remark 1 Observe that the expression given in the theorem can also be expressed as

T -1 n=1 m Xn = K i=2 m S T -1 (i) i , where S n (i) := #{1 ≤ k ≤ n : X k = i} is the time spent in habitat i by time n, that is S T -1 (i)
is the number of times the random disperser has visited patches of habitat type i strictly before time T . The advantage of this alternative formulation is that it carries over to models expressed in continuous time.

Let us prove this result and then give more challenging applications.

Proof. Let a (ancestor) be some individual placed at time 0 in patch 1. Define Y 1 as the number of offspring of a staying put in patch 1. Now for any integer n ≥ 2, let Y n denote the number of descendants of a at generation n living in patch 1 and whose ancestors at generations 1, 2, . . . , n -1 have all lived outside patch 1. Then set

Y := n≥1 Y n ,
that can be seen as the total number of descendants of a who live in patch 1 for the first time in their lineage (except a).

In the theory of random trees, this set of individuals belonging to one of the Y n individuals for some n, is called a stopping line. It is known [START_REF] Chauvin | Sur la propriété de branchement[END_REF] that a stopping line enjoys the extended branching property, in the sense that all the subtrees descending from distinct elements of a stopping line are i.i.d. copies of the tree (conditional on their types). Then the total numbers of descendants of each of the individuals of this stopping line who live in patch 1 for the first time in their lineage, are independent and all follow the same law as Y . In addition, any individual in the tree is either an ancestor or a descendant of some element of the stopping line. Therefore, the total number of descendants of a living in patch 1 is finite iff the branching process with offspring number distributed as Y is finite, which is equivalent to the a.s. extinction of Y . The bottomline is that the population persists in habitat 1 with positive probability iff E(Y ) > 1. Indeed, we have excluded the critical case when P(N 1 = 1) = 1. But this local persistence in habitat 1 is equivalent to the global persistence since the graph is irreducible.

Let us then compute E(Y ) to conclude. We first note that for every i = 1, . . . , K,

E(Z (i) n ) = K j=1 E(Z (j) n-1 )m j d ji ,
where

Z (i)
n denotes the number of individuals located in patch i at generation n. We prove easily by induction that the number of individuals

Y (i) n in patch i at generation n which have avoided patch 1 at generations k = 1, 2, • • • , n -1 satisfies E(Y (i) n ) = K j=2 E(Y (j) n-1 )m j d ji = j 1 ,...,j n-1 ∈{2,...K} d 1j 1 d j 1 j 2 . . . d j n-2 j n-1 d j n-1 i m 1 m j 1 • • • m j n-2 m j n-1 . As Y n = Y (1)
n , we get

E(Y n ) = j 1 ,...,j n-1 ∈{2,...K} d 1j 1 d j 1 j 2 . . . d j n-2 j n-1 d j n-1 1 m 1 m j 1 m j 2 • • • m j n-2 m j n-1 = m 1 E(1 T =n m X 1 m X 2 • • • m X n-2 m X n-1 ). Adding that Y = n≥1 Y (1)
n , we have

E(Y ) = m 1 n≥1 E(1 T =n m X 1 • • • m X n-1 ) = m 1 E(m X 1 • • • m X T -1 ).
This yields the result. 2

Case of two habitat types

Let us focus now on the special case when there are 2 habitat types and the source is solely connected to sinks:

M := m 1 > 1, m := m 2 = • • • = m K < 1.
We denote by

p = K j=2 d 1j
the probability of dispersing for an individual living in patch 1. The per capita mean offspring number sent out from a source at each generation is M p. Let σ be the time of first visit of a sink by the random disperser σ := inf{n ≥ 0 : X n = 1}, so that σ is a geometric random variable with success probability p. Next, let S denote the waiting time (after σ) before the random disperser visits a source (this source might or might not be the initial source patch X 0 )

S := inf{n ≥ 0 : X σ+n = 1}.
The duration S can be seen as the time spent in sinks between two consecutive visits of sources. By using the first transition of the random disperser, we get

E T -1 i=1 m Xn = 1 -p + p E m S ,
so that the previous theorem reads as follows.

Proposition 3.2 We assume (A1, A2, A4). Then the population persists with positive probability iff

M (1 -p) + eM p > 1, ( 1
)
where e is the depleting rate due to the sink habitat in the graph, defined as

e := E m S = k≥1 m k P(S = k).
Remark 2 If the average time spent in sinks has

E(S) < M -1 M p(1 -m) ,
then the population persists with positive probability. Indeed, the mapping f :

x → E(x S ) is convex so e = f (m) ≥ 1 + f ′ (1)(m -1) = 1 -(1 -m)E(S) > 1 - M -1 M p = 1 -M (1 -p) M p , which yields eM p + M (1 -p) > 1.
Let us check, in the simple case when there are only one source and one sink (K = 2 vertices), that criterion (1) is equivalent to the condition that the maximal eigenvalue ρ of A exceeds unity. Here the mean offspring matrix A is

A = M (1 -p) M p mq m(1 -q) .
The characteristic polynomial C of this square matrix is

C(x) = (M (1 -p) -x)(m(1 -q) -x) -M mpq.
Either M (1 -p) > 1 and the population living in the source ensures the persistence. Or M (1 -p) ≤ 1 and the quadratic polynomial is convex and has non negative derivative at 1. Thus, its leading eigenvalue is greater than 1 iff C(1) < 0, which reads

M p 1 -M (1 -p) > 1 -m(1 -q) mq .
We recover (1) since here S is geometric with success probability q, which yields

e = k≥1 q(1 -q) k-1 m k = mq 1 -m(1 -q) .
Notice that even in this simple case where A is a 2 × 2 matrix, the computation of the leading eigenvalue is cumbersome, and we have used a trick to explicitly specify the persistence criterion.

Example with pipes of identical sinks

Assume that the source is a vertex of degree 2 in the graph, connected to a left sink and a right sink. The probability of staying put on a source is still 1 -p, the probability of dispersing onto a left sink is pL, and the probability of dispersing onto a right sink is pR (so that L + R = 1). The sinks form a pipeline of n adjacent sinks linking adjacent sources. The probability of staying put on a sink is always s, the probability of dispersing from a sink onto one of its two neighboring sinks is r in the left-to-right direction of the pipe, and l in the right-to-left direction of the pipe (so that q = l + r = 1 -s). See Figure 5 for an example. This example will be directly extended in the last Section to infinite graphs, where pipelines of n sinks periodically connect sources (see Figure 6). We can compute exactly the depleting rate e of the above kind of the corresponding infinite graph (finitely transitive). Let λ > 1 > µ be the two ordered solutions to

mrx 2 -(1 -ms)x + ml = 0 x ≥ 0.
Then λµ = l r and λ + µ = 1 -ms mr .

Proposition 3.3

The depleting rate e is equal to

e = λ n -µ n λ n+1 -µ n+1 (L + Rλµ) + λ -µ λ n+1 -µ n+1 (R + L(λµ) n ).
Remark 3 In the two-patch case case (n = 1), we recover e = (1 -s)m 1 -ms .

In the case when dispersal is isotropic (l = r), we get

e = λ n -µ n + λ -µ λ n+1 -µ n+1 .
Notice that in both previous cases, the depleting rate does not depend on L or R. In the case of one single source and a large number of sinks (n → ∞), we get e = Lλ -1 + Rµ.

In the case of one single source and isotropic displacement, we then get e = λ -1 = µ.

Proof. Consider a random walk Y on {0, 1, . . . , n + 1}, with displacement at each time step being -1 with probability l, 0 with probability s, and +1 with probability r. Let T i denote the first hitting time of i by Y , and set T := min(T 0 , T n+1 ) as well as

a k := E(m T | Y 0 = k).
Returning to the random disperser X on the graph, it is easily seen that

E(m S | X σ is a left neighbour) = a 1 , while E(m S | X σ is a right neighbour) = a n , so that e = La 1 + Ra n .
Computations of a 1 and a n rely on the following recurrence relationship

a k = msa k + mla k-1 + mra k+1 k ∈ {1, . . . , n},
with boundary conditions a 0 = a n+1 = 1. This relation is obtained easily by considering the first transition of the walk X. 2

4 Growth rate and habitat occupation frequencies

General case

If u denotes an individual in generation n, we define H k (u) as the patch occupied by the ancestor of u in generation k ≤ n. Then, for every i ∈ {1, . . . , K},

F i (u) := 1 n #{0 ≤ k ≤ n : H k (u) = i}
is the occupancy frequency of patch i by the ancestral line of u. We further denote by U n an individual chosen randomly in the surviving population at generation n. We will see that the dispersal history of U n , as described by (F i (U n ) : i = 1, . . . , K) can be very different from that of a random disperser, since the ancestors of surviving individuals have better chance of having spent more time in sources than in sinks.

It is known that the growth rate ρ of the metapopulation is equal to the maximal eigenvalue of A (see [START_REF] Asmussen | Branching processes, volume 3 of Progress in Probability and Statistics[END_REF][START_REF] Athreya | Branching processes[END_REF][START_REF] Kurtz | A conceptual proof of the Kesten-Stigum theorem for multi-type branching processes[END_REF]). Moreover, the asymptotic occupancy frequencies of U n are deterministic and can be expressed as the product of the right and left maximal eigenvectors associated to ρ. We refer to [START_REF] Jagers | Stabilities and instabilities in population dynamics[END_REF][START_REF] Jagers | The asymptotic composition of supercritical multi-type branching populations[END_REF][START_REF] Hermisson | Mutation-selection balance: Ancestry, load, and maximum principle[END_REF], to Theorem 3.1 and 3.2 in [START_REF] Georgii | Supercritical multitype branching processes: the ancestral types of typical individuals[END_REF] in continuous time and to [START_REF] Kurtz | A conceptual proof of the Kesten-Stigum theorem for multi-type branching processes[END_REF] in discrete time. In this section, we want to give an alternative characterization to these quantities in terms of the random disperser and show an application.

To that purpose, we use the transition matrix D of the random disperser X on the graph. We denote by F i (X n ) the occupancy frequency of patch i by the random disperser X by time n

F i (X n ) := 1 n #{0 ≤ k ≤ n : X k = i}.
By assumptions (A4, A5), the random disperser has a stationary probability on V (stochastic equilibrium), that we denote by u = (u i : i = 1, . . . , K), which is the unique positive solution to uD = u. By the ergodic theorem, we also know that with probability 1,

F i (X n ) n→∞ -→ u i .
A typical single disperser will therefore occupy patch i with asymptotic frequency u i . This may not be the case of the ancestors of surviving individuals, whose paths must have favoured source patches. There is a trade-off between the pay-off in terms of fitness, gained by visiting source patches, and the cost in terms of likelihood, paid by deviating from the typical dispersal behavior. This trade-off is particularly obvious if we consider the case of a perfectly unsuitable habitat patch where the mean offspring is zero. In this case, the path followed by the ancestors of a surviving individual will necessarily have avoided this patch. Nevertheless, the asymptotic occupancy of this patch by a random disperser must be nonzero by the irreducibility assumption (A4).

There is a way of quantifying both the cost and pay-off of deviating from the typical dispersal behavior, that is, of having asymptotic occupancy frequencies f = (f i : i = 1, . . . , K), where f is a given element of the set F of non-negative frequencies on the graph

F :=    f = (f i : i = 1, . . . , K) : f i ≥ 0, i=1,...,K f i = 1    .
First, the probability that a random disperser has occupancy frequencies close to some given f by time n decreases exponentially with n at rate I(f ), which can thus be interpreted as the cost of the f -occupancy scheme:

I(f ) := sup    i=1,...,K f i log(v i /(vD) i ) : v ≫ 0    , (2)
where v ≫ 0 denotes a positive row vector, that is, v i > 0 for each i = 1, . . . , K. Indeed, large deviations theory [START_REF] Dembo | Large deviations techniques and applications[END_REF][START_REF] Hollander | Large deviations[END_REF] ensures that for any ǫ ≪ 1, as n → ∞, we have

P(f i -ǫ ≤ F i (X n ) ≤ f i + ǫ for all i = 1, • • • , K) ≍ exp(-nI(f )). (3)
We refer the reader to Section A for a more rigorous formulation. Taking v = u in (2) shows that I(f ) is of course always non-negative. This function is also convex. When f = u, one can easily check that each partial derivative of I is zero and it can be proved that the supremum in ( 2) is attained for v = u, so that I(u) = 0. This was indeed expected, since f = u is the natural occupancy scheme of the random disperser. Second, the reproductive pay-off of f can be defined as the fitness of a non-random disperser with given f -occupancy scheme, that is

R(f ) := i=1,...,K f i log(m i ).
Indeed, the total size of a population of individuals all adopting this dispersal behavior can be seen to grow like

K i=1 m nf i i = exp(nR(f )).
Thus the cost (in terms of likelihood) for a population to follow some occupancy scheme is quantified by I and the reproductive pay off by R. The best strategy (regarding the growth of the population) is to have an asymptotic occupancy frequency ϕ which maximizes the difference R -I. If this optimum is positive, then the population survives with positive probability. In addition, the ancestral line of a randomly chosen surviving individual will have visited patch i with frequency ϕ i . These results are stated below. The last assertion indicates that this optimal occupancy scheme ϕ = (ϕ i , i = 1 . . . , K) is always different from the natural occupancy scheme u of one single random disperser, except when all habitat types have the same quality.

Theorem 4.1 We assume (A2', A3, A4, A5). The growth rate ρ of the metapopulation is given by

log(ρ) = max {R(f ) -I(f ) : f ∈ F} .
In addition, if ρ > 0, for any patch i = 1, . . . , K, conditional on the population being alive at time n, the occupancy frequency of patch i by the ancestral line of a randomly chosen individual U n in the surviving population at time n, converges to ϕ i in probability :

F i (U n ) n→∞ -→ ϕ i ,
where the frequency vector ϕ = (ϕ i ) ∈ F is uniquely characterized by

log(ρ) = R(ϕ) -I(ϕ).
The occupancy frequency ϕ coincides with the stationary distribution u of X (if and) only if

m 1 = m 2 = . . . = m K .
In the same vein, we refer to Theorem 3.3 in [START_REF] Georgii | Supercritical multitype branching processes: the ancestral types of typical individuals[END_REF] for a description of the lineage of surviving individuals for multitype branching processes in continuous time. The irreducibility assumption (A4) is required to use Sanov's large deviation theorem. The first result is actually a consequence of the functional version of Sanov's theorem. The aperiodicity assumption (A5) and the exclusion of the degenerated cases by (A3) are used for the two additional results. The assumption (A2) would be enough for the first part but (A2') is required for the additional results.

The proof is deferred to Section A.

The fully mixing case

In the case of a fully mixing metapopulation, D has all its columns equal to some column vector δ and we have already seen that ρ = K j=1 δ j m j . It is also easy to see that δ is the stable geographic distribution of the population. As a first observation, note that here the stable geographic distribution δ is also the stationary distribution of the random disperser. Then the last part of the previous theorem ensures that, except when all m j 's are equal, the stable geographic distribution is different from the occupation distribution ϕ of random long-lived lineages.

We now use the approach developed in the previous subsection for an alternative computation of ρ. We will also determine the occupation frequency ϕ of ancestral lineages. We first compute the functional I. Here, for any row vector v ≫ 0,

K i=1 f i log(v i /(vD) i ) = K i=1 f i log(v i /|v|δ i ),
where |v| = K j=1 v j . Then, differentiating this functional with respect to v j yields (f j /v j ) -1/|v|. As a result,

I(f ) = K i=1 f i log(f i /δ i ).
Then differentiating R -I with respect to f j (j = 1 for example, and

f 1 = 1 -K j=2 f j ) yields log(δ j m j ) -1 -log(f j ). The bottomline is ϕ j = δ j m j K i=1 δ i m i .
Plugging this as the argument of R -I yields

log(ρ) = R(ϕ) -I(ϕ) = log K i=1 δ i m i ,
which was the expected result.

5 Fluctuating environments

General setting

We now enrich our model with a fluctuating environment. The environment is embodied by a certain value w which belongs to a finite set of states. We assume that the environment affects simultaneously all patches, but not necessarily in the same way. We keep on assuming a simple asexual life cycle with discrete non-overlapping generations and no density-dependence. Now the environment is assumed to affect reproduction and survival, but not dispersal. Specifically, at each time step, conditional on the state w of the environment, individuals reproduce independently according to some distribution which depends on the habitat type of their dwelling patch. We denote by m i (w) the mean offspring number of individuals dwelling in patch i when the environment is in state w.

Except in the last subsection, we will assume that the environment alternates periodically at each time step between two states (circadian cycle, seasons). Actually, the same method would allow to deal with any finite number of environmental states varying periodically.

We call e 1 and e 2 the two possible states of the environment, so now we have 2K habitat qualities m i (e j ), for i = 1, . . . , K and j = 1, 2. The Markov chain Z = (Z (i) n , i = 1, . . . , K; n ≥ 0) is no longer time-homogeneous and is called multitype branching process in varying environment [START_REF] Haccou | Branching processes: Variation, growth, and extinction of populations[END_REF]. However, restricting the observation of the metapopulation to times when the environment is in the same state allows to adapt the arguments of the previous section. Indeed, (Z (i) 2n , i = 1, . . . , K; n ≥ 0) is a multitype branching process with mean offspring matrix A with generic element a ij = k=1,...,K m i (e 1 )d ik m k (e 2 )d kj i, j = 1, . . . , K.

This amounts to changing the stopping line of the previous section, which was made of descendants returning to the ancestor patch for the first time in their lineage, for the stopping line of descendants returning to the ancestor patch for the first even time in their lineage. In the following subsection, we treat the case of two patches and determine the global persistence criterion. We then handle the general case using the random disperser.

Example with two patches and two periodic environments

For convenience, even if the environment is now variable, the two patches are still called respectively the source (patch 1) and the sink (patch 2). The mean number of offspring in the source are denoted by M 1 = m 1 (e 1 ) and M 2 = m 1 (e 2 ). In the sink, they are denoted by m 1 = m 2 (e 1 ) and m 2 = m 2 (e 2 ).

Theorem 5.1 A necessary and sufficient condition for global persistence is

M 1 M 2 (1 -p) 2 + (M 1 m 2 + m 1 M 2 )pq + m 1 m 2 (1 -q) 2 > min 2, 1 + M 1 M 2 m 1 m 2 (1 -p -q) 2 .
Remark 4 It is easy to find examples where both patches are sinks on average but the metapopulation survives with positive probability thanks to dispersal. Indeed, each patch is a sink if (and only if ) M 1 M 2 ≤ 1 and m 1 m 2 ≤ 1. Assuming for example that p = q = 1/2 and m 1 = m 2 = m, the global survival criterion becomes

M 1 M 2 + m(M 1 + M 2 ) + m 2 > 4, which holds as soon as m(M 1 + M 2 ) > 4.
Proof. The mean offspring matrix of (Z 2n ; n ≥ 0) is given by

a(1, 1) = M 1 M 2 (1 -p) 2 + M 1 m 2 pq a(1, 2) = M 1 m 2 p(1 -q) + M 1 M 2 (1 -p)p a(2, 1) = m 1 M 2 q(1 -p) + m 1 m 2 (1 -q)q a(2, 2) = m 1 m 2 (1 -q) 2 + m 1 M 2 qp.
The maximum eigenvector of the matrix A = (a(i, j) :

1 ≤ i, j ≤ 2) is the largest root of the polynomial x 2 -(a(1, 1) + a(2, 2))x + a(1, 1)a(2, 2) -a(1, 2)a(2, 1).
So it is less than 1 iff

a(1, 1) + a(2, 2) + (a(1, 1) -a(2, 2)) 2 + 4a(1, 2)a(2, 1) ≤ 2
Then the criterion for a.s. extinction of (Z 2n :

n ∈ N) is a(1, 1) + a(2, 2) ≤ 2 and (a(1, 1) -a(2, 2)) 2 + 4a(1, 2)a(2, 1) ≤ (2 -a(1, 1) -a(2, 2)) 2 .
The second inequality becomes a(1, 1) + a(2, 2) ≤ 1 + a(1, 1)a(2, 2) -a(1, 2)a(2, 1), which gives

M 1 M 2 (1 -p) 2 + M 1 m 2 pq + m 1 m 2 (1 -q) 2 + m 1 M 2 qp ≤ 1 + M 1 M 2 m 1 m 2 (1 -p -q) 2 .
This completes the proof. 2

Global persistence for more than two patches

Here, we extend the previous result to the case of a general, finite graph. We want to state a global survival criterion which generalizes Theorem 3.1 to periodic environments. Assume again that the random disperser starts at time 0 in patch 1 and set T the first even time when the random disperser goes back to habitat 1

T := min{n ≥ 1 : X n = 1 and n is even}.

By a direct adaptation of the proof of Theorem 3.1 replacing Z n with Z 2n , we get the following statement.

Theorem 5.2 We assume (A1) holds for at least one environment, (A2) holds for both environments and (A4) holds. Then the population persists with positive probability iff

m 1 E T -1 n=1 m Xn (w n ) > 1,
where the sequence (w n : n ≥ 1) can take one of the two values (e 1 , e 2 , e 1 , . . .) or (e 2 , e 1 , e 2 , . . .), depending whether the initial environment is e 1 or e 2 .

Rate of growth and habitat occupation frequency

The generalization to periodic environments of the results of the previous section can be achieved by changing the state-space {1, . . . , K} of the random disperser to the state-space of oriented edges of the graph, i.e., ordered pairs of vertices

E := {1, . . . , K} 2 .
Denote by B the transition matrix of the Markov chain (X 2n , X 2n+1 ; n ≥ 0), which indeed takes values in E. Then denote by F the set of frequencies indexed by E

F := (f E , E ∈ E) : f E ≥ 0, E∈E f E = 1 ,
and define the new cost function I : F → R as

I(f ) := sup E∈E f E log(v E /(vB) E ) : v ≫ 0 ,
where v denotes a non-negative vector indexed by E, such that v ≫ 0, that is, v E > 0 for all E ∈ E. Also define the new pay-off function R :

F → R as R(f ) := E=(i,j)∈E f E log(m i (e 1 )m j (e 2 )).
We can also provide an expression of I in terms of the entropy function using Theorem 3.1.13 in [START_REF] Dembo | Large deviations techniques and applications[END_REF]. The generalization of Theorem 4.1 can be stated as follows.

Theorem 5.3 We assume that (A2', A3) hold for both environments and (A4, A5) hold.

The growth rate ρ of the metapopulation is given by

2 log(ρ) = max{R(f ) -I(f ) : f ∈ F}.
In addition, for any patch i ∈ {1, . . . , K}, conditional on the population being alive at time n, the frequencies of occupation of patch i by the ancestral line of a randomly chosen individual U n in the surviving population at time n, converges in probability :

F j (U n ) n→∞ -→ i∈{1,...,K} ϕ i,j ,
where the vector (ϕ i,j : (i, j) ∈ E) is characterized by

2 log(ρ) = R(ϕ) -I(ϕ).
The proof follows that of Theorem 4.1, with now

E(|Z 2n+1 |) = E n k=0 m X 2k m X 2k+1 = E   K i=1 K j=1 m i (e 1 ) S (1) n (i) m j (e 2 ) S (2) n (j)   ,
where |Z n | is the total number of individuals in source patches at generation n and

S (1) n (i) = #{k ≤ n : X 2k = i}, S (2) n (i) = #{k ≤ n : X 2k+1 = i}.

Example with fully mixing patches

We extend the computations of the previous section to periodic environments. We focus on the fully mixing population :

B (i,j)(k,l) = d jk d kl = δ k δ l .
The mean matrix associated to the Galton Watson process Z 2n is

A ij = m i (e 1 )δ j K k=1 δ k m k (e 2 ) .
Thus the right and left eigenvectors are still given by δ and (m 1 (e 1 ), • • • m K (e 1 ) and the spectral approach given previously can be followed readily.

Let us now focus on the approach given in the last Theorem. As (vB

) E = |v|δ E 1 δ E 2 , the differentiation of E∈E f E log(v E /(vB) E ) = E∈E f E log(v E ) -log(|v E |δ E 1 δ E 2 )
with respect to v E yields the minimum. As in the previous section, we get

I(f ) = E∈E f E log f E δ E 1 δ E 2 . Then R(f ) -I(f ) = E∈E f E log (m E 1 (e 1 )m E 2 (e 2 )δ E 1 δ E 2 /f E ).
We arbitrarily choose E 0 ∈ E, so we can write

f E 0 = 1 - E∈E 0 f E
and make all other partial derivatives of R -I equal 0 when evaluated at ϕ. We get for every

E = E 0 : log (m E 1 (e 1 )m E 2 (e 2 )δ E 1 δ E 2 )+1-log(ϕ E )-log m E 0 1 (e 1 )m E 2 0 (e 2 )δ E 0 1 δ E 0 2 ) + 1 -log(ϕ(E 0 )) = 0,
which gives the habitat occupation frequencies

ϕ E = δ E 1 δ E 2 m E 1 (e 1 )m E 2 (e 2 ) E∈E δ E 1 δ E 2 m E 1 (e 1 )m E 2 (e 2 ) -1
.

We can also now deduce the growth rate ρ. We get

R(ϕ) -I(ϕ) = log( E∈E δ E 1 δ E 2 m E 1 (e 1 )m E 2 (e 2 ))
and

ρ = 1 2 [R(ϕ) -I(ϕ)] = K i=1 δ i m i (e 1 ) 1/2 K i=1 δ i m i (e 2 ) 1/2
. This is the same growth rate as the one computed in [START_REF] Jansen | Populations can persist in an environment consisting of sink habitats only[END_REF] for large populations (with two patches).

Some comments on random environments

A more natural way of modeling fluctuating environment in ecology is to assume random rather than periodic environment. The approach developed for periodic environments cannot be extended to random environments directly. Indeed, since the environment affects the whole metapopulation simultaneously, the randomness of environments correlates reproduction success in different patches. The process (Z

(i)
n , i = 1, . . . , K; n ≥ 0) counting the population sizes on each patch is now a multitype branching process in random environment (MBPRE) [START_REF] Athreya | On branching processes with random environments, I: Extinction probabilities[END_REF][START_REF] Athreya | Branching processes with random environments, II: Limit theorems[END_REF].

Let us denote by A(w) the mean offspring matrix (involving dispersal) in environment w. Specifically, the generic element a ij (w) of A(w) is the mean offspring number of a typical individual dwelling in patch i sent out to patch j by dispersal, when the environment is w, so that

A ij (w) = m i (w) d ij .
We will now assume that the state-space of environments is finite and that the sequence (w n : n ≥ 0) of environment states through time is a stationary, ergodic sequence, possibly autocorrelated, in the sense that the states need not be independent. Under this assumption, it is proved in [START_REF] Furstenberg | Products of random matrices[END_REF] (under the further assumption E(log + A(w 0 ) ) < ∞, where expectation is taken w.r.t. the environment) that the limit γ of the sequence

1 n log A(w n )A(w n-1 ) . . . A(w 0 )
exists with probability 1 and is deterministic, where . denotes the maximum row sum of the matrix. This is interesting to us because it is further shown in [START_REF] Athreya | On branching processes with random environments, I: Extinction probabilities[END_REF][START_REF] Kaplan | Some results about multidimensional branching processes with random environments[END_REF][START_REF] Tanny | Limit theorems for branching processes in a random environment[END_REF] (again under some further assumptions, see Section B), that the extinction criterion and the growth rate of this MBPRE are respectively given by the sign and the value of γ, more specifically, γ = log ρ.

Unfortunately, this does not give a very explicit condition for global persistence. But again using the random disperser, we can give some sufficient conditions for survival. For simplicity, we turn our attention to the example of two patches and two environments e 1 and e 2 . At any time step, the probability that the environment is in state e 1 is denoted by ν ∈ (0, 1) (so the probability that the environment is in state e 2 is 1-ν). We show again that the population may survive in sinks only. As in the case of periodic environments, the mean number of offspring in the first patch (the source) is denoted by M 1 = m 1 (e 1 ) and M 2 = m 2 (e 1 ). In the second patch (sink), they are denoted by m 1 = m 2 (e 1 ) and m 2 = m 2 (e 2 ).

For the sake of simplicity, we state the results for the special case when the sequence is a Markov chain. We denote by α the transition from e 1 to e 2 and by β the transition from e 2 to e 1 . Then it is well-known that ν = β/(α + β) is the asymptotic fraction of time spent in state e 1 . The case of independent environments is recovered when α + β = 1. Note that as soon as α + β = 1, the sequence of environment states is auto-correlated. Proposition 5. [START_REF] Athreya | Branching processes[END_REF] We have the following lower bound for the growth rate of the metapopulation.

log(ρ) ≥ ν log(M 1 ) + (1 -ν) log(m 2 ) + να log(pq) + ν(1 -α) log(1 -p) + (1 -ν)(1 -β) log(1 -q).
Remark 5 Observe that this lower bound does not depend on M 2 and m 1 . Again one can display examples where both patches are sinks but the metapopulation survives with positive probability in the presence of dispersal. Each patch is a sink if (and only if )

M ν 1 M 1-ν 2 ≤ 1 and m ν 1 m 1-ν 2 < 1.
Actually one can manage to keep γ > 0 while M ν 1 M 1-ν 2 < 1, m 1 < 1 and m 2 < 1, for example with M 2 small and M 1 large for some fixed p, q, m 1 , m 2 . This corresponds to e 2 being a catastrophic environment in the source patch but the population survives in patch 2 when a catastrophe occurs.

Proof. We consider only the subpopulation avoiding patch 1 when the environment is equal to e 2 . This means that this population reproduces with mean offspring number M 1 (1 -p) in patch 1 while the environment is e 1 . Each time the environment e 2 occurs, we consider the part of this population which has dispersed to patch 2. This corresponds to a mean offspring number of M 1 p. This population then stays in patch 2 and reproduces with mean offspring number m 2 (1 -q) until the environment is again equal to e 1 . We then consider the part of this population which goes back to patch 1. This corresponds to a mean offspring number of m 2 q.

Thus the patch of the ancestors of the individuals we keep is equal to 1 (resp. 2) if it lived in environment e 1 (resp. e 2 ). Then at time n, the mean size of the population we consider is equal to

M N 1 (n) 1 m N 2 (n) 2 (1 -p) N 11 (n) (1 -q) N 22 (n) p N 12 (n) q N 21 (n)
where

N i (n) (i ∈ {1, 2}
) is the number of times before generation n when the environment is equal to e i and N ij (n) (i, j ∈ {1, 2}) is the number of one-step transitions of the environment from e i to type e j until time n. By ergodicity, we know that these quantities have deterministic frequencies asymptotically. In the case of a Markovian sequence of environments, as n → ∞,

N 1 (n) ∼ νn, N 2 (n) ∼ (1 -ν)n,
and

N 11 (n) ∼ ν(1 -α)n, N 12 (n) ∼ ναn, N 21 (n) ∼ (1 -ν)βn, N 22 (n) = (1 -ν)(1 -β)n.
Using the growth of this particular part of the whole population directly gives us a lower bound for γ:

γ ≥ ν log M 1 +(1-ν) log m 2 +να log p+(1-ν)β log q+ν(1-α) log(1-p)+(1-ν)(1-β) log(1-q).
Noticing that (1 -ν)β = να completes the proof. 2

Remark 6 We could improve these results by considering more sophisticated strategies. For example, we could consider the subpopulation which stays in patch 1 if (and only if ) the number of consecutive catastrophes is less than k and then optimize over k.

Actually, the proof relies on a stochastic coupling. Roughly speaking, the subpopulation we consider avoids the bad patches at the bad times and follows a (one type) branching process in random environment e 11 , e 12 , e 21 and e 22 respectively with stationary probabilities ν(1 -α), να,

(1 -ν)β and (1 -ν)(1 -β)and mean offspring M 1 (1 -p), M 1 p, m 2 q and m 2 (1 -q).
Observe also that we can derive a lower bound using the permanent of the mean matrix M of the MBPRE from Proposition 2 in [START_REF] Benaïm | Persistence of structured populations in random environments[END_REF]. But this lower bound is not relevant for understanding the survival event in sinks only.

Metapopulation on infinite graphs

We now turn our attention to infinite graphs labeled by a countable set F . Each patch P ∈ F has a type i = (P ) ∈ N which gives its habitat quality, that is, the mean number of offspring in patch P is equal to m (P ) .

To generalize all the previous results, we require that the infinite graph has a finite motif. Let us first provide the reader with some examples satisfying this assumption, before giving rigorous definitions. These examples are chosen among source-sink metapopulations with two habitat qualities (one source type and one sink type).

• An infinite linear periodic array of patches (see Figure 2 for an example). The patches can then be labeled by P ∈ F = Z and the type of patch P is equal to (P ), for some integers N, K > 0 and a function  : Z → [1, K] such that (P + N ) = (P ) for every P ∈ Z. The motif is a line of lenght N with one source.

• The chessboard (see Figure 7). The motif is composed of one source and one sink.

• Star sources with 2d pipelines of sinks (see Figure 8 for d = 2). The motif is built by a source with d pipelines of sinks of the same length.

Let us now specify mathematically these definitions. The oriented edges from P to Q are weighed by d P Q . A mapping T of the graph is called an isomorphism if it conserves the types of the vertices as well as the weights of the oriented edges: T is a bijection of F such that for all P, Q ∈ F ,

d T (P )T (Q) = d P Q , (T (P )) = (P ).
The associated equivalence relation ∼ between the patches of the graph is defined by P ∼ P ′ iff there exists an isomorphism T of the graph such that T (P ) = P ′ .

The class of a patch P is defined as the equivalence class Cl(P ) = {P ′ : P ′ ∼ P }. Every patch of this class has the type of P . A graph for which there exists an isomorphism which is not the identity is called transitive.

The collection of the distinct classes (Cl(i) : i ∈ V ) of a transitive graph form a partition of the patches of the graph. Such subsets V of patches are called motifs. With a slight abuse of notation, the transition probabilities on a motif V are denoted by (d P Q : P ∈ V, Q ∈ V ) where

d P Q = Q ′ ∈Cl(Q) d P Q ′
is constant in the same equivalence class. ?

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q a) b)

Figure 7: a) The chessboard (with periodic arrow labels -not represented) is a graph that can be collapsed into: b) a two-vertex graph (loop edges are not represented).

The initial graph can be seen as a family of copies of a motif properly connected. Observe that not all graphs with a regular structure are finitely transitive. For example, the cases illustrated by Figure 3 can be collapsed into a finite motif, but those given in Figure 4 (a diagonal or a horizontal array of sources in Z 2 ) cannot. We can easily extend the results of Section 3 to infinite graphs whose motifs have a finite number of sources. Thus, we consider the random disperser X n on the graph which follows the dispersal probabilities (d P Q : P, Q ∈ F 2 ). To that purpose, we assume that the random disperser starts in a patch of type 1 ((X 0 ) = 1) and we denote by T the first return time of the random disperser into a (possibly different) patch of type 1.

T := min{n ≥ 1 : (X n ) = 1}. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q ? ? a) b) probability of dispersal from patch P to patch Q m i mean number of offspring in a patch with habitat type i Table 2: Notation. Theorem 6.1 Assume that the metapopulation graph is transitive and that there is a motif with a finite number of sources. Assume (A1, A2, A4) hold for this motif. Then the population persists with positive probability iff

m 1 E T -1 n=1 m (Xn) > 1.
We provide the idea of the proof. The number of individuals located in patch P in generation n is denoted by N n , P ∈ V ; n ≥ 0) is a multitype Galton-Watson process (with possibly an infinite number of types). We see that m (P ) d P Q is equal to the mean number of offspring of an individual living in patch P which will land into patch Q in one time step, and therefore we call mean offspring matrix A := (m (P

) d P Q : P, Q ∈ V ),
Then, everything happens as if the metapopulation evolves on a graph given by the motif. We follow then the proof of Theorem 4 to prove that the criterion m 1 E( T -1 n=1 m (Xn) ) > 1 is the persistence criterion in the source habitat 1. But using the facts that the number of sources in a motif is finite and that a motif is irreducible ensures that it is a global persistence criterion.

When the motif is finite, we get exactly the model described in the introduction and derive the following counterpart of Theorem 4.1. We denote now by D the transition matrix of the Markov chain (X n ), and still use for f = (f P : P ∈ V ) Theorem 6.2 Assume that the metapopulation graph is transitive with a finite motif. Assume (A1, A2', A3, A4, A5) hold for this motif. The growth rate ρ of the metapopulation is given by log(ρ) = max R(f ) -I(f ) :

I(f ) := sup P ∈V f P log(v P /(vD) P ) : v ≫ 0 ,
P ∈V f P = 1, f P ≥ 0 .
In addition, for any P ∈ V , conditional on the population being alive at time n, the occupancy frequency of a patch with type (P ) by the ancestral line of a randomly chosen individual U n in the surviving population at time n, converges to ϕ i in probability

F i (U n ) n→∞ -→ ϕ i ,
where the vector ϕ belongs to F and is characterized by

log(ρ) = R(ϕ) -I(ϕ).
The occupancy frequency ϕ coincides with the stationary distribution u of X only if m (P ) = m (P ′ ) for all patches P, P ′ .

The results obtained for the periodic and random environments can be derived similarly for transitive graphs.

Discussion

We have proposed two new techniques to derive simple criteria of persistence (by using the path of a random disperser) and to characterize the long-term growth rate of a single species in a source-sink metapopulation (by using large deviations for the path of the random disperser), as well as the occupancy frequencies of long-lived ancestral lineages. The expressions obtained thanks to these techniques decouple the contributions of reproduction/survival vs dispersal to the past, growth and persistence of these metapopulations. These techniques apply to a general class of stochastic, individual-based, stepping stone models of source-sink metapopulations, even when habitat quality is (not only variable in space but) variable through time, and even for (some) infinite metapopulations.

Growth rate and stable geographic distribution could have as well been computed respectively as the maximal eigenvalue and associated eigenvector of the mean offspring matrix (encompassing dispersal). This well-known spectral approach can still more efficiently be used in general for the numerical computation of these quantities than the techniques presented here. However, the solution presented here has the advantage of giving a clear biological interpretation of the contributions of demography and dispersal to growth and stable distribution. In particular, we hope that this new presentation will help researchers in conservation biology to make informed decisions in planning reintroduction strategies or in designing successfully protecting areas. Indeed, the benefit of our approach is to enable the ecologist to separate the effect on source-sink dynamics of changing the suitability of habitats (reproduction/survival scheme) or of changing the connectivity of the metapopulation (dispersal scheme).

Last, we want to indicate possible extensions of our ideas. First, we explain how migration-induced mortality has been encompassed in the growth phase, so that habitat patches i with m i > 1 are sources in the strict sense of the word. Indeed, let m ′ i denote the mean number of offspring produced in patch i before dispersal, and d ′ ij denote the probability of migrating from i to j and to survive to this migration event. Then

p i = K i=1 d ′ ij
is the probability of survival to migration starting in patch i. It is easily seen that m i = p i m ′ i is the mean offspring number including survival to migration and that

d ij = p -1 i d ′ ij
are the dispersal probabilities of surviving individuals, so that considering only individuals surviving migration, everything happens as if growth with mean m i preceded dispersal with weights d ij . These parameters m i and d ij are the ones we have used throughout the paper. Note that the matrix A ′ with generic element m ′ i d ′ ij is obviously equal to A. On the other hand, we could as well have sticked to the framework where mortality-induced migration is not encompassed in the growth phase, in which case the dispersal matrix D ′ is now sub-stochastic. Then the associated random disperser X ′ would die with probability 1 -p at each time step, and the criterion for persistence would have remained

m ′ 1 E T ′ -1 n=1 m ′ X ′ n > 1,
where T ′ is the first return time to patch 1 of the killed random walk X ′ , with the convention that the term inside the expectation is set to 0 when T ′ is infinite (death of the walker before returning home).

Second, we can relax the assumption that dispersal behaviors of siblings are independent, provided there is no correlation between dispersal behaviors of different groups of siblings. Indeed, in this case, we just need to consider the mean number m ij of offspring of an individual living in patch i which go into patch j. Setting

m i = K j=1 m ij , d ij = m ij /m i
allows us to come back to our framework. In the case of finite multitype branching processes, the questions handled here are only linked to the mean offspring matrix. It makes such a procedure valid. Third, if the model was expressed in continuous time then the disperser X would be a time-continuous random walk, and the criterion for persistence to generalize is the criterion involving occupation times (see Remark 1).

Fourth, we mention that our ideas could also be adapted to age-structured populations, by considering the stopping line of descendants of a focal juvenile ancestor who are the first descending juveniles to be born in the ancestor patch. Similarly, we can consider different (but a finite number) phenotypes or genotypes just by increasing the state space of the type of the branching process. Roughly speaking, habitat types are now replaced with a new composite type encompassing habitat, age, phenotype, genotype...

There is one further question that our methods could possibly solve. In infinite metapopulations, persistence can occur with the population failing to fill out the whole space, as in cases where dispersers always follow the same direction. This phenomenon is known as a dichotomy between local extinction and local exponential growth (conditional on global persistence). Our prediction of global persistence in infinite metapopulations relies on the study of the path of a random disperser started in a source patch, until the first time when it returns to a (possibly different) source patch. By making the difference between cases when the final patch is the same source patch or another source patch, we could display criteria for local persistence. Indeed, we know that there is global persistence in a source-transitive patch iff

m 1 E T -1 n=1 m Xn > 1,
where T is the first return time to a source patch, but we conjecture that in this case there is local persistence only if we also have

m 1 E τ -1 n=1 m Xn > 1,
where τ is the first return time to the very same source patch as initially.

A Proof of Theorem 4.1

We assume (A2', A3, A4, A5). Thus the Markov chain X is irreducible and evolves on a finite state space. The real number I(f ) is the cost for the habitat occupation frequencies associated with the random walk X to equal f . It gives the geometric decrease of the probability that the portion of time spent in habitat i until generation n is close to f i :

I(f 1 , . . . , f K ) = lim ǫ→0 lim n→∞ - 1 n log P(f k -ǫ ≤ F k (X n ) ≤ f k + ǫ). (4)
This result holds when I is finite. It is guaranteed by Sanov's theorem (see e.g. [9, Theorem 3.1.6 page 62]), which also ensures that I is convex continuous. The function I is called the rate function associated to the path of the random walk X.

Finally, by (A2'), the offspring distribution N i for an individual living in patch i satisfies E(N i log + N i ) < ∞, which ensures that for all i, j = 1, . . . , K, Expression of the growth rate ρ and habitat occupation frequencies. We start with one individual in patch 1. Recalling that for every i = 1, . . . , K, E(Z

E(Z j 1 log + Z (j) 1 | Z (i) 0 = 1, Z (k
(i) n ) = K j=1 E(Z (j)
n-1 )m j d ji , we get by induction

E(Z (i) n ) = j 0 =1, jn=i, 1≤j 1 ,...,j n-1 ≤K n-1 k=0 m j k d j k j k+1 = E 1 1 Xn=i n-1 i=0 m X i .
This yields

E(|Z n |) = E n-1 i=0 m X i
Denoting by S n (i) the number of visits of the random disperser X in patch i = 1, . . . V, :

S n (i) = #{k ≤ n -1 : (X k ) = i} = nF i (X n ),
we deduce

E(|Z n |) = E K i=1 m Sn(i) i = F exp(n K i=1 f i log(m i ))P(F 1 (X n ) ∈ df 1 , . . . , F K (X n ) ∈ df K ).
Using (4) and Laplace method, we get

log(E(|Z n |) 1/n ) n→∞ -→ max{ K i=1 f i log(m i ) -I(f 1 , . . . , f K ) : f ∈ F}.
This proves the first part of the result.

The maximum of h := R -I is reached for a unique frequency ϕ, which means that there is a unique ϕ such that log(ρ) = R(ϕ) -I(ϕ). This is the object of Proposition A.1. Moreover the partial derivatives of h at ϕ are zero. As

f K = 1 -f 1 -. . . f n-1 , for every 1 ≤ i ≤ K -1, log(m i ) -log(m K ) - ∂ ∂f i I |f =ϕ = 0.
If ϕ = p, then the partial derivatives of I at ϕ are zero. This can be directly computed or deduced from I ≥ 0 and I(p) = 0. This ensures that for every 1 ≤ i ≤ K -1, log(m i )log(m K ) = 0, i.e. m 1 = m 2 = . . . = m K . This proves the third part of the theorem. Finally, let us prove that the habitat occupation frequencies of a typical individual are given by the vector ϕ. Let ǫ > 0 and i ∈ {1, . . . , K}. The individuals alive in generation n are labelled by u k , k = 1, . . . , Z n .

Following the first part of the proof,

E(

Zn k=1

1 |F i (u k )-ϕ i |≥ǫ ) = E 1 |F i (Xn)-ϕ i |≥ǫ n i=0 m X i = E 1 |F i (Xn)-ϕ i |≥ǫ K i=1 m Sn(i) i = F 1 |f i -ϕ i |≥ǫ exp(n K i=1 f i log(m i ))P(F 1 (X n ) ∈ df 1 , . . . , F K (X n ) ∈ df K ).
Using again (4) and the Laplace method, we get

1 n log E Zn k=1 1 |F i (u k )-ϕ i |≥ǫ n→∞ -→ C i,ǫ , with C i,ǫ = max{ K i=1 f i log(m i ) -I(f 1 , . . . , f K ) : f ∈ F, |f i -ϕ i | ≥ ǫ}.
The uniqueness of the argmax ϕ ensures that C i,ǫ < C i,0 . Moreover the growth rate C i,0 is equal to log(ρ) and f j log(u j /(uD) j ) : u ∈ R K , u ≫ 0}, is reached for a unique unit positive vector. This means that there exists a unique u(f ) = (u 1 , . . . , u K ) such that

I(f ) = K j=1 f j log(u j (f )/(u(f )D) j ), u 1 (f ) + • • • + u K (f ) = 1, u 1 (f ) > 0, . . . , u K (f ) > 0.
Indeed this vector u(f ) realizes a maximum for u ∈ R K , u ≫ 0 and thus satisfies for j = 1, . . . , K, f j u j -

K i=1 d ji f i (uD) i = 0. ( 5
)
This equation characterizes u, see Exercise IV.9 page 46 in [START_REF] Hollander | Large deviations[END_REF], which ensures that u(f ) is uniquely defined. Note also that if f i is the stationary distribution, u j = f j satisfies this equation since (f D) i = f i , so that I = 0.

Proposition A.1 There exists a unique ϕ ∈ F such that log(ρ) = R(ϕ) -I(ϕ).

Proof. We observe that f → u(f ) can be extended from F to [0, ∞) K \ {0} and can satisfy Observe that the right hand side does not depend on i, so that for every i = 1, . . . , K,

(uD ′ ) i = αu i ,
where D ′ ji = d ji m i and α is a positive constant. Then u is left eigenvector of D ′ with positive entries. Moreover D ′ is strongly irreducible since D is strongly irreducible and m i > 0 for every i by assumption (recall that D is strongly irreducible if it is both irreducible and aperiodic, that is, if there exists n 0 ≥ 1 such that all the coefficients of D n 0 are positive). Now Perron-Frobenius theory ensures that there is a unique left positive eigenvector u of D ′ such that K i=1 u i = 1. Indeed, it is known [START_REF] Horn | Topics in matrix analysis[END_REF] that two positive eigenvectors of a primitive matrix are colinear. This actually comes from the classical decomposition of A n using the maximum eigenvalue and the associated left and right eigenvectors. Moreover, following the literature on large deviations [START_REF] Hollander | Large deviations[END_REF], (5) reads

f j = K i=1 D ′′ ji f i ,
with D ′′ ji = u j d ji /(uD) i . Here again D ′′ is strongly irreducible since D is strongly irreducible and both u and uD are positive vectors. Using again Perron-Frobenius theory guarantees the uniqueness of the solution f such that K i=1 f i = 1. This ensures the uniqueness of the argmax of R -I in the interior of [0, ∞) K . We complete the proof by adding that there is at least one argmax in the interior of [0, ∞) K since we recall that the frequency occupation is the product of the right and left eigenvectors of A associate to ρ, which are both positive (using again Perron Frobenius theory with assumptions (A3, A4, A5)). If ϕ 1 and ϕ 2 realize the max of R -I, the concavity of this function ensures that so do all elements in the segment [ϕ 1 , ϕ 2 ], which is in contradiction with the uniqueness in the interior of [0, ∞) K and completes the proof. 2

B Classification Theorem for MBPRE

We consider here a multitype branching process in random environment Z n = (Z (i)

n : i = 1, . . . , K) whose mean offspring matrix is denoted by A = A(w), where w = (w 0 , w 1 , . . .) is the environment. We introduce the extinction probability vector in environment w starting from one individual in habitat i: q i (w) = lim Proposition B.1 ([40], Theorems 9.6 and 9.10) Assuming that P(q i (w) < 1 : i = 1, . . . , K) = 1 or P(q i (w) = 1 : i = 1, . . . , K) = 1, ( * )

we have

• If γ < 0, then the probability of extinction is equal to 1 for almost every w.

• If γ = 1, then

Figure 1 :

 1 Figure 1: Two patches of different qualities. The filled circle is a source and the empty circle is a sink. The arrow labels are the dispersal probabilities.

Figure 2 :

 2 Figure 2: Three examples where sources are only connected to sinks and two adjacent sources are connected by n sinks; a) cyclic graph with one source, n = 7; b) two sources, n = 3; c) infinite periodic array, n = 2, arrow labels (not represented) are also assumed periodic.

Figure 5 :Figure 6 :

 56 Figure 5: A pipeline where n identical sinks (n = 7) connect the source to itself.

Figure 8 :

 8 Figure 8: Star sources with 2d pipelines of sinks a) (and periodic arrow labels -not represented, here d = 2) form a graph that can be collapsed into b) (a single source with d pipelines, loop edges are not represented). Notation Interpretation vertex P patch set V of vertices motif type of P (P ) = i habitat quality of patch P is i oriented edge with weight d P Q probability of dispersal from patch P to patch Q m i mean number of offspring in a patch with habitat type i

=P

  the graph into some motif V , we denote byZ (P ) n ′ ∈Cl(P ) N (P ′ )n the total number of individuals in some habitat of the class of patch P ∈ V . Then Z = (Z (P )

  and R(f) := P ∈V f P log(m (P ) ).

  ) 0 = 0 for k = i) < ∞.Using (A3, A4, A5), we have both A and D strongly irreducible. Then Z n /ρ n converges to a non degenerate variable W , see e.g. [4, Chapter 5, section 6, Theorem 1]. It is positive and finite on the survival event.

1111

  |F i (u k )-ϕ i |≥ǫ -1 n log ρ n n→∞ -→ C i,ǫ -C i,0 < 0. |F i (u k )-ϕ i |≥ǫ n→∞ -→ 0.In other words Zn k=1 1 |F i (u k )-ϕ i |≥ǫ /ρ n goes to 0 in probability. Adding thatZ n ∼ W ρ n a.s. as n → ∞ and {W > 0} = {∀n ∈ N, Z n > 0} a.s. ensures |F i (u k )-ϕ i |≥ǫ n→∞ |F i (u k )-ϕ i |≥ǫ n→∞ -→ 0.Then, conditionally on Z n > 0, denoting by U n an individual chosen uniformly in generation n,P (|F i (U n ) -ϕ i | ≥ ǫ, Z n > 0) n→∞ -→ 0.This proves thatF i (U n ) n→∞ -→ ϕ i inprobability and completes the proof. Study of I and uniqueness of argmax R -I. The supremum I defined by I(f ) = I(f 1 , . . . , f K ) := sup{ K j=1

( 5 ) 1 KK k=1 f k ] 2 .

 512 on [0, ∞) K \ {0} by settingu(f ) = u(f / f ), where f = K i=1 f i . Then R -I can be extended to [0, ∞) K \ {0} with (R -I)(f ) = log(m j (u(f )D) j /u j (f )), and (R -I)(λf ) = (R -I)(f ) for every λ ∈ (0, ∞).Consider a vector f which realizes the maximum of R -I and does not belong to the boundary of [0, ∞] K . Then the partial derivatives are zero and for every i = 1, . . . , K,k=1 f k . log(m i (uD) i /u i ) + log(m j (uD) j /u j ) [ i (uD) i /u i ) = K j=1 f j log(m j (uD) j /u j ) K k=1 f k .

n→∞P(

  |Z n | = 0 | w, Z (i) 0 = 1, Z j 0 = 0 if j = i).

Table 1 :

 1 Notation. and if ρ > 1, the metapopulation can survive with positive probability, in which case

Acknowledgement. We warmly thank Sebastian J. Schreiber and another (anonymous) reviewer for their feedback. This work was funded by project MANEGE 'Modèles Aléatoires en Écologie, Génétique et Évolution' 09-BLAN-0215 of ANR (French national research agency).

either for every m ≥ 1, w-a.s., there exists 1 ≤ i ≤ K such that

or q i (w) = 1 w-a.s..

Assuming that there exist integers N, L > 0 such that

• If γ > 0, then w-a.s q i (w) < 1 for every i = 1, • • • , K, and

Moreover thanks to [START_REF] Tanny | Limit theorems for branching processes in a random environment[END_REF]Theorem 9.11], if all the coefficients of the matrix A are positive and bounded, i.e.,


Introduction

We study in this paper the propagation of an optical beam inside the Output Mode Cleaner (OMC) cavity in order to test the tolerancing of construction specications. The main goal of this problem is to check if the beam stay inside polished mirrors during the propagation assuming random angular errors for mirror positions. The diameters of mirrors, assume to be 5mm, and the maximum allowed angular error, which we would determine, constitute the set of two relevant parameters in this problem if we exclude conguration parameters of the OMC. To solve the problem we will progress in two stpng, rst we will establish the equations of the beam propagation inside the cavity using the geometrical optics and second we will do the numerical resolution with a C++ code of those equations and nd the maximum allowed angular error.

Equations

2.1

Introduction to the mathematical problem

Conguration of the cavity

We can nd characteristics of the OMC in the Advanced Virgo Technical Design Report. We put a number for each mirror following the gure 1 and all variables in relation of a mirror will take the same subscript number. Angles are oriented in the counterclockwise. The axis going from the center O 1 to the center O 2 with the associated base vector e x will be the reference axis for all of this study. The origin of frames will be O 1 . We will use unitary vectors u i parallel to the mirror i (see Fig. 1 ). The pair of vectors (e x , u i ) make a base which is not orthogonal.

Finally, mirrors 1 and 4 are plane and mirrors 2 and 3 are spherical. We will consider beams only inside the cavity, input and output beams are deduced from corresponding inside beam with the refraction law.

2.1.2

Position of mirrors and beams

The position of mirrors and beams is done in the following way (see gure. 2):

• γ i is the angle between e x and u i which represent the angular position of the mirror i from the x axis.

• α i is the angle between an orthogonal axis to the mirror i and the incident beam.

• α ir is the angle between an orthogonal axis to the mirror i and the reected beam.

• L i is the distance between the center of the mirror i and the impact point of the beam. , -l sin(γ 3 )

), ( l tan(γ 4 )

, -l sin(γ 4 )

) . The angle α j of the incident beam on the mirror j is calculated as a function of the angle α ir by getting the angle θ. This angle represent the angular position of the beam from the x axis as shown on the gure 4. The angle θ is:

θ = γ i + π 2 + α ir = γ j + π 2 + α j (1)
Then we get:

α j = α ir + γ i -γ j (2)
Comment : In our case of closed cavity, with the equation 2 we nd the beam after reections on the four mirrors has to respect the following equation:

α 1 = -α 1r + 2 (γ 2 -γ 4 + γ 3 -γ 1 ) (3)
We observe it exist a solution for the beam to come back with the same angular position using plane mirrors only in the particular case of (γ 2γ 4 + γ 3γ 1 ) = 0

Beam position from center of mirrors

We will specify in this section the position of the incident beam on the mirror j from the center O j as a function of the reected beam characteristics on the mirror i. We will represent the propagation of the beam from i to j with the vector v ij . In order to obtain searched parameters we will make projection of all quantities on base vectors (e x , u j ).

First, we need to formulate u i in the chosen base as:

u i = a e x + b u j (4) VIR-0266B-12 -October 17, 2012
We know u i • e x = cos(γ i ), u j • e x = cos(γ j ) and u i • u j = cos(γ iγ j ), from that scalar products we obtain :

u i = sin(γ j -γ i ) sin(γ j ) e x + sin(γ i ) sin(γ j ) u j (5)
Second, we formulate in the same way the vector v ij as :

v ij = a e x + b u j (6)
We need following scalar products :

v ij • e x = -µ ij cos(γ i + π 2 + α ir ) = µ ij sin(γ i + α ir ) and v ij • u j = µ ij cos( π 2 -α j ) = µ ij sin(α j ).
We specify the direction of propagation by using the constant µ ij which is µ 12 = µ 43 = 1 and

µ 24 = µ 31 = -1
After few trigonometrical calculus, we obtain:

v ij = µ ij cos(α j ) sin(γ j ) e x -µ ij cos(α ir + γ i ) sin(γ j ) u j (7)
Now, we will express the coordinates (x O j , L j ) in the frame (e x , u j ) as a function of the coordinates (x O i , L i ). For this, we will write the equation of the line describing the propagation of the beam with the vector v ij , the origin (x O i , L i ) and the path parameter r:

(x O i e x + L i u i ) + rv ij (8)
The arriving impact point T on the mirror j correspond to the intersection between this line and the line representing the mirror j which has the coordinate x O j along the x axis. This intersection correspond to the following relation :

x O i e x + L i u i + r ij v ij = x O j e x + L j u j (9)
Since the vectors e x and u j form a base, the previous equation is veried when we have :

x O i + L i sin(γ j -γ i ) sin(γ j ) + r ij µ ij cos(α j ) sin(γ j ) = x O j (10)
With the equation 10 we get r ij :

r ij = µ ij x O j -x O i -L i sin(γ j -γ i ) sin(γ j ) sin(γ j ) cos(α j ) (11)
Then we get the other coordinate L j from r ij :

L j = L i sin(γ i ) sin(γ j ) -µ ij x O j -x O i -L i sin(γ j -γ i ) sin(γ j ) sin(γ j ) cos(α j ) µ ij cos(α ir + γ i ) sin(γ j ) (12)
VIR-0266B-12 -October 17, 2012

Curved mirrors 6/13

Let rewrite this expression as :

L j = L i sin(γ i ) sin(γ j ) -x O j -x O i -L i sin(γ j -γ i ) sin(γ j ) cos(α ir + γ i ) cos(α j ) (13)
With equations 2 and 13 we can compute all reexions inside the cavity taking care to put this additional equation for the reection law:

α ir = -α i ( 14
)
2.3 Curved mirrors

Characterization of spherical mirrors

At this point, we will take into account the curvature radius of mirrors 2 and 3. First, we need to determine the coordinate position (x K 2 , L K 2 ) and (x K 3 , L K 3 ) of respective sphere centers K 2 and K 3 in the corresponding vector basis. As shown on the gure 5 we get:

(x K 2 , L K 2 ) =   L - R cos(γ 2 - π 2 ) , -R tan(γ 2 - π 2 )   = L - R sin(γ 2 ) , R tan(γ 2 ) (15)
With the gure 6 we get:

(x K 3 , L K 3 ) = L + l tan(γ 3 ) - R sin(γ 3 ) , - l sin(γ 3 ) + R tan(γ 3 ) (16)
VIR-0266B-12 -October 17, 2012 We will determine the intersection points T ′ 2 and T ′ 3 of incident beams with circle associated to the curved mirror 2 and 3.

The equation of such circles is:

(x ′ O j -x K j )e x + (L ′ j -L K j )u j ) 2 = R 2 (17)
This give:

(x ′ O j -x K j ) 2 + (L ′ j -L K j ) 2 + 2(x ′ O j -x K j )((L ′ j -L K j )u j • e x = R 2 (18)
We get nally:

(x ′ O j -x K j ) 2 + (L ′ j -L K j ) 2 + 2(x ′ O j -x K j )(L ′ j -L K j ) cos(γ j ) = R 2 (19)
We formulate the line equation corresponding to the beam propagation using vector basis (e x , u j ) : The upper part show the dierent iteration of the program with the iteration numbers in the rst column, the angle α 1 in the second one, the position L 1 for the third, and dierences δL 1 and δα 1 for the two last column. The lower part give results for each mirrors respect to the commented column.

x O i e x + L i u i + r v ij = x O i e x + L i sin(γ j -γ i ) sin(γ j ) e x + sin(γ i ) sin(γ j ) u j + rµ ij cos(α j ) sin(γ j ) e x - cos(α ir + γ i ) sin(γ j ) u j (20)

Results

The main idea is to test the tolerancing of construction specications for the OMC, thus it is essential to take the worst case of angular error conguration. We can nd this condition by doing several tests and we get the following sign conguration for the angular error ε γ :

γ 1 = γ 10 -ε γ or + ε γ γ 2 = γ 20 + ε γ or -ε γ γ 3 = γ 30 + ε γ or -ε γ γ 4 = γ 40 -ε γ or + ε γ ( 31
)
where γ i0 = 90 ± 8.876 • correspond to the ideal conguration of OMC. The two possibility of sign congurations are equivalent and give identical results. The set of parameters taken in the Advanced VIRGO TDR are:

γ i0 = 90 ± 8.876 • L = 60 mm (32) l = 19.20 mm (33) RoC = 1499 mm (34)
The precision chosen for nding roots is a dierence of δL 1 and δα 1 less than 10 -6 . Within this conguration we can plot position oset L 1reso and tilt α 1reso in comparison to the ideal OMC.

VIR-0266B-12 -October 17, 2012 The rst observation of the curve in gure 8 is the oset on each mirror is quite close and strongly linear dependent on the angle error. Assuming the equation L 1reso = a ε γ the slope for the two lines are a = 52.3 mm • ( • ) -1 for the upper one and a = 51.2 mm • ( • ) -1 for the lower one. This linear behaviour is not so surprising since we are dealing with very small angles ε γ and this has the consequence to make this problem mainly dependent to the rst order of ε γ .

With ε γ = 0.03 • , which was the specication to be checked, we nd an oset L 1reso = 1.57 mm. Then, it is usual to take 2.5 r for calculate the edge of mirrors. Assuming the beam has a radius around r = 0.3 mm (the beam waist is 0.256 mm), this constraint corresponds to 1.57 + 2.5 r = 2.32 mm, which is less than the radius of the polished surface mirror (2.5 mm). These results conrm that an error of 0.03 • is acceptable for the actual conguration of OMC even in the worst case. Similar results had been obtain by Romain Gouaty with a naive model of cavity which conrm the coherence of such results with the code.

In the same idea we can check the angle oset (α iresoγ i0 ) as a function of angle error ε γ . We observe in gure 9 that the angular oset is also enough small to does not have any VIR-0266B-12 -October 17, 2012

Conclusion

The tolerancing ε γ = ±0.03 • conrm to us that the resonant beam does not exit from polished surfaces and does not have any impact on the optical set-up.

The code can be used for other congurations and can be easily modied to adapt it to other kind of cavity such as one spherical mirror for example.

The equations developed in the rst part can be used in other contexts and other programming languages.

VIR-0266B-12 -October 17, 2012

Figure 1 :

 1 Figure 1: Conguration of the cavity

Figure 2 :

 2 Figure 2: Variables of position for mirrors and beams

Figure 3 :VIR

 3 Figure 3: Positions of mirrors 3 and 4

Figure 4 :

 4 Figure 4: Beam positions

Figure 5 :

 5 Figure 5: Center position of the spherical mirror 2

Figure 6 :

 6 Figure 6: Center position of the spherical mirror 3

Figure 8 :

 8 Figure 8: Oset position L ireso as a function of angle error ε γ .

Figure 9 :

 9 Figure 9: Angle oset (α iresoγ i0 ) as a function of angle error ε γ .

VIR-0266B-12 -October 17, 2012

-8.906 0 -0.3776254233884 0.003930800733225

-8.9058 0.001 -0.37871115559 0.003932149038560

-9.173909486653 -0.347783936083 0.7576565219381 0.004552065445260 VIR-0266B-12 -October 17, 2012

Including equation 20 in the equation 19 we get the equation giving the parameter r ′ ij for the intersection points T ′ j :

We rewrite this equation using quantities a, b, c, d in order to simplify the equation:

Then we obtain:

The solution which we are seeking, is at the "right" side of the circle:

Reection on spherical mirrors

The radius joining the center K j with the point T ′ j represent the orthogonal to the tangent of the spherical mirror at this point. With the coordinates of these points we get unit vector u ′⊥ j associated to this orthogonal direction as following:

From the gure 7 we calculate the angular correction γ ′ j to make transition from plane mirror to spherical mirror.

VIR-0266B-12 -October 17, 2012 9/13 2 EQUATIONS

Figure 7: Reections on the spherical mirror j

We deduce from the gure 7 the equation for reection law including angular correction due to sphericity:

Variables in an equivalent system of mirror parallel to a reference plane mirror

The position along x axis of the equivalent mirror j = 2, 3 which pass at the point T ′ j is:

In the same way we nd the position of the incident beam on the equivalent mirror along the u j axis

VIR-0266B-12 -October 17, 2012 10/13

3 Numerical resolution

3.1

The method and the code

The cavity is assumed to be at the resonance when the propagating beam comes back with the same position and the same incident angle. For the ideal conguration of OMC, the position L 1reso and the angle α 1reso , which make resonant the beam inside the cavity, are obvious, i. e. L 1reso = 0 and α 1reso = -γ 1 . But in the case of random errors for the angular position of mirrors there are no obvious solutions. Since analytical solutions for the beam propagation inside the cavity with four mirrors are too complicated and not relevant for our problem, we will proceed by using an algorithm to nd the solution (L 1reso , α 1reso ). To help us, previous equations had been established in order to be used iteratively and therefore easily implementable in any codes. Initialization of the algorithm is done by taking the value (L 1 , α 1 ) assuming the position of mirrors 1 would be the ideal conguration for the beam, i. e. (L

where ε is the angular error. The procedure of the algorithm consist to nd roots of the two functions δL 1 (L 1 , α 1 ) and δα 1 (L 1 , α 1 ) which are the dierences of value L 1 and α 1 before and after a round trip of the beam inside the cavity. Then at each iteration we correct the couple of parameters (L 1 , α 1 ) to have null dierences δL 1 and δα 1 . The algorithm has been optimized by using the "good Broyden's method" for nding the two roots of two variables functions.

The program works with an input le of parameters and create an output le with results. Below is shown an example of input le, "par_cavite.dat", with a set of parameters , all of them are commented. For simplicity of input le the sign of angular error implemented for each mirror is hard written in the code but can be changed easily before compiling the code again. The compilation is done by writing in a command terminal " make version11" after going in the corresponding le, such as "version11_recherche_resonnance_optimisee_fonction_erreur_angulaire/". The code is started by typing in a command terminal " ./version11 ", then we obtain this output le "donnees_sortie_v11.dat":