Thierry E Huillet

ON A MARKOV CHAIN MODEL FOR POPULATION GROWTH SUBJECT TO RARE CATASTROPHIC EVENTS

Keywords: Population growth, Markov chain, catastrophic events, height and length of excursions, scaling PACS 87.23.Cc, 02.50.Ey MSC primary 60J10, secondary 42C05

We consider a Markov chain model for population growth subject to rare catastrophic events. In this model, the moves of the process are getting algebraically rare (as from x -λ ) when the process visits large heights x, and given a move occurs and the height is large, the chain grows by one unit with large probability or undergoes a rare catastrophic event with small complementary probability ∼ γ/x. We assume pure reflection at the origin. This chain is irreducible and aperiodic; it is always recurrent, either positive or null recurrent.

Estimates are obtained for first-return time probabilities to the origin (excursion length), eventual return (contact) probability and excursion height. All exhibit power-law decay in some range of the parameters (γ, λ). We show a scaling relationship between heights and lengths of the excursions. From this, the mean and median of both the empirical average and sample maximum are shown to grow algebraically with exponents being identified in terms of (γ, λ).

Introduction

We study a particular instance of a discrete-time Markov chain subject to rare catastrophic events or disasters. This model is in the class of the "house-of-cards" processes as the authors of ( [START_REF] Fernandez | Coupling, renewal and perfect simulation of chains of infinite order[END_REF], p. 47) and ( [START_REF] Ferrari | Construction of stochastic processes, coupling and regeneration[END_REF], p. 9) suggestively call them. In contrast to most studies on similar subjects in the literature (also often dealing with continuous times), we limit ourselves to the case where only birth or holding events can occur at times when disasters are ruled out. See [START_REF] Brockwell | Catastrophe processes with continuous statespace[END_REF], [START_REF] Cairns | Extinction times for a general birth, death and catastrophe process[END_REF], [START_REF] Gripenberg | A stationary distribution for the growth of a population subject to random catastrophes[END_REF] and [START_REF] Stirzaker | Processes with catastrophes[END_REF].

Our model can precisely be described as follows.

Let λ ≥ 0, γ > 0. Consider the following Markov chain X on Z + = {0, 1, ...} : Given X n = x ∈ {1, 2, ...}, the increment of X n is +1 with probability : p x = x -λ e -γ/x 0 with probability : r x = 1 -x -λ -x with probability : q x = x -λ 1 -e -γ/x ∼ γ/x λ+1 . λ > 0 accounts for the fact that moves of X get rare when x gets large, and given a move occurs and x is large, X grows by one unit with large probability e -γ/x or 1 undergoes a rare catastrophic event with small complementary probability ∼ γ/x. We assume pure reflection at the origin. This chain is irreducible and aperiodic; it is always recurrent, either positive or null recurrent.

The birth and holding probabilities depend on the current state of the population in the specific way just described. A model with a similar flavor can be found in [START_REF] Karlin | Occupation time laws for birth and death processes[END_REF] although the Authors deal here with continuous-time birth and death processes, which is of course not our case. The study of stochastic stability of the model with long-range-dependence and power-law crash sizes developed in [START_REF] Robert | Stochastic stability of some state-dependent growth-collapse processes[END_REF] seems to be the closest to our own purposes.

Our study includes first-return time probabilities to the origin (excursion length), eventual return (contact) probability to the origin, excursion height, time to failure and estimates of the growth of the mean and median of both the empirical average and sample maximum. It is based on a scaling relationship between excursions heights and lengths. Because catastrophic events are assumed to be rare, the population growth is shown to be slow and in fact algebraic in the number of steps. We identify the algebraic exponents of growth in terms of (γ, λ) .

Motivations and relation to similar works in the literature

Before we proceed with the mathematical study of our specific model, let us briefly give some motivations together with the mentioning of some relations to the vast existing literature on growth-collapse models. Growth-collapse processes where long periods of growth are interrupted by rare crash events occur in a large variety of systems. For instance in applications related to Transmission Control Protocol (see [START_REF] Robert | Stochastic stability of some state-dependent growth-collapse processes[END_REF] for a precise description) or, as we will explain it in more details in the next Section, in the context of machine replacement (see [START_REF] Martin | Batch queues, reversibility and first-passage percolation[END_REF]) in Queuing Systems. But also in the Physics of Self-Organized-Criticality such as sandpile growth processes where a long accumulation of grains is brutally interrupted by the occurrence of a critical avalanche (see [START_REF] Bak | How nature works. The science of self-organized criticality[END_REF]). Self-Organized-Criticality where an energy accumulation occurs preceding a rapid energy redistribution (or stress-release) is also suspected to explain power-law distributions of earthquakes sizes (see [START_REF] Bak | Unified scaling law for earthquakes[END_REF], [START_REF] Carlson | Dynamics of earthquake faults[END_REF] and [START_REF] Corral | Long-term clustering, scaling and universality in the temporal occurrence of earthquakes[END_REF]). And also, as we do here typically, growth-collapse models occur as models for population growth subject to rare catastrophic extinction events. Not to briefly mention crashes in financial mathematics.

As we mentioned in the Introduction, most works deal with continuous times, see for instance [START_REF] Boxma | A Markovian growth-collapse model[END_REF]. In several studies, the growing part of the process is more complicate than ours; for instance in [START_REF] Eliazar | A growth-collapse model: Lévy inflow, geometric crashes, and generalized Ornstein-Uhlenbeck dynamics[END_REF] and [START_REF] Kella | On growth-collapse processes with stationary structure and their shot-noise counterparts[END_REF], the Authors consider it as a one-sided Lévy subordinator with times and sizes of the crashes depending on the system's state. In [START_REF] Eliazar | Stochastic Ornstein-Uhlenbeck capacitors[END_REF], similar ideas are applied to model the charge and discharge of a capacitor. In [START_REF] Eliazar | Growth-collapse and decay-surge evolutions, and geometric Langevin equations[END_REF], the inflow process undergoes a smooth deterministic growth/decay evolution with perturbations designed in such a way that the overall process obeys a geometric Langevin equation driven by a state-dependent noise. On the other hand in [START_REF] Boxma | A Markovian growth-collapse model[END_REF], the inflow process is constant but still with state-dependent times and magnitudes of the crashes. In all these papers, the main statistical characteristics of the systems at stationarity are computed.

To a large extent, the Markov chain model under study here is the simplest possible in the vast family of growth-collapse models. Using ideas stemming from excursion theory, we make precise here to what extent growth is (algebraically) slow when catastrophic events are rare under our over-simplistic model hypothesis. It is hoped that models in the same class of universality could share similar behaviors as one reasonably can expect some sort of robustness.

3. The Model: A special Markov chain with catastrophe on Z + Before proceeding with the MC model with holding probabilities, let us first investigate the simpler case where holding probabilities are absent. Let γ > 0. We shall first consider the following discrete-time homogeneous Markov chain X := (X n ; n ≥ 0) with state-space Z + = {0, 1, ...} and transition probabilities characterized by:

• given X n = x ∈ {1, 2, ...}, the increment of X n is +1 with probability : p x = e -γ/x -x with probability : q x = 1 -e -γ/x ∼ γ/x.

• given X n = 0, the increment of X n is +1 with probability p 0 = 1.

In other words, with (U n , n ≥ 1) a sequence of independent identically distributed (i.i.d.) uniform random variables (rvs), the dynamics of X n reads

X n+1 = (X n + 1) 1 (U n+1 > q Xn ) .
This irreducible Markov chain (MC) is in the class of general MCs with catastrophes (whose transition probabilities are state-dependent). It is reflected at the origin. In this model, the walker X n is occasionally bounced back to the origin. When x is large, the drift of this MC is of order 1 -γ -γ/x. So when γ > 1, the walker is attracted to the origin: The strength of the attraction goes like γ -1 for large x. For γ < 1, the walker is repelled from the origin correspondingly. When γ = 1, its drift is still attracting but of order -1/x, and the drift that the walker feels vanishes when x approaches ∞. We will see that while crossing the critical value γ = 1 from above, the process X switches from positive recurrent to null-recurrent.

In the context of machine replacement in reliability theory, one may classically interpret this MC as follows: at time n = 0, a machine is put into service. This machine has a random lifetime, say τ 0,0 . After τ 0,0 , a new machine (with lifetime a copy of τ 0,0 ) is installed to replace the old defective one and so on. The failure epochs of the successive machines constitute a renewal process on N 0 generated by τ 0,0 -1. In this context, X n clearly represents the age of the machine currently in action (the time till the last machine failed before it was replaced by the new one in action at time n). Given the age of the current machine is x, there is an age-dependent probability p x that the machine will survive one more time unit and a probability q x that it will not. Clearly, the law of τ 0,0 is P (τ 0,0 = k + 1) = k-1 y=0 p y q k , k ≥ 1, in terms of the ps and the τ 0,0 s are seen to be the times between consecutive visits to 0 for X. While assuming, as we do, that p 0 = 1, we suppose that the new machine is put into service as soon as the old one fails (there is no latency waiting time).

One may also interpret this MC as follows: X n represents the size of some growing population at time n. Given X n = x, at the next time unit, a new element presents itself for possible integration to the clan. The population can then grow by one unit with size-dependent probability p x , integrating normally the new element, so: X n → X n+1 = X n + 1 with probability p x . But there is a 'chance' q x that the new element is a black sheep (a "terminator") at the contact of which the whole population will get decimated, so: X n → X n+1 = 0 with probability q x ; the black sheep is then the only survivor in the next generation. Clearly the occurrence of a black sheep is a catastrophic renewal event.

The contamination (or collapse) probability q x may be a decreasing or an increasing function of the current size x. In the former case, large populations are getting more and more immune to black sheep and one expects that large population sizes will get stable. This is the case under study here with q x = 1 -e -γ/x ∼ x→∞ γx -1 . We will see that such a MC X is always recurrent. In case for instance

q x ∼ x→∞ γx -α (α, γ > 0), we would check that X is transient if α > 1, positive recurrent if α < 1, whereas if α = 1, X is positive recurrent if γ > 1, null recurrent if γ ≤ 1.
In the latter case, large populations are more susceptible and vulnerable to black sheep and so quite unlikely to grow large and develop too much. Think of a forest fire occasionally destroying completely an otherwise regularly growing forest. This would be the case for a model with q x = 1 -e -γx . If one also thinks of the process of building a house of cards, clearly adding a new card to a house of cards of size x is more likely to lead to a collapse of the whole structure if x is already large.

Including holding probabilities. An important ingredient of our model is the opportunity to include holding probabilities, making the chain's changes algebraically rare. Let λ ≥ 0. We shall in fact study the following more general Markov chain X, which may be viewed as a time-changed of X if λ > 0 :

• given X n = x ∈ {1, 2, ...}, the increment of X n is +1 with probability : p x = x -λ e -γ/x

(1) 0 with probability : r x = 1 -x -λ -x with probability :

q x = x -λ 1 -e -γ/x ∼ γ/x λ+1 .
• given X n = 0, the increment of X n still is +1 with probability p 0 = 1.

Clearly now, the dynamics of X n reads

X n+1 = X n 1 (U n+1 ≤ r Xn ) + X n + 1 1 r Xn < U n+1 ≤ p Xn + r Xn .
When x is large, the holding probability r x is close to 1 if λ > 0 so the chain X remains steady with large probability. Whenever X moves, with large probability e -γ/x , it moves up by one unit and only exceptionally, with probability 1 -e -γ/x ∼ γ/x, does the catastrophic event take place. In this context, X n is just the size of the cluster since the last catastrophe. This model therefore accounts for the lazy growth of stable large populations when the events leading either to growth or to collapse are themselves getting rare.

First properties of the general catastrophe Markov chain

In this Section, we first supply a detailed study of the general catastrophe Markov chain, without specifying the catastrophe probabilities q x . 4.1. Without holding probabilities. Consider a general catastrophe process X n reflected at the origin (p 0 = 1 and q 0 = 0) for which both p x and q x > 0, for all x ≥ 1, with p x + q x = 1 and so with associated stochastic transition matrix: P = [P (x, y)], (x, y) ∈ Z 2 + with P (x, 0) = q x and P (x, x + 1) = p x , x ≥ 1.

Existence and shape of the invariant measure. Let π ≡ (π 0 , π 1 , ..) be the row-vector of the invariant measure, whenever it exists. Then π should solve π = π P, whose formal solution is:

(2)

π 0 = x≥1 π x q x and π x = π 0 x-1 y=0 p y , x ≥ 1. Let u x =
x-1 y=0 p y . Using the second equation, the first equation is satisfied whenever

x≥1 q x x-1 y=0 p y = x≥1 (u x -u x+1 ) = 1, so also when u ∞ = ∞ y=1 p y = 0 which is fulfilled iff S 2 ≡ ∞ y=1 q y = ∞.
We first conclude that there exists an invariant measure iff S 2 = ∞.

If in addition, S 1 ≡ x≥1

x-1 y=0 p y < ∞, then π 0 = 1 1+S1 ∈ (0, 1) and the invariant measure is unique and is a proper invariant probability measure. In this case, with the empty product being 1, we have

π x = x-1 y=0 p y 1 + S 1 , x ≥ 0.
When S 1 = ∞, the measure solution to (2) exists but it is not a probability measure as its total mass π 0 (1 + S 1 ) sums to infinity.

Return time to the origin. Let X n be the Markov chain with transition probability matrix P and state-space N 0 . Starting from X 0 = x ≥ 1, the walker moves one step up with probability p x or goes back to the origin (the catastrophic event) with probability q x and then the process starts afresh from 0. Clearly of interest are the times τ 0,0 between consecutive visits to 0 (the first return times to 0) because if X 0 = 0, X n represents the backward recurrence time (the time separating n to the previous visit to 0) of a discrete renewal process generated by τ 0,0 ≥ 2. We have

P (τ 0,0 = k + 1) = k-1 y=0 p y q k = u k -u k+1 , k ≥ 1,
which is also

P (τ 0,0 > k) = k-1 y=0 p y = u k .
Note that p x = P (τ 0,0 > x + 1) /P (τ 0,0 > x) and q x = P (τ 0,0 = x + 1) /P (τ 0,0 > x) : If the law of the lifetime τ 0,0 is known in the first place, this gives the survival probability p x given the age (backward recurrence time) X n of the current machine is x.

Note also that P (τ 0,0 < ∞) iff u ∞ = 0 (S 2 = ∞) which is the recurrence condition for X n .

We conclude that:

If S 2 ≡ ∞ y=1 q y < ∞, X n is transient. If S 2 = ∞, X n is recurrent. If S 2 = ∞ and S 1 ≡ x≥1 x-1 y=0 p y < ∞, X n is positive recurrent with µ := E (τ 0,0 ) = 1/π 0 = 1 + S 1 < ∞. If S 2 = S 1 = ∞, X n is null recurrent with τ 0,0 < ∞ almost surely (a.s.
) and E (τ 0,0 ) = ∞. In the positive recurrent case,

π x = x-1 y=0 p y 1 + S 1 = P (τ 0,0 > x) µ , x ≥ 0
which is the law of the limiting (n → ∞) backward recurrence time X ∞ of the discrete renewal process generated by τ 0,0 -1 ≥ 1. Let indeed τ n be the length of the interval to which n belongs. Then

τ n d → τ (convergence in distribution) as n → ∞ with P (τ = k) = kP (τ 0,0 -1 = k) µ , k ≥ 1 which is the size-biased version of the law of τ 0,0 -1. If Y ∞ is the limiting (n → ∞) forward recurrence time of the discrete renewal process generated by τ 0,0 -1 ≥ 1, with X ∞ + Y ∞ = τ , we have X ∞ d = Y ∞ (distributional equality) and (see e.g. [20], Lemma 9, p. 447) (X ∞ , Y ∞ ) d = (U • τ , (1 -U ) • τ ) .
Here, with U a uniform random variable independent of τ , U • τ is the U -thinning of τ : U • τ = τ l=1 B l (U ) where, given U, B l (U ) ; l ≥ 1 are mutually independent and independent of τ rvs with law Bernoulli(U ) . Indeed

P (U • τ = x) = 1 0 P (u • τ = x) du and P (u • τ = x) = k≥x k x u k (1 -u) k-x P (τ = k) .
Thus,

P (U • τ = x) = k≥x k x 1 0 u x (1 -u) k-x duP (τ = k) = k≥x 1 k P (τ = k) = 1 µ P (τ 0,0 > x) = π x .
So there is a clear connection of our model with the discrete theory of renewal processes.

Time-reversal. It is of importance to check whether or not detailed balance holds for the MC under study here. Assume X n is recurrent. The catastrophe MC is not time-reversible as detailed balance does not hold. Let Q = P be the transition matrix of the process Y n which is X n backward in time. With denoting matrix transposition and D π =diag(π 0 , π 1 , ...), we have

Q = D π P D -1 π .
It can be checked that the only non-null entries of Q are its first row with Q 0,0 = Q 0,1 = 0 and Q 0,k = P (τ 0,0 = k) if k ≥ 2 and the lower diagonal whose entries are ones. Starting from Y 0 = y, the process Y n decays linearly till it hits 0 and once at state 0, Y n jumps abruptly upward, undergoing a jump of amplitude k ≥ 2 with probability P (τ 0,0 = k) before diminishing again and again to 0. Clearly, in case of positive recurrence, Y n models the forward recurrence time of the original process with

Y n d → Y ∞ as n → ∞. Note that π = π Q (π is also the invariant measure for Q) and so: X ∞ d = Y ∞ ∼ π x , as required.
The scale (or harmonic) function.

In the recurrent case, the sample paths of X n are made of i.i.d. excursions (the pieces of the sample paths between consecutive visits to 0). The lengths of the excursions are τ 0,0 . Let us look at their heights H. Clearly H d = τ 0,0 -1, because X grows linearly between consecutive visits to 0. Let us rapidly check this with the use of the scale function. The scale function idea will appear more useful when we shall deal next with the time-changed version of the catastrophe MC because here the relationship between height and length of an excursion is far from trivial. Assume X 0 = x. Let X n∧τ x,0 stopping X n when it first hits 0. Let us define the scale (or harmonic) function ϕ of X n as the function which makes Y n ≡ ϕ X n∧τ x,0 a martingale. The function ϕ is important because, as is well-known, for all 0 < x < x * , with τ x the first hitting time of {0, x * } starting from x

P (X τ x = x * ) = P (τ x,x * < τ x,0 ) = ϕ (x) ϕ (x * )
.

Using this remark, the event H = h is realized when τ 1,h < τ 1,0 and τ h,h+1 > τ h,0 , the latter two events being independent. Thus:

(3)

P (H = h) = ϕ (1) ϕ (h) 1 - ϕ (h) ϕ (h + 1) , h ≥ 1.
We clearly have h≥1 P (H = h) = 1 because partial sums are part of a telescoping series. Clearly, this is also P (H ≥ h) = 1/ϕ (h). It remains to compute ϕ with ϕ (0) = 0. We wish to have:

E x (Y n+1 | Y n = y) = y, leading to ϕ (x) = p x ϕ (x + 1) + q x ϕ (0) = p x ϕ (x + 1) , x ≥ 1.
Thus, the searched 'harmonic' function is

(4) ϕ (x) = 1 x-1 y=0 p y , x ≥ 1, ϕ (0) ≡ 0.
Note ϕ (1) = 1. Equations ( 3) and ( 4) characterize the law of the excursion height of the random walker. Note

P (H ≥ h) = 1/ϕ (h) = h-1 y=0 p y = P (τ 0,0 > h) .
showing, as expected from the beginning, that

H d = τ 0,0 -1.
Probability of extinction. Consider now the same Markov chain but assume now that p 0 = 0, q 0 = 1. In this case, the state 0 is absorbing. Consider then the restriction P of matrix P to the states {1, 2, ...}. Let φ x , x ≥ 1 be the probabilities that state 0 is hit in finite time given the chain started originally at x. Let φ ≡ (φ 1 , φ 2 , ..) be the column-vector of these absorption probabilities. Let q ≡ (q 1 , q 2 , ..) . Then φ is the smallest non-negative solution to φ = q + P φ whose formal solution is: φ = I -P -1 q. All φ x can therefore be expressed in terms of φ 1 , leading simply to:

1 -φ x = 1 x-1 y=1 p y (1 -φ 1 )
.

If S 2 ≡ y≥1 q y = ∞, then u x =
x-1 y=1 p y → 0 : the restriction φ x ∈ [0, 1] forces φ 1 = 1 and so φ x = 1 for all x < ∞: The state 0 is hit with probability 1, starting from x, for all x < ∞. The Markov chain is recurrent.

But, if S 2 < ∞, then we can take φ 1 < 1 so long as φ x ≥ 0 for all x ≥ 1. The minimal solution occurs when 1 -φ 1 = P (τ 0,0 = ∞) = y≥1 p y > 0, leading to:

(5)

φ x = 1 - y≥1 p y x-1 y=1 p y = 1 - y≥x p y .
In this case, φ x < 1 for x ≥ 1 and the absorbed random walker started at x avoids 0 with positive probability (a transience case for the original reflected Markov chain). Note that x < x ⇒ φ x > φ x .

To summarize, we have:

Proposition 1. (i) If S 2 < ∞, the MC is transient and, with τ x,0 = inf (n ≥ 1 : X n = 0 | X 0 = x) , P (τ x,0 = ∞) = y≥x p y > 0. (ii) If S 2 = ∞, the MC is recurrent with P (τ x,0 = ∞) = 0. Moreover, it is: null recurrent if S 1 = ∞, positive recurrent if S 1 < ∞.
Due to irreducibility (because p x and q x > 0, for all x ≥ 1), states are either all transient or recurrent.

Times to collapse. How long does it take, starting from x ≥ 1, to first hit 0?

We give here some insight on the way to compute the law of these first times to collapse. With x ≥ 1, let thus τ x,0 be the time it takes to first hit 0, starting from X 0 = x ≥ 1. With τ x+1,0 a statistical copy of τ x+1,0 , from first-step analysis, we clearly have:

τ x,0 d = (1 -B x ) • 1 + B x • 1 + τ x+1,0 ,
where B x is a Bernoulli random variable with P (B x = 1) = p x . Therefore with φ x (z) = E (z τ x,0 ) , φ x (z) obeys the recurrence φ x (z) = q x z + p x zφ x+1 (z), with initial condition φ 1 (z): again, φ x (z) can easily be deduced once φ 1 (z) is known. With φ 0 (z) = E (z τ 0,0 ) the probability generating function of the first return time to 0, τ 1,0

d = τ 0,0 -1 entails φ 1 (z) = φ 0 (z) /z. Because P (τ 0,0 = k + 1) = k-1 y=0 p y q k = u k -u k+1 , with u (z) = k≥1 z k u k , we have φ 1 (z) = k≥1 z k k-1 y=0 p y q k = k≥1 z k (u k -u k+1 ) = z -1 z u (z) + 1,
fixing the initial condition. The full probability generating function (p.g.f.) of τ x,0 follows.

With φ (z) = (φ 1 (z) , φ 2 (z) , ...) the column-vector of the φ x (z), and q = (q 1 , q 2 , ...) the column-vector of the q x , φ (z) solves:

(6) φ (z) = zq + zP φ (z) ,
whose formal solution is φ (z) = z I -zP -1 q, involving the resolvent of P . When z = 1, φ x (1) = φ x = P (τ x,0 < ∞) are the absorption probabilities already computed.

Transience versus recurrence. We here discuss the criterion for recurrence or transience of the general catastrophe Markov chain. When S 2 = ∞, the recurrent chain started at x first hits 0 with probability 1 and returns infinitely often to 0. Given X 0 = x, with N x,y ≡ n≥0 1 (X n = y) , the number of visits to state y, then N x,y = ∞, P x -almost surely. If τ x,x is the first return time at x, then

P (τ x,x < ∞) = 1. Furthermore, with N x,y ≡ τ x,x
n=0 1 (X n = y) the number of visits to state y before the first return time to state x, then: E (N x,y ) = πy πx and by the Chacon-Ornstein limit ratio ergodic theorem:

N n=0 1 (X n = y) N n=0 1 (X n = x) → N ∞ π y π x , P x -almost surely.
Starting in particular from x = 0, a recurrent chain is made of infinitely many independent and identically distributed (i.i.d.) excursions which are the sample paths of (X n ; n ≥ 0) between consecutive visits to state 0. We have:

E (N 0,x ) = πx π0 =
x-1 y=0 p y . When the chain is positive recurrent (S 1 < ∞) the expected time elapsed between consecutive visits to 0 is finite and equal to E (τ 0,0 ) ≡ µ = 1/π 0 = 1 + S 1 , whereas this expected time is infinite when the chain is null recurrent.

When S 2 = ∞, the state x ≥ 0 is transient. Thus, N x,x < ∞, P x -almost surely and

P (N x,x = k) = (1 -α x ) α k-1 x where α x = P (τ x,x < ∞) < 1.

4.2.

Including holding probabilities and time change. We now proceed with including the opportunity of a holding probability in the transition matrix of the catastrophe process.

Let ρ x ∈ (0, 1), x ≥ 1. Consider now a general catastrophe MC X n still reflected at the origin ( p 0 = 1 and q 0 = 0) and for which p x = p x ρ x , q x = q x ρ x and r x = 1 -ρ x > 0, for all x ≥ 1, with p x + q x + r x = 1. So X n is generated by the stochastic transition matrix: P = P (x, y) , (x, y) ∈ Z 2 + with P (x, 0) = q x , P (x, x) = r x and P (x, x + 1) = p x , x ≥ 1. There is now a positive holding probability r x that given X n = x ≥ 1, X n+1 = x. The Markov chain X is lazy.

In the context of machine replacement, the times at which X n = X n+1 = x are idle times, whereas in the population growth image, these times correspond to periods where no new candidate present themselves for possible integration to the current cluster (so in particular with no risk of appearance of the black sheep): the collapse is delayed. X n accounts for the population size since the last catastrophe.

Invariant measure. We now investigate the way the invariant measure is modified by the adjunction of holding probabilities. Let D ρ =diag(ρ 0 := 1, ρ 1 , ρ 2 , ...) . We have P = I + D ρ (P -I) . Let π be the invariant measure associated to P , when it exists. It should solve π = π P , and we get:

π 0 = x≥1 π x q x ρ x and π x = π 0 x-1 y=0 p y ρ x , x ≥ 1.
Using the second equation, the first equation is satisfied whenever

x≥1 q x ρ x x-1 y=0 p y ρ x = x≥1 (u x -u x+1 ) = 1, so again when u x = x-1 y=1 p y → x→∞ 0 which is fulfilled iff S 2 ≡
∞ y=1 q y = ∞ (the recurrence condition for X): The time change leading from X to X does not change the road map of X so the recurrence criteria are identical for both X and X. However, the criteria for positive recurrence are modified. Indeed, if in addition S 1 ≡ x≥1 ρ -1 x x-1 y=0 p y < ∞, then π 0 = 1 1+ S1 ∈ (0, 1) and then the MC is positive recurrent with invariant probability measure

π x = π 0 x-1 y=0 p y ρ x , x ≥ 0.
Else if S 1 = ∞, the MC is null recurrent.

Scale function. X n is again the size at time n of the current population before the last catastrophic event. In the recurrent case (S 2 = ∞), X is again made of i.i.d. excursions, whose height H has the same law as the one for X : indeed, one can easily check that the scale function ϕ of X solving P ϕ = ϕ, ϕ 0 = 0, coincide with the scale function ϕ of X solving P ϕ = ϕ, ϕ 0 = 0. Because the scale function determines the height's law, we have the claimed statement that the height's law is left unchanged by time substitution.

Excursions lengths. However, because of the time change, the times τ 0,0 between consecutive visits to 0 (the excursions lengths of X) are of course very different from τ 0,0 (statistically longer). With X 0 = 0, both times are related by

τ 0,0 = τ 0,0 n=1 1 X n = X n-1 .
We thus have

P (τ 0,0 = k) = l≥k P X l-1 = k, τ 0,0 = l ,
where, on the event τ 0,0 = l, X l-1 is the height H l of an excursion of length l for X.

Given the height of an excursion is h, we have

P τ 0,0 = l | H = h = q h P ( τ 0,h = l -1 | τ 0,h < τ 0,0 ) = q h P 0 X l-1 = h | τ 1,h < τ 1,0 .
The p.g.f. of τ 0,0 can be obtained as follows. With x ≥ 1, let τ x,0 be the time it takes for X to first hit 0, starting from X 0 = x ≥ 1. Therefore with φ x (z) = E z τ x,0 , φ x (z) obeys the recurrence φ x (z) = q x z + r x z φ x (z) + p x z φ x+1 (z), with initial condition φ 1 (z): again, φ x (z) can easily be deduced once φ 1 (z) is known. With φ (z) = φ 1 (z) , φ 2 (z) , ... the column-vector of the φ x (z), and q = ( q 1 , q 2 , ...)

the column-vector of the q x , φ (z) solves:

φ (z) = z q + z P φ (z) , whose formal solution is φ (z) = z I -z P -1
q, involving the resolvent of P . In particular, with e 1 = (1, 0, 0, ...), the exact p.g.f. of τ 0,0 is:

φ 0 (z) = E z τ 0,0 = z φ 1 (z) = z 2 e 1 I -z P -1
q.

Note that P is an upper-diagonal matrix with upper-diagonal entries p x and diagonal entries r x . So this matrix is diagonable with eigenvalues r x = 1 -ρ x . From this fact, in principle

φ 0 (z) = x≥1 A x z 2 (1 -r x ) 1 -z r x ,
showing that τ 0,0 is an infinite mixture of (shifted) geometrically distributed rvs with success probabilities r x .

Time reversal. We now show how to compute the transition matrix of the timereversed process corresponding to the catastrophe MC with holding probabilities. When X n is recurrent, this catastrophe MC is not time-reversible either. Let Q = P be the transition matrix of the process Y n which is X n backward in time. Y now models a population that shrinks slowly and continuously to 0 before undergoing an instantaneous rebirth allowing the process to start afresh. Thus Y n represents the size of the disintegrating cluster at time n till the next forthcoming rebirth event.

With denoting matrix transposition and D π =diag( π 0 , π 1 , ...), we indeed have

Q = D π P D -1 π . The only non-null entries of Q are its first row with Q 0,k = Q 0,k if k ≥ 1 ( Q 0,0 = 0), the lower diagonal whose entries are ρ k , k ≥ 1 and the diagonal terms Q k,k = 1-ρ k , k ≥ 1.
Starting from Y 0 = y, the process Y n decays linearly but with steady periods, till it hits 0 and once at state 0, Y n jumps abruptly upward, undergoing a jump of amplitude k ≥ 1 with probability Q 0,k = Q 0,k before shrinking again to 0.

The special catastrophe Markov chain

We now specify the results of the latter Section to the special MC defined in [START_REF] Bak | How nature works. The science of self-organized criticality[END_REF]. The special MC deserves special interest in particular because it is, to a large extent, amenable to exact analytic computations. We shall limit ourselves to the study of X with the two parameters γ > 0, λ ≥ 0 because the study of X (without holding probabilities) can be obtained from the one of the time-changed X while plugging λ = 0.

Consider then the special MC X, so determined by: p x = x -λ e -γ/x and q x = x -λ 1 -e -γ/x , x ≥ 1 and satisfying p 0 = 1 and q 0 = 0 (reflection at the origin). Given X n = x, the local drift at x of this MC is f (x) = p x -x q x . For all 0 < γ = 1, λ > 0, this drift tends to 0 as x tends to ∞ as a result of f (x) ∼ -(γ -1) x -λ . The asymptotic drift is < 0 (respectively > 0) when γ > 1 (respectively γ < 1). In the former (latter) case, X n is asymptotically attracted (repelled) by the wall. When γ = 1, f (x) ∼ -x -(λ+1) .

Its local variance at

x is g 2 (x) = p x + x 2 q x -( p x -x q x ) 2 , with g 2 (x) ∼ γx 1-λ . If λ < 1 (respectively λ > 1)
, the volatility at large x diverges (vanishes) whereas

g 2 (x) ∼ γ if λ = 1.
Both MCs X and X are always recurrent because, with q y ∼ γ/y, S 2 = y q y = ∞.

There is no transience case for the special MC. The invariant measure of X satisfies

π x ∼ π 0 x -(γ-λ) ,
suggesting that π ≡ ( π x ; x ∈ Z + ) is a proper (summable) probability measure if and only if γ -λ > 1. Using the results of the latter Section, we thus get Proposition 2. When γ -λ ≤ 1, the MC X is null recurrent.

When γ -λ > 1, the MC X is positive recurrent; it has a unique and summable invariant measure satisfying π x ∼

x ∞

x -(γ-λ) . It is given explicitly by [START_REF] Corral | Long-term clustering, scaling and universality in the temporal occurrence of earthquakes[END_REF] π 0 = 1

1 + S 1 = 1 1 + x≥1 x λ x-1 y=0 e -γ/y , and (8)
π x = π 0 x-1 y=0 p y ρ x = π 0 x λ x-1 y=0 e -γ/y , x ≥ 1. When γ -λ > 2, (9) E X ∞ = x≥1 x π x < ∞,
so that the invariant measure has a finite explicit mean.

Note that when γ -λ > 1 (positive recurrence), the invariant probability measure is heavy-tailed with power-law exponent γ -λ > 1. So X ∞ can get large which is a hint that large population sizes are stable as a result of q x ∼ γ/x, decreasing with x.

Excursion heights. The scale or 'harmonic' function of X is

(10) ϕ (x) = 1 x-1 y=0 p y ∼ x γ , x ≥ 1, ϕ (0) ≡ 0.
Therefore, the excursion heights of X have power-law tails

P H ≥ h = 1/ϕ (h) ∼ h -γ .
Excursion lengths. In order to estimate the lengths of a big excursion, we suggest the following heuristics: Consider an excursion whose length k and height h are both assumed to be very large. We wish to estimate the height h = h (k) for which the conditional probability P τ 0,0 = k | H = h is maximal, with:

P τ 0,0 = k | H = h = q h P ( τ 0,h = k -1 | τ 0,h < τ 0,0 ) = q h P 0 X k-1 = h | τ 1,h < τ 1,0 .
We have

P τ 0,0 = k | H = h ≈ q h k-1 l=1 (1 -q x l )
where x l is the position at time l of the walker X which can roughly be estimated to be of order x l ≈ l h k . For large x, q x is small and of order γx -(λ+1) .Thus

P τ 0,0 = k | H = h ≈ h -(λ+1) 1 -γ k l=1 (x l ) -(λ+1) ≈ h -(λ+1) 1 -γ h k -(λ+1) k l=1 l -(λ+1) ≈ h -(λ+1) 1 -γkh -(λ+1)
Maximizing this over h we get h (k) ≈ k 1/(λ+1) suggesting that the height of a big excursion should scale like its length raised to the power 1/ (λ + 1) . Plugging h (k) in the latter expression

P τ 0,0 = k | H ≥ h (k) ≈ k -1 .
Thus

P ( τ 0,0 = k) ≈ P τ 0,0 = k, H ≥ h (k) ≈ k -1 h (k) -γ = k -(γ/(λ+1)+1) .
Scaling properties. We now describe the scaling relationship between height and length of an excursion. Because P H ≥ h ∼ h -γ , we conclude that the excursions lengths should also display power-law tails with

P ( τ 0,0 > k) ∼ k -γ/(λ+1) .
Thus τ 0,0 and H λ+1 are expected to be tail equivalent with P ( τ 0,0 ≥ k) ≈ P H λ+1 ≥ k , meaning that the ratio P ( τ 0,0 ≥ k) /P H λ+1 ≥ k should tend to some constant when k → ∞.

Recall that the exact p.g.f. of τ 0,0 is:

φ 0 (z) = E z τ 0,0 = z φ 1 (z) = z 2 e 1 I -z P -1 q,
where P is an upper-diagonal substochastic matrix with upper-diagonal entries p x and diagonal entries r x . So this matrix is diagonable with eigenvalues r x = 1 -ρ x .

The spectral gap here is ρ x = x -λ → 0 as x → ∞, algebraically fast, so the tails of τ 0,0 are not expected to be geometric, rather they are power-law. We have

k z k P ( τ 0,0 > k) =: Φ 0 (z) = 1 -φ 0 (z) 1 -z = 1 z (1 -z) h (z)
,

where h (z) ∼ z→1 (1 -z) -a and a = γ/ (λ + 1) . Thus Φ 0 (z) ∼ z→1 (1 -z) -(1-a) .
From the Karamata theorem (see [START_REF] Karamata | Sur un mode de croissance régulière. Théorèmes fondamentaux[END_REF]), therefore

z k Φ 0 (z) = P ( τ 0,0 > k) ∼ k→∞ k -a /Γ (1 -a) . (11)
Note that in the positive recurrent case (γ -λ > 1), a > 1 so τ 0,0 has finite mean, whereas in the null recurrent case, a < 1, and E ( τ 0,0 ) = ∞. We shall return to this crucial point later.

First passage times. Let τ x,y be the first passage time at y = x when the process is started at x. We wish here to briefly derive an exact formal formula for the law of τ x,y , making use of the Green function of a MC. Let

φ x,y (z) ≡ ∞ k=1 z k P ( τ x,y = k)
be the generating function of the law of τ x,y . Then, with

g z (x, y) ≡ ∞ n=0 z n P x X n = y = ∞ n=0 z n P n (x, y) = I -z P -1 (x, y)
the generating function of P n (x, y) (the Green potential function of the chain), using P n (x, y) = n m=0 P ( τ x,y = m) P n-m (y, y), we easily get the expression:

φ x,y (z) = g z (x, y) g z (y, y) .
In particular, φ x,0 (z) = g z (x, 0) g z (0, 0) and φ 0,x (z) = g z (0, x) g z (x, x) are the generating functions of τ x,0 and τ 0,x .

The p.g.f. φ x,x (z) first-return time τ x,x to state x satisfies

φ x,x (z) = g z (x, x) -1 g z (x, x) = 1 - 1 g z (x, x) where g z (x, x) = ∞ n=0 z n P x X n = x = ∞ n=0 z n P n (x, x) is the Green function at x.

First return times to 0 and contact probability at 0

First return times to 0 and contact probability at 0 are classically related to one another by renewal arguments. Let us here compute the contact probability at 0 making use of the first return time probability to the origin. With |z| < 1, let [START_REF] Ferrari | Construction of stochastic processes, coupling and regeneration[END_REF] φ 0 (z) = Because the Markov chain X n is always recurrent, u 0 (1) = 1+ n≥1 P 0 X n = 0 = ∞.

Recalling P ( τ 0,0 > k) ∼ k -γ/(λ+1) , by the Karamata's theorem, we get φ 0 (z)

∼ z↓1 1 -(1 -z) γ/(λ+1) . Thus u 0 (z) ∼ z↓1 (1 -z) -γ/(λ+1)
We conclude from singularity analysis and Karamata's theorem applied to u 0 (z) that

• If 0 < γ -λ < 1 (null recurrence): P 0 X n = 0 ∼ n→∞ n -(1-γ/(λ+1))
with 1 -γ/ (λ + 1) ∈ (0, 1) . The contact probability at 0 decays algebraically to 0, with decay exponent 1 -γ/ (λ + 1).

• If γ -λ > 1 (positive recurrence):

P 0 X n = 0 → n→∞ π 0 < ∞,
the mass at 0 of the invariant probability measure. In the positive recurrent setup, the return time to the origin only occurs in finite time ( τ 0,0 < ∞ with probability 1). Recall also that E ( τ 0,0 ) < ∞. In fact, a more detailed study of the singularities of φ 0 (z) shows that τ 0,0 has no moment of order larger or equal than γ λ+1 .

Excursion statistics of the special Markov chain: extreme value analysis heuristic

In this Section, we study the way the mean and median of both the empirical average and sample maximum grow, depending on the values of (γ, λ). Because growth is rare and regularly obstructed by catastrophic events, we expect growth to be slow. In fact, we will show that growth is algebraic only, with exponents that we can identify in terms of the model parameters (γ, λ).

Recall the height of an excursion obeys ( 15)

P H ≥ h = 1/ϕ (h) ∼ h -γ .
Therefore, E H = ∞ when γ ∈ (0, 1], E H < ∞ as soon as γ > 1 and

E H 2 < ∞ as soon as γ > 2.
For the excursions lengths, we have ( 16) λ+1) .

P ( τ 0,0 > k) ∼ k -γ/(
Therefore E ( τ 0,0 ) = ∞ when γ -λ ≤ 1, E ( τ 0,0 ) < ∞ as soon as γ -λ > 1 and E τ 2 0,0 < ∞ as soon as γ -λ > 2. We will summarize these results as follows: Proposition 3. With the scale function given by [START_REF] Eliazar | Growth-collapse and decay-surge evolutions, and geometric Langevin equations[END_REF], the law of the excursion height of the special walker is exactly given by:

(17) P H ≥ h = 1/ϕ (h) , satisfying P H ≥ h ≈ h -γ .
(i) when γ -λ ∈ (0, 1] and γ ∈ (0, 1], the special MC is null recurrent with both

E ( τ 0,0 ) = ∞ and E H = ∞.
(ii) when γ-λ ∈ (0, 1] and γ > 1, the special MC is null recurrent with E ( τ 0,0 ) = ∞ but with E H < ∞.

(iii) In the positive recurrent case (γ -λ > 1) the special MC sample paths are made of infinitely many i.i.d. excursions with both E ( τ 0,0 ) and E H < ∞.

H (n)
where H (n) d = H are i.i.d. with law governed by [START_REF] Karamata | Sur un mode de croissance régulière. Théorèmes fondamentaux[END_REF]. Due to [START_REF] Karlin | Occupation time laws for birth and death processes[END_REF], there exists a sequence h N such that N µ P H > h N → N ∞ α for some α > 0, say α = log 2.

From [START_REF] Karamata | Sur un mode de croissance régulière. Théorèmes fondamentaux[END_REF], we obtain:

(18) h N ≈ N 1/γ .
Using this, we have

P (H N ≤ h N ) = 1 -P H > h N N/µ → N ∞ e -α = 1/2,
and therefore (when γ -λ > 1), the typical (median) maximal height that the MC reaches by time N, grows like N 1/γ . With M denoting the median value, with M N := max X 1 , .., X N , we therefore have:

M (M N ) = h N ≈ N 1/γ .
Clearly, the number of

H n , n = 1, .., [N/µ] exceeding h N converges to a Poisson(α) distributed random variable. Furthermore, H N /h N d → N ∞ F where F is Fréchet (α, γ) distributed with P (F ≤ x) = e -αx -γ .
(Length). Similarly for the lengths, let

τ N = max n=1,..,[N/µ] τ (n) 0,0
be the length of the largest excursion, with τ (n) 0,0 d = τ 0,0 i.i.d. with law governed by [START_REF] Kella | On growth-collapse processes with stationary structure and their shot-noise counterparts[END_REF]. Due to [START_REF] Kella | On growth-collapse processes with stationary structure and their shot-noise counterparts[END_REF], there exists a sequence of time lags k N such that

N µ P ( τ 0,0 > k N ) → N ∞ α > 0. Since P ( τ 0,0 > k) ≈ k -γ/(λ+1) , we get: (19) k N ≈ N (λ+1)/γ .
We have

P (τ N ≤ k N ) = (1 -P ( τ 0,0 > k N )) [N/µ] → N ∞ e -α = 1/2.
Therefore, when γ -λ > 1, the typical length of the excursion with maximal length by time N grows like N (λ+1)/γ N. Note that (20)

h N ≈ k 1/(λ+1) N
, so that the typical height of the largest excursion scales like its length raised to the power 1/ (λ + 1).

Note also that, observing

E (M k N ) ≈ M (M N ), we obtain E (M N ) ≈ M max X 1 , .., X k -1 N = N γ/(λ+1) 1/γ = N 1/(λ+1)
with M (M N ) E (M N ) in the positive-recurrent regime under study γ -λ > 1.

Let τ 0,h be the first time at which some excursion height exceeds the level h (the time between failure at h). We have

P (τ 0,h > N ) = P (H N ≤ h) = 1 -P H > h [N/µ]
.

Due to [START_REF] Karamata | Sur un mode de croissance régulière. Théorèmes fondamentaux[END_REF], for all α > 0, assuming h large, we get:

(21) P h -γ τ 0,h > α ∼ h large 1 -h -γ α µ h γ → e -α µ ,
showing that τ 0,h is of order M (τ 0,h ) ≈ h γ with an exponential limit law. This point is in accordance with [START_REF] Martin | Batch queues, reversibility and first-passage percolation[END_REF]. We can summarize these results as follows:

Proposition 4. Assume γ -λ > 1 (positive recurrence of the special walker). Then, (i) the typical (median) height h N of its largest excursion satisfies h N ≈ N 1/γ .

(ii) the typical (median) length k N of its largest excursion satisfies k N ≈ N (λ+1)/γ , so with h N ≈ k 1/(λ+1) N .

(iii) the typical (median) time to failure at level h satisfies: M (τ 0,h ) ≈ h γ .

7.2. The null recurrent case. To some extent, this situation extends to the range γ -λ < 1, although it deserves a special treatment. In this null recurrent case indeed, µ = ∞ and so one deals with very large and therefore rare excursions (see [START_REF] Huillet | On rare and extreme events[END_REF]). By renewal arguments, the expected number of such excursions by time N (large) now is of order N a /c (with a = γ/ (λ + 1) < 1 and c > 0), so much smaller than N . The typical length k N of the largest excursion by time N is now given by:

N a c P ( τ 0,0 > k N ) → N ∞ α > 0,
leading to k N ≈ N (with no dependence of the scaling exponent on γ and λ): In this null-recurrent regime, the size of a typical excursion is the largest possible, corresponding to a single big excursion (or perhaps just a few of them).

Similarly, the maximal height H N reached by time N in this null recurrent case is again given by:

H N = max n=1,..,[N a /c] H (n) ,
where H (n) are i.i.d. with law governed by [START_REF] Karamata | Sur un mode de croissance régulière. Théorèmes fondamentaux[END_REF]. Due to [START_REF] Eliazar | Growth-collapse and decay-surge evolutions, and geometric Langevin equations[END_REF], there exists a sequence h N such that N a c P H > h N → N ∞ α for some α > 0, leading to

(22) h N ≈ N 1/(λ+1) . Thus, when γ -λ < 1, h N ≈ k 1/(λ+1) N
. Although when γ -λ < 1, h N and k N are not individually of the same order of magnitude as when γ -λ > 1, the typical height of the largest excursion continues to scale like its length raised to the power 1/ (λ + 1).

In the null recurrent case, we also have

P (τ 0,h > N ) = P (H N ≤ h) = 1 -P H > h [N a /c] .
Due to [START_REF] Karamata | Sur un mode de croissance régulière. Théorèmes fondamentaux[END_REF], for all α > 0, when h is large, we get:

(23) P h -(λ+1) τ 0,h > α ∼ h large 1 -h -γ α a h γ /c → e -α a /c ,
showing that the time to failure τ 0,h at h is now of order M (τ 0,h ) ≈ h λ+1 with a Weibull limit law. We shall summarize these results as follows:

Proposition 5. Assume γ -λ < 1 (null recurrence of the special walker). Then, (i) the typical (median) height h N of its largest excursion satisfies h N ≈ N 1/(λ+1) .

(ii) the typical (median) length k N of its largest excursion satisfies k N ≈ N, so still with h N ≈ k 1/(λ+1) N .

(iii) the typical (median) time to failure at level h satisfies: M (τ 0,h ) ≈ h λ+1 .

From this result concerning the null recurrent case and the previous one in the positive recurrent case, we observe that the typical height of the largest excursion always scales like its length raised to the power 1/ (λ + 1).

7.3. Height and length of excursions. Scaling properties in the regime 1 < γ -λ < 2. The previous study in the positive recurrent case, suggests that, at least when 1 < γ -λ < 2, it could always be true that the expected height of any excursion (not only the largest) should scale like its length raised to the power 1/ (λ + 1). If this were to be the case, defining the empirical mean of the walk as (24)

X N ≡ 1 N N n=1 X n ,
then,

E X N = 1 N N k=1 k N µ P (τ 0,0 = k) k 1/(λ+1) ,
because there are approximately N µ P (τ 0,0 = k) excursions of length k whose contribution to the empirical mean is k N µ P (τ 0,0 = k), each with typical height k 1/(λ+1) . Using [START_REF] Fernandez | Coupling, renewal and perfect simulation of chains of infinite order[END_REF], we therefore expect that in the range 1 < γ -λ < 2 Note that we get a sub-diffusive (super-diffusive) regime in the domain {1 < γ -λ < 2}∩ {γ < (λ + 3) /2} (respectively {1 < γ -λ < 2} ∩ {γ > (λ + 3) /2}). In particular if λ > 1, the regime is always super-diffusive.

Rather considering the median value M X N of X N , still when 1 < γ -λ < 2, we get

M X N ≈ k N k=1 k 1+1/(λ+1) P (τ 0,0 = k) = k N k=1 k -γ-1 λ+1 ≈ k 1-γ-1 λ+1 N ≈ N λ-γ+2 γ
, with M X N E X N h N = M (M N ) for the parameter range 1 < γ -λ < 2 under concern. The median height M X N of X N can also be estimated in the following way:

(27) M X N ≈ h N x=1 xP X N = x ≈ h N x=1 xx -(γ-λ) ≈ h 2+λ-γ N ≈ N (2+λ-γ)/γ ,
where we used that when N is large, P X N = x is 'close' to the invariant measure π which scales for large x like x -(γ-λ) . Indeed, with P -

Q T V ≡ sup A {P (A) -Q (A)}
defining the total variation distance between two probability measures P and Q,

P X N = • -π • T V ≈ P (X N ≥ h N ) ≈ h 1+λ-γ N
≈ N -(γ-λ-1)/γ , from [START_REF] Martin | Batch queues, reversibility and first-passage percolation[END_REF]. The latter quantity N -(γ-λ-1)/γ tends to 0 because γ -λ > 1. This heuristic approach also suggests that the convergence to equilibrium of P X N = . should be algebraically slow when 1 < γ -λ < 2. To summarize, we obtained: Proposition 6. Let X N in (24) define the empirical mean height of the walker by time N .

(i) Assume 1 < γ -λ < 2. Then, the expected value E X N of X N satisfies: λ+1) . Its typical (median) value M X N satisfies: M X N ≈ N (2+λ-γ)/γ with, since M (M N ) ≈ N 1/γ and E (M N ) ≈ N 1/(λ+1) :

E X N ≈ N (2+λ-γ)/(
M X N E X N M (M N ) E (M N ).
(ii) When γ -λ > 2, E X N → N ∞ E X ∞ , the finite mean of the invariant measure π. However, M (M N ) ≈ N 1/γ and E (M N ) ≈ N 1/(λ+1) still holds true in this parameter range.

In (i) , it can be checked that the scaling exponent (2 + λ -γ) / (λ + 1) of E X N is indeed smaller than the one 1/γ of M (M N ) in the regime under study 1 < γ -λ < 2. Note also that one expects that, normalizing the maximum by its mean, M N /N 1/(λ+1) should admit a limit law as N → ∞.

(ii) is clear because, when γ -λ > 2, X N →

N →∞ E X ∞ < ∞ almost surely.

1 1

 1 k≥1z k P ( τ 0,0 = k) be the generating function of the first return time to zero probability. Let also[START_REF] Gripenberg | A stationary distribution for the growth of a population subject to random catastrophes[END_REF] u 0 (z) = 1 + n≥1 z n P 0 X n = 0 be the Green potential function of the chain at state x = 0. As can easily be checked by renewal arguments, u 0 (z) = 1 + u 0 (z) φ 0 (z), showing that (14) u 0 (z) = -φ 0 (z) and φ 0 (z) = 1 -1 u 0 (z) .

7. 1 .

 1 Height and length of the largest excursion: the positive recurrent case. (Height). The preceding informations allow to derive the following qualitative result about the maximal height H N reached by time N in the positive recurrent case (γ -λ > 1): by time N , with µ = E ( τ 0,0 ) < ∞ there are indeed N/µ i.i.d. excursions on average. Thus, H N = max n=1,..,[N/µ]

  the expected height of the special MC X at time N grows like N λ-γ+2 λ+1 .


G Salagnac

S Yovine

D Garbervetsky

Fast Escape Analysis for Region-based Memory Management

Keywords:

We present an algorithm for escape analysis inspired by, but more precise than, the one proposed by Gay and Steensgaard [11]. The primary purpose of our algorithm is to produce useful information to allocate memory using a region-based memory manager. The algorithm combines intraprocedural variable-based and interprocedural points-to analyses. This is a work in progress towards achieving an application-oriented trade-off between precision and scalability. We illustrate the algorithm on several typical programming patterns, and show experimental results of a first prototype on a few benchmarks.

Introduction

Garbage collection (GC) [START_REF] Jones | Garbage collection. Algorithms for automatic dynamic memory management[END_REF] is not used in real-time embedded systems. The reason is that temporal behavior of dynamic memory reclaiming is extremely difficult to predict. Several GC algorithms have been proposed for real-time embedded applications (e.g., [START_REF] Henriksson | Scheduling garbage collection in embedded systems[END_REF][START_REF] Siebert | Eliminating external fragmentation in a non-moving garbage collector for Java[END_REF][START_REF] Ritzau | Decreasing memory over-head in hard real-time garbage collection[END_REF][START_REF] Higuera | Memory management for real-time Java: an efficient solution using hardware support[END_REF]). However, these approaches are not portable (as they impose restrictive conditions on the underlying execution platform), do require additional memory, and/or do not really ensure predictable execution times.

An appealing solution to overcome the drawbacks of GC algorithms, is to allocate objects in regions (e.g., [START_REF] Tofte | Region-based memory management[END_REF]) which are associated with the lifetime of a computation unit (typically a thread or a method). Regions are freed when the corresponding unit finishes its execution. This approach is adopted, for instance, by the Real-Time Specification for Java (RTSJ) [START_REF] Bollella | The Real-Time Specification for Java[END_REF], where regions can be associated to runnables, and by [START_REF] Gay | Language support for regions[END_REF], which implements a library and a compiler for C. These region-based approaches define APIs which can be used to explicitly and manually handle allocation and deallocation of objects within a program. However, care must be taken when objects are mapped to regions in order to avoid dangling references. Thus, programming using such APIs is error-prone, mostly because determining objects' lifetime is difficult.

An alternative to programming memory management directly using an API consists in automatically transforming a program so as (a) to replace (whenever possible) "new" statements by calls to the region-based memory allocator, and (b) to place appropriate calls (i.e., guaranteeing absence of dangling references) to the deallocator. Such an approach requires to analyze the program behavior to determine the lifetime of dynamically allocated objects. In [START_REF] Deters | Automated discovery of scoped memory regions for real-time Java[END_REF], analysis is based on profiling, while [START_REF] Garbervetsky | Program instrumentation and run-time analysis of scoped memory in java[END_REF][START_REF] Cherem | Region analysis and transformation for Java programs[END_REF] rely on static (points-to and escape) analysis.

Escape analysis aims at conservatively determining if an object escapes from or is captured by a method. Intuitively, an object escapes a method when its lifetime is longer than the method's lifetime, so it can not be collected when the method finishes its execution. An object is captured by the method when it can be safely collected at the end of its execution.

Several approaches to escape analysis for Java have been proposed, most of which aim at allocating objects on the stack, and removing unnecessary synchronizations. [START_REF] Blanchet | Escape Analysis for Java(TM). Theory and Practice[END_REF] works on the bytecode, which brings in an additionnal complexity due to the stack-based model. [START_REF] Choi | Stack Allocation and Synchronization Optimizations for Java Using Escape Analysis[END_REF][START_REF] Whaley | Compositional pointer and escape analysis for Java programs[END_REF] use points-to analysis to determine if an object escapes a method through a path in the points-to graph. [START_REF] Gay | Fast escape analysis and stack allocation for object-based programs[END_REF] proposes a fast but very conservative escape analysis, based on solving a simple system of linear constraints obtained from a Static Single Assignment (SSA) form [START_REF] Cytron | Efficiently computing static single assignment form and the control dependence graph[END_REF] of the program.

For region-based allocation in Java, we are aware of two works. [START_REF] Garbervetsky | Program instrumentation and run-time analysis of scoped memory in java[END_REF] exploits method-call chains and escape analysis to dynamically map allocation sites to regions associated with methods. [START_REF] Cherem | Region analysis and transformation for Java programs[END_REF] defines a points-to analysis to determine regions of objects with similar lifetimes (with instruction-level resolution, as opposed to method-level).

In this paper, we present an algorithm for escape analysis inspired by, but more precise than, the one proposed in [START_REF] Gay | Fast escape analysis and stack allocation for object-based programs[END_REF]. The primary purpose of our algorithm is to produce useful information to allocate memory using a region-based memory manager. The algorithm combines intraprocedural variable-based and interprocedural points-to analyses. This is a work in progress towards achieving an application-oriented trade-off between precision and scalability. We illustrate the algorithm on several typical programming patterns, and show experimental results of a first prototype on a few benchmarks.

The algorithm

In this section we describe our escape analysis algorithm in detail. We assume the program is in static single assignment form (SSA) [START_REF] Cytron | Efficiently computing static single assignment form and the control dependence graph[END_REF], that is, every variable is assigned only once in the program. The transformation of the program into SSA comes at a cost, but gives to a flow-insensitive analysis the power of a flow-sensitive one. Our algorithm is mainly based on local variables, instead of on a complex points-to graph, which would be much more expensive to build and to work with. The analysis is based on abstract interpretation [START_REF] Cousot | Abstract interpretation and application to logic programs[END_REF] and computes several properties for local variables and methods.

Properties

escape

For each local variable v of a method, escape(v)∈Escape, where Escape is the lattice in figure 1(a), says whether v may escape from its method, that is, if an object pointed to by v is referenced in a way such that its lifetime may exceed the method.

A variable v escapes because it is returned (escape(v)=RETURNED) or it is copied into a global variable (escape(v)=STATIC). When a variable is stored into an object field (escape(v)=FIELD), v may escape through a chain of references. Determining whether v escapes in this case, requires further analysis that will be explained later. The value stands for variables that escape by several ways, or when the analysis cannot compute a tighter information (e.g., when v is used as a parameter in a non-analyzed method). For example, in the program shown on figure 1 Notice that escape(v) = ⊥ is not sufficient to say that the object pointed to by v is local to the method. It only means that the method does not create any new reference path from the outside of the method to the object, but the object may already be reachable from outside. This is the case for variable b in figure 1

mfresh

Let MFresh be the lattice: ⊥ ≤ RETURNED ≤ . For each method m, mfresh(m) ∈ MFresh describes how objects returned by m escape: mfresh(m)=⊥ when m does not return any object (it may be void, or return some primitive type value); mfresh(m)= when returned values are already known to escape from m in a different way; mfresh(m)=RETURNED when m returns an object (or several objects) which does (do) not escape otherwise. If there is no other path leading to this object (see section 2.2.2), the caller of m can capture it.

sites

Let Sites be P(AllocationSites∪{UNKNOWN}), where AllocationSites is the set of all allocation sites in the program. For each local variable v, sites(v)∈Sites contains all allocation sites that can create an object referenced by v. sites(v) can always be computed at the unique (thanks to SSA) statement where v is defined. To be conservative, if we cannot determine all the sites that v can point to (e.g., because of a not analyzed method call), a "fake" allocation site UNKNOWN is added to sites(v). In the program of figure 1 Notice that, if mfresh(m) = RETURNED, then objects from msites(m) are possibly captured by callers of m, but it is not certain. In some complex situations, there can still be a path of references leading to these objects. For example in the program shown on figure 6(a), the e variable is not captured by m0.

2.1.5 isdereferenced isdereferenced(v) is true iff v, or one of its aliases, is dereferenced in m. That is, is v.f appears in the right-hand side of an assignment.

2.1.6 usedasparameter usedasparameter(v) is true iff v, or one of its aliases, is used as a concrete parameter in a method call.

def

For each variable v, def(v) says how v was defined.

2.1.8 fielduse fielduse shows reference relations between local variables. For each v in m, fielduse(v) is the set of variables u in m such that v may be an alias of u.f (for some field f). fielduse is mainly useful when a variable v escapes by a FIELD: for example, if escape(v)=FIELD, but all variables of fielduse(v) are captured by m, then so is v.

the mrefs graph

When objects are passed through several methods, knowledge about local variables is often not sufficient to determine objects' lifetimes, that's why a reference graph is needed. Our reference graph is very simple, in order to minimize the algorithmic cost of the analysis. mrefs is a subset of AllocationSites× F ields × AllocationSites, where (α, f, β) ∈ mrefs means: "an object created in α, may point, with its f field, to an object created in β".

side

The main goal of our analysis is to determine in which regions to allocate objects. Each method m has an associated region, containing objects which do not escape m. To determine the region, we compute for each variable v of m, where objects pointed to by v live, namely, side(v):

• side(v)=INSIDE, when objects pointed to by v are captured by m. If they are created by m, they can be allocated in m's region. If they are created by callees, m can ask for them to be allocated in its region, as is described in [START_REF] Garbervetsky | Program instrumentation and run-time analysis of scoped memory in java[END_REF];

• side(v)=OUTSIDE, when objects pointed to by v live longer than m. If they are created by m, they must be allocated outside its stack frame. But such an object may be captured by a caller n of m, in this case m can allocate the object in n's region.

An example is presented on figure 6(a): the RefObject allocated by m2 is captured by m1. Our analysis detects this situation by computing side(a)=OUTSIDE and side(c)=INSIDE.

The rules

The algorithm works in two phases. First, it determines for each variable the values of escape, sites, isdereferenced, usedasparameter, fielduse, def, it builds the mrefsgraph and computes msites and mfresh values. To compute these values, the algorithm solves the least fixpoint in Figures 2 and3.

In a second phase, the algorithm uses these values to compute, for each variable, its side value, as presented on figure 4. It is the combination of side and sites that will enable us to instrument the bytecode in order to use a region memory allocator for captured sites.

First phase

Most of these rules are simple. They are only intraprocedural information propagation. The only complicated rule is the one on figure 3, which handles method calls. This is not trivial, because we do not want to perform a full

α: v := new α ∈ sites(v) v:= ϕ(v 1 ..v n ) def(v) = PHI ∀i = 1..n sites(v) ⊇ sites(v i ) escape(v) escape(v i ) escape(v i ) escape(v) isdereferenced(v) ≥ isdereferenced(v i ) isdereferenced(v i ) ≥ isdereferenced(v) usedasparameter(v) ≥ usedasparameter(v i ) usedasparameter(v i ) ≥ usedasparameter(v) fielduse(v) ⊇ fielduse(v i ) fielduse(v i ) ⊇ fielduse(v) v := v 1 def(v) = COPY other properties: similar to ϕ-expression v 1 .f := v escape(v) FIELD fielduse(v) v 1 mrefs ⊇ {s 1 f -→ s 2 , s 1 ∈ sites(v 1 ), s 2 ∈ sites(v 2 )} s := v escape(v) STATIC mrefs ⊇ {UNKNOWN -→ s, s ∈ sites(v i )} v := s def(v) = STATIC sites(v) UNKNOWN v := p def(v) = PARAM sites(v) UNKNOWN other properties: similar to ϕ-expression v := constant def(v) = CONSTANT sites(v) UNKNOWN v := v 1 .f def(v) = FIELD isdereferenced(v 1 ) ≥ true sites(v) ⊇ {s |∃ s ∈ sites(v 1 ), s f -→ s } If UNKNOWN ∈ sites(v 1 ) sites(v) UNKNOWN return m v escape(v) RETURNED mfresh(m) escape(v) msites(m) ⊇ sites(v)
Fig. 2. Escape analysis rules variable's field (e.g. by a u.f=v), v cannot be captured unless u is.

v := v 0 .m(v 1 ..v n ) ∀ m that may be invoked here If istobeprocessed(m) sites(v) ⊇ msites(m) If mfresh(m) = RETURNED escape(v) mfresh(m) ∀i = 0..n usedasparameter(v i ) ≥ true Let p i the i-th formal parameter of m isdereferenced(v i ) ≥ isdereferenced(p i ) If ¬escape(p i ) ∈ {RETURNED, ⊥} escape(v i ) mrefs ⊇ {UNKNOWN -→ s, s ∈ sites(v i )} If isdereferenced(p i ) = true mrefs ⊇ {UNKNOWN -→ s |∃ s ∈ sites(v i ), s -→ s)} else sites(v) UNKNOWN ∀i = 0..n usedasparameter(v i ) ≥ true isdereferenced(v i ) ≥ true escape(v i ) mrefs ⊇ {UNKNOWN -→ s, s ∈ sites(v i )}

Examples

Let us consider the example presented in fig. 5(a). First, m0 builds a small chained structure, then it calls m1 which makes the last element (t3) escape.

As shown on fig. 5(b), the analysis of m0 understands the behavior of m0, but as we can only match x with t1, and not a with t2, we cannot keep track of m1. Nevertheless, to stay conservative, we put an edge from UNKNOWN to the site of t2 because x is dereferenced in m1. Notice that, t2 and t3 are usedasparameter, because they are the this parameter of their constructor.

Empirical results

We have implemented a prototype version of this algorithm using the Soot framework [START_REF] Sundaresan | Soot -A java optimization framework[END_REF] v.2.2.1. Table 1 presents the results of our algorithm on the Jolden benchmarks [START_REF] Cahoon | Data flow analysis for software prefetching linked data structures in java controller[END_REF]. The first two column are the size of the program in lines, and the number of allocation sites. The next three columns present the time spent by our escape analysis, in seconds, not including Soot's phases:

  (b) escape(a)=STATIC, escape(b)=⊥, escape(c)=RETURNED, escape(d)= .

Fig. 1 .

 1 Fig. 1. The Escape lattice and the Test01 program

  (b), sites(a) = sites(c) = {[m1:c=new Object]}, and sites(b) = sites(d) = {[m2:d=new Object]}. 2.1.4 msites For each method m, msites(m) is an element of Sites, saying where objects returned by m come from. In the program of figure 1(b), msites(m0) = ∅, msites(m1) = {[m1:c=new Object]}, and msites(m2) = {[m2:d=new Object]}.

Fig. 3 .

 3 Fig. 3. Escape analysis rules (cont)

Fig. 5 .

 5 Fig. 5. The Test25 program

  (b) which is an alias of the static variable s.

			class Test01 {
			static Object s,t;
			void m0() {
				Object a=m1();
				s=a;
	↑		}	Object b=m2();
			Object m1() {
	FIELD RETURNED STATIC			Object c=new Object(); return c;
	↑		} Object m2() {
				Object d=new Object();
	⊥			s=d; return d;
	(a) The escape lattice	}	}
		(b) the Test01 program

  That is why the only captured site is [m0:t1 = new RefObject].

	class RefObject {
	Object f;
	}
	class Test25 {
	void m0() {
		RefObject t1=new RefObject();
		RefObject t2=new RefObject(); Object t3=new Object();		UNKNOWN
		t1.f=t2; t2.f=t3;						t2=new	f	t3=new
		m1(t1);
	} void m1(RefObject x) static Object s;				t1=new	f
	{	RefObject a=(RefObject)x.f;				(b) mrefs graph
		Object b=a.f;
		s=b;
	}
	}
	(a) The Test25 program
		escape mfresh	def	IsD	uP	fielduse		sites msites	side
	m0	⊥						∅
	t1	⊥	NEW	true	true	[]	[m0:t1 = new RefObject]	INSIDE
	t2	FIELD	NEW	false	true	[t1]	[m0:t2 = new RefObject]	OUTSIDE
	t3	FIELD	NEW	false	true	[t2]	[m0:t3 = new java.lang.Object]	OUTSIDE
	m1	⊥						∅
	x	⊥	PARAM	true	false	[]		[UNKNOWN]	OUTSIDE
	a	⊥	FIELD	true	false	[]		[UNKNOWN]	OUTSIDE
	b STATIC	FIELD	false	false	[]		[UNKNOWN]	OUTSIDE
					(c) analysis results

Ackgnoledgements

We thank Chaker Nakhli and the anonymous referees for their helpful remarks.

 [START_REF] Cahoon | Data flow analysis for software prefetching linked data structures in java controller[END_REF]Partially supported by projects DYNAMO (Min. Research, France) and MADEJA (Rhône-Alpes, France). 4 Partially supported by projects ANCyT grant PICT 11738 and IBM Eclipse Innovation Grants.

there are typographic errors in the AIOOL'05 version

points-to analysis, neither to be too conservative about method calls.

Our analysis is designed to process arbitrary portions of an application. That is why we have an istobeprocessed predicate, that tells if a method must be analyzed or not. If not, for example because the method is native, or unavailable, we must be conservative about it.

For a not analyzed method, we assume that all parameters escape, and are referenced by the UNKNOWN site.

On the other hand, if the method is analyzed, then we can be more precise. Obviously, we have sites(v) ⊇ msites(m), that is, v will point to any object returned by m. If these objects have escaped (mfresh(m) = RETURNED), then the return value is not capturable either. (escape(v) mfresh(m))

To process the parameters of m, we match the formal parameters (p i ) with the concrete ones (v i ): if p i escapes from m, v i is considered as escaping from the current method, and we put an edge from UNKNOWN to all sites pointed to by v i . If p i does not escape but isdereferenced in m, then we cannot be precise about those references without performing a points-to analysis. In this case, we conservatively consider that all children of v i escape.

Second phase

Once the fixed point is reached, the algorithm computes side(v) for each variable using rules shown in figure 4. This is not a one-pass computation, but a second least fixpoint, because of the (1) and (2) rules:

• The (1) rule says that, if a variable may alias another, then those two variables cannot have different side values.

• Similarly, the (2) rules says that if a variable v is referenced by another class loading, transformation from bytecode to Jimple (Soot's three-address stackless code), and transformation into SSA form. The last three columns give the number of INSIDE variables and allocation sites, as computed by our algorithm, and the number of stackable variables, as computed by our implementation of G&S's analysis [START_REF] Gay | Fast escape analysis and stack allocation for object-based programs[END_REF]. Our analysis is more precise than [START_REF] Gay | Fast escape analysis and stack allocation for object-based programs[END_REF] as it subsumes all its rules. That is, all stackable variables in the sense of [START_REF] Gay | Fast escape analysis and stack allocation for object-based programs[END_REF] are INSIDE variables, but the converse is not true. In our experiments, we did not use any inlining of analyzed code. It is interesting to remark that without inlining, [START_REF] Gay | Fast escape analysis and stack allocation for object-based programs[END_REF] does not find any stackable variable in the programs of figures 5 and 6. As noted in [START_REF] Gay | Fast escape analysis and stack allocation for object-based programs[END_REF], both analyses will benefit from method inlining.

We did not have enough time to use the computed information to actually instrument the benchmarks as described in [START_REF] Garbervetsky | Program instrumentation and run-time analysis of scoped memory in java[END_REF]. We count on doing this soon. Anyway, a preliminary implementation on another test program revealed a gain of 20% of total utilized memory, when using GC together with regionbased manager, w.r.t. GC only, even if the actual region-allocated memory is about 5%.

Besides, only a subgraph of the whole call graph has been analyzed for each test case. The subgraph contains all application methods and a subset of library methods transitively invoked by the program. This explains why there are only a few allocation sites. Nevertheless, these results are interesting, because an important fraction of analyzed allocation sites are indeed computed to be captured. Our algorithm is parameterized by the set of classes to be analyzed. This allows the user to fine-tune the analysis trading precision against performance according to specific application behaviors.