We consider the recursive equation "x(n + 1) = A(n) ⊗ x(n)" where x(n + 1) and x(n) are R k -valued vectors and A(n) is an irreducible random matrix of size k × k. The matrix-vector multiplication in the (max,+) algebra is defined by

This type of equation can be used to represent the evolution of Stochastic Event Graphs which include cyclic Jackson Networks, some manufacturing models and models with general blocking (such as Kanban). Let us assume that the sequence {A(n), n ∈ N} is i.i.d or more generally stationary and ergodic. The main result of the paper states that the system couples in finite time with a unique stationary regime if and only if there exists a set of matrices C such that P {A(0) ∈ C} > 0 and the matrices C ∈ C have a unique periodic regime.

Introduction

Let us consider the following recursive equation:

x i (n + 1) = max j (A ij (n) + x j (n)) x i (0) = (x 0 ) i (1)
The sequences {A ij (n)} are given (exogenous data). The process we want to study is the sequence of vectors {x(n) = (x 1 (n), . . . , x i (n), . . . ) ′ }. The vector x 0 is the initial condition.

Because of the generality of Equation ( 1), it appears in many different types of applications.

In fact, this equation appears in statistical mechanics in the study of crystal structures, see Griffiths [START_REF] Griffiths | Frenkel-Kontorova models of commensurate-incommensurate phase transitions[END_REF]. It is also common under one form or another in economic and control literature, see Yakovenko & Kontorer [START_REF] Yakovenko | Nonlinear semigroups and infinite horizon optimization[END_REF]. In fact, it is the basic Bellman equation of dynamic optimization in discrete time, for a finite state space. Recently, this kind of equation has become very popular in the study of Discrete Events Dynamic Systems (DEDS), see for example the recent textbooks of Baccelli, Cohen, Olsder & Quadrat [START_REF] Baccelli | Synchronization and Linearity[END_REF] and Glasserman & Yao [START_REF] Glasserman | Monotone Structure in Discrete-Event Systems[END_REF]. Let us insist on the last point.

A large class of computer or communication networks accepts a description as DEDS. Different approaches have been proposed to model DEDS. Petri Nets is one of the most common formalism. More precisely, a sub-class of Petri Nets, Event Graphs, appears to be very efficient in describing models with synchronization, blocking and/or fork-join properties. We can mention Job-Shop models (see Cohen, Dubois, Quadrat & Viot [START_REF] Cohen | A linear system-theoretic view of discreteevent processes and its use for performance evaluation in manufacturing[END_REF] or [START_REF] Baccelli | Synchronization and Linearity[END_REF]), cyclic Jackson Networks (see Section 2.2) or asymmetric exclusion models as examples. On the other hand, Event Graphs cannot be used to model systems with routings. We can describe the evolution of an Event Graph by the daters associated with the transitions (nodes) of the graph. It is well known, see for example [START_REF] Baccelli | Synchronization and Linearity[END_REF], that the evolution of the daters of an Event Graph can be represented in the form of Equation [START_REF] Anantharam | Stationary solutions of stochastic recursions describing discrete event systems[END_REF]. Generalized Semi-Markov Processes (GSMP) are another common formalism for DEDS. It is shown in [START_REF] Glasserman | Monotone Structure in Discrete-Event Systems[END_REF] that a GSMP with convex and homogeneous structural properties admits a representation of the form [START_REF] Anantharam | Stationary solutions of stochastic recursions describing discrete event systems[END_REF]. In order to give a flavor of the modeling power of Equation [START_REF] Anantharam | Stationary solutions of stochastic recursions describing discrete event systems[END_REF] and in order to motivate the practical interest of this work, we will present two different models in Section 2.

• The first model appears in the modeling and the analysis of parallel programs and architectures. It is a task graph with 'and' synchronizations also known as PERT Network.

• The second model is a closed cyclic Jackson Network. We will consider both infinite and finite buffers with various blocking modes.

It is very fruitful to use a matrix-vector notation for Equation [START_REF] Anantharam | Stationary solutions of stochastic recursions describing discrete event systems[END_REF]. We define the following "(max,+)" notations: ∀x, y ∈ R ∪ {-∞}, x ⊕ y = max(x, y), x ⊗ y = x + y .

We define also the k × k matrix A(n) = {A ij (n), i, j = 1, . . . , k} and the column vector x(n) = (x 1 (n), . . . , x k (n)) ′ . With these notations, the basic Equation (1) takes a very simple and convenient form. In fact, it can be rewritten as:

x(n + 1) = A(n) ⊗ x(n) .
(2)

The matrix-vector product is defined in a natural way just by replacing + and × by ⊕ and ⊗, i.e. (A ⊗ x) i = j A ij ⊗ x j = max j (A ij + x j ).

Historically, the first approach has been to consider deterministic systems where A ij (n) ≡ A ij .

It is natural to consider a stochastic extension where {A ij (n)} is a sequence of random matrices. As a consequence, here is an equivalent way of introducing our subject: it is a counterpart of the classical theory of products of random matrices (see Furstenberg & Kesten [START_REF] Furstenberg | Products of random matrices[END_REF] or Bougerol & Lacroix [START_REF] Bougerol | Products of Random Matrices with Applications to Schrödinger Operators[END_REF]) but in another algebraic structure, the (max,+) algebra.

For systems described by Equation ( 2), we will consider two kinds of asymptotic results.

• First order limits, on ratios:

lim n x(n) ∞ n , lim n x i (n) n . (3)
• Second order limits, on differences:

lim n x i (n + 1) -x i (n), ∀i, lim n x j (n) -x i (n), ∀i = j . (4)
A first order limit is a cycle time or equivalently the inverse of a throughput. Second order limits are related to waiting and idle times, queue length and frequency of occupation. More insights on the relations between these limits and quantities of interest for the system will be provided in Section 2. Our goal is to find stationary regimes for second order limits. Multiple stationary regimes will mean multiple possible regimes for waiting times or queue lengths, depending on the initial condition.

Among the systems modeled by Equation ( 2), we can distinguish two classes: the open (or nonautonomous) systems and the closed (or autonomous) ones. Open systems have been exhaustively treated by Baccelli [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF] [4] (for both first and second order limits). Problems of existence and uniqueness of first order limits for closed systems have been solved by Cohen [START_REF] Cohen | Subadditivity, generalized product of random matrices and operations research[END_REF] (see also [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF]). These results are recalled in §6.1. This paper deals with the open question of existence and uniqueness of second order limits for closed systems. These problems were considered in several earlier papers (Resing, de Vries, Hooghiemstra, Keane & Olsder [START_REF] Resing | Asymptotic behavior of random discrete event systems[END_REF] and [START_REF] Olsder | Discrete event systems with stochastic processing times[END_REF], Baccelli [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF]) but only sufficient conditions of uniqueness were known. The approach we use is new and exploits completely the common hidden algebraic structure of the different models we consider. It enables us to obtain necessary and sufficient conditions for stability (in some cases) together with simple proofs.

The conditions we give are based on the structure of deterministic matrices chosen in the support of the random matrix A(0). One of the main results states that the system couples in finite time with a unique stationary regime if and only if there exists a set of matrices C such that P {A(0) ∈ C} > 0 and the matrices C ∈ C have a unique periodic regime. The proof makes use of Borovkov's theory of renovating events, see Borovkov & Foss [10] [11]. This theory appears to be much more tractable than classical Harris regeneration due to the specific form of our recursive equations.

The paper is organized as follows. We introduce two models in Section 2, cyclic Jackson Networks and task graphs with random precedences. Sections 3, 4 and 5 are presenting the tools that we are using in the paper. They can be skipped by people knowing the subject. Section 3 is devoted to the (max,+) algebra, Section 4 to the spectral theory in this algebra and Section 5 to Borovkov's theory of renovating events. Section 6 presents the main results. In 6.1, we recall some results from [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF] and [START_REF] Cohen | Subadditivity, generalized product of random matrices and operations research[END_REF]. In 6.2, we state some preliminary results. In 6.3 and 6.4, we give sufficient conditions for the stability of discrete and general models respectively. In Section 7, we establish the converses of the results of the previous section. In Section 8, we weaken the assumptions under which some of our results apply. Finally, for convenience, some of the proofs are given in Appendix.

2 Two Motivating Models

Task graphs

We consider a parallel program executed on several identical processors. We model it by its precedence graph τ . If we consider a system of k processors, the graph τ has a set of nodes which is k × N. The node (i, n) represents the n-th task to be executed at processor i. The arcs between nodes represent the synchronization constraints. There is an arc between the node (i, n) and the node (j, m) (notation : (i, n) → (j, m)) if the n-th task at processor i has to be completed in order for the m-th task at processor j to be enabled. The execution of a task begins as soon as all the tasks of its incoming arcs are completed. Each task has a duration which may depend on the processor.

Let us consider a task graph with synchronizations only between consecutive levels n (i.e. nodes (1, n), . . . , (k, n)) and n + 1. We assume that the synchronizations depend on n. We denote by L(i, n) the set of nodes j such that (i, n) → (j, n + 1). We suppose that for all i, there exists a probability law P i on the subsets of (1, . . . , k) such that L(i, n) = (j 1 , . . . , j p ) with probability P i {(j 1 , . . . , j p )}. We denote by x i (n) the date of completion of task n at processor i, and by A ji (n) the duration of the synchronization constraint between nodes (i, n) and (j, n + 1) (it may include a transmission time as well as the execution time at processor j). We adopt the convention that

A ji (n) = -∞ if j ∈ L(i, n).
It is easy to check that such a model, we could call it a task graph with random precedences, verifies Equation [START_REF] Asmussen | On coupling and weak convergence to stationarity[END_REF].

A Queuing Network model studied by Baccelli & Liu [START_REF] Baccelli | On a class of stochastic recursive equations arising in queueing theory[END_REF] corresponds to this model. The task resource models studied in [START_REF] Brilman | Synchronisation by resources sharing : a performance analysis[END_REF] or [START_REF] Gaubert | Task resource models and (max,+) automata[END_REF] also have this kind of structure.

Cyclic Jackson network

We consider a closed Jackson Network. The study of such closed networks can be traced back to Gordon and Newell, [START_REF] Gordon | Closed queuing systems with exponential servers[END_REF]. In their original model, there is a given number of indistinguishable customers. The routing of the customers leaving a given queue is provided by a sequence of i.i.d. Bernouilli random variables. All the service times are exponential. They prove the existence of an explicit product form for the unique stationary distribution.

A natural generalization of the basic model is to consider i. We consider a restriction of the previous model. There are k queues and all customers have the same cyclic route (1, 2, . . . , k, 1), see Figure 1. We will denote this model by CJN for Cyclic Jackson Network, following the terminology of [START_REF] Kaspi | Regenerative closed queueing networks[END_REF].

Queue 1 Queue 2

./GI/1/∞ FIFO ./GI/1/∞ FIFO ./GI/1/∞ FIFO Queue k Figure 1: A Cyclic Jackson Network consisting of k queues.
In the following, the numbering of queues has to be understood modulo [k], for example queue (k + 2) is queue 2. We denote by {σ j (n), n ∈ N}, the sequence of service times at queue j.

Instead of describing the system by the workload or the queue length process, as is usually done, we propose to study this model by introducing the following variables. With each queue j, we associate a dater {x j (n), n ∈ N}. The variable x j (n) represents the date of completion of the n-th service at queue j. All variables of interest for the network can be derived from these daters and from the sequences of service times. More precisely, we have:

• Asymptotic throughput at queue j:

γ j = lim n-→+∞ n x j (n) .
• Idle time of queue j before the arrival of the n-th customer to visit queue j.

I j (n) = x j (n) -σ j (n) -x j (n -1)
.

• Workload at queue j at the instant of the arrival of the n-th customer to visit queue j. This customer comes from queue j -1. We suppose that it was the n ′ -th customer to visit queue j -1.

W j (n) = x j (n) -σ j (n) -x j-1 (n ′ ) .
The variables (γ i ) which are obtained as ratios of daters will be called first order variables. The ones (I j , W j ) which are obtained as differences of daters will be called second order variables. We are concerned with the problem of deriving necessary and sufficient conditions under which there is a unique stationary regime for both first and second order variables. In such a case, we say that our model is stable.

Suppose for the moment that there are exactly k customers. We suppose also that there is initially one customer in each queue. These assumptions together with the FIFO service discipline at each queue yields the following property. The n-th customer to visit queue j will be, at the next step of its route, the (n + 1)-th customer to visit queue j + 1. As a consequence, we have

           x 1 (n + 1) = max(x 1 (n), x k (n)) + σ 1 (n) x 2 (n + 1) = max(x 2 (n), x 1 (n)) + σ 2 (n) • • • x j (n + 1) = max(x j (n), x j-1 (n)) + σ j (n) • • • .
Using the (max,+) notation, this can be rewritten as:

x(n + 1) = A(n) ⊗ x(n), where A(n) =          σ 1 (n) -∞ • • • -∞ σ 1 (n) σ 2 (n) σ 2 (n) . . . -∞ -∞ . . . . . . -∞ . . . . . . . . . . . . . . . -∞ -∞ • • • -∞ σ k (n) σ k (n)          . ( 5
)
The initial condition is x(0) 0, where x i (0) is the remaining service time of the customer being served at queue i at time 0.

When the service times are deterministic, it is possible to obtain many asymptotic behaviors, depending on the initial condition x(0). In fact, initial delays between customers might never vanish. Therefore, it is possible to have several stable regimes for second order quantities (I j , W j , . . . ) including periodic ones. For stochastic systems, when the service times are random variables, it is still possible to have several stationary regimes if the system is not "stochastic enough". As an application of the results presented in this paper, we obtain the necessary and sufficient conditions for the existence of a unique stationary regime for this CJN. This model will be used as an illustration of the results throughout the paper (Examples 4.8, 6.9, 6.16, 7.7).

When there are less than k customers in the network, the system can be represented in the same way as previously. The only difference is that the structure of matrices {A(n)} is more complicated. When there are more than k customers, the trick consists in splitting queues. Each queue which has originally (p > 1) customers in its buffer is transformed into p queues with one customer per buffer. This is done by creating p -1 fictive queues with service times identically equal to zero. By doing this, one gets back to the previous case. The main difference is that we have represented our model by a (max,+) linear system of dimension greater than the original number of queues. For more details on these transformations, see [START_REF] Mairesse | Stabilité des systèmes à événements discrets stochastiques[END_REF].

We can also model CJN with finite buffers (the ./G/1/L case). Finite buffers imply the blocking of some customers. There are different possible types of blocking.

1. Blocking before service. The service begins at queue i only when the buffer at queue (i+1) is not full.

2. Blocking after service (of communication type). Service at queue i begins as soon as a customer is available. After completion of the service, if the buffer of queue (i + 1) is full, the customer starts another service at queue i.

3. Blocking after service (of manufacturing type). Service at queue i begins as soon as a customer is available. After completion of the service, the customer has to wait in queue i if the buffer of queue (i + 1) is full. It prevents another customer from being served at queue i.

4. Blocking after service (of Kanban type). The mechanism is the same as previously but there exists a finite intermediate buffer between queue i and (i+1). A customer completing its service at queue i enters this intermediate buffer if the buffer of queue (i + 1) is full and the intermediate buffer is not. It enables to serve a new customer at queue i.

Excepting the blocking after service of communication type, all these types of blocking can be considered. We can even consider different types of blocking for the different queues of the network. In all cases, we obtain a (max,+) linear representation for the network.

In the case of a CJN with i.i.d. general service times (./GI/1 servers), the classical method for studying the network is to consider the Markov chain formed by queue lengths and remaining service times and to apply Harris regeneration techniques. This method was first introduced for closed acyclic Jackson Networks by Borovkov [START_REF] Borovkov | Limit theorems for queueing networks[END_REF], [START_REF] Borovkov | Limit theorems for queueing networks[END_REF]. For closed Cyclic Jackson Networks, it is used by Bambos [START_REF] Bambos | On closed ring queueing networks[END_REF] and Kaspi & Mandelbaum [27] [28]. The method provides sufficient conditions of stability. Our approach is completely different and provides necessary and sufficient conditions of stability for CJN.

For a much more complete description of the systems modeled by Equation ( 2), the reader is referred to the textbook of Baccelli, Cohen, Olsder and Quadrat [START_REF] Baccelli | Synchronization and Linearity[END_REF].

(max,+) Algebra

Definition 3.1 ((max,+) algebra). We consider the semiring (R ∪ {-∞}, ⊕, ⊗). The law ⊕ is "max" and ⊗ is the usual addition. We set ε = -∞ and e = 0. The element ε is neutral for the operation ⊕ and absorbing for ⊗. The element e is neutral for ⊗. The law ⊕ is idempotent, i.e. a ⊕ a = a. (R ∪ {ε}, ⊕, ⊗) is an idempotent semiring, called a dioid. It is usually referred to as the (max,+) algebra. We shall denote it by R max .

We define the product spaces R k max , R k×k max . We define the product of a vector by a scalar:

a ∈ R max , u ∈ R k max , (a ⊗ u) i = a ⊗ u i .
Matrix product is defined in a natural way, replacing + and × by ⊕ and ⊗ respectively. Let

A, B ∈ R k×k max , (A ⊗ B) ij = max l (A il + B lj ) = l A il ⊗ B lj .
Matrix-vector product is defined in a similar way, see Section 1. In the rest of the paper, the notations "+,×" will stand for the usual addition and multiplication. Nevertheless, we will write ab for a ⊗ b whenever there is no possible confusion. For example, for

A ∈ R k×k max , A d = A ⊗d = A ⊗ • • • ⊗ A d factors .
Let us recall some definitions adapted to the R max algebra. Definition 3.2. The graph of a square matrix A is a directed graph having a number of nodes equal to the size of A. This graph contains an arc from i to j iff A ji = ε. The valuation of this arc is A ji .

Definition 3.3. A square matrix A is irreducible if: ∀i, j ∃n 0 | (A n ) ij > ε (or equivalently if its graph is strongly connected). Definition 3.4. A square matrix A is aperiodic if: ∃N, ∀n N, ∀i, j, (A n ) ij > ε.
Definition 3.5. Let (Ω, F, P ) be a probability space. A stochastic matrix {A(ω), ω ∈ Ω} has a fixed structure if

P (A ij = ε) = 1 or P (A ij = ε) = 0, ∀i, j. Definition 3.6 (PR k ). The projective space PR k is defined as the quotient of R k by the paral- lelism relation: u, v ∈ R k u ≃ v ⇐⇒ ∃a ∈ R such that u = a ⊗ v .
Let π be the canonical projection of R k into PR k .

For example (e, -1) ′ and (2, 1) ′ = (e + 2, -1 + 2) ′ are in the same parallelism class, i.e. are two representatives of the same element of PR k . We define in the same way PR k max , PR k×k max and PR k×k . We use the same notation π for the different canonical projections. We define a norm and a distance on PR k which we call the projective norm and distance. Definition 3.7. Let x ∈ PR k and u ∈ R k be a representative of x, i.e. π(u) = x. We define:

|x| P = max i u i -min i u i .
Let x, y ∈ PR k and u, v ∈ R k be two representatives of x and y respectively. We define:

d(x, y) = d(u, v) = |x -y| P = i (u i -v i ) ⊗ i (v i -u i ) .
We write either d(x, y) or d(u, v) with some abuse of notation.

The space (PR k , |.| P ) is an Euclidean space. In particular, it is complete. The norm |.| P corresponds to the L ∞ norm1 on the projective space PR k . Proposition 3.8. Let A ∈ R k×k max be an irreducible matrix. Let u, v be two vectors of R k . We have:

d(Au, Av) d(u, v) .
Proof. By definition, we have:

d(Au, Av) = i ((Au) i -(Av) i ) ⊗ i ((Av) i -(Au) i ) .
We define j(i) such that (Au

) i = j A ij ⊗ u j = A ij(i) ⊗ u j(i) .
Note that j(i) depends on A and u. We have:

i ((Au) i -(Av) i ) = i   ( j A ij ⊗ u j ) -( j A ij ⊗ v j )   = i   (A ij(i) ⊗ u j(i) ) -( j A ij ⊗ v j )   i (A ij(i) ⊗ u j(i) -A ij(i) ⊗ v j(i) ) = i u j(i) -v j(i) i u i -v i We obtain d(Au, Av) i (u i -v i ) ⊗ i (v i -u i ) i.e. d(Au, Av) d(u, v). ⊓ ⊔
There is no simple criterion to get a strict inequality. This monotonicity has to be interpreted as a synchronization property.

Definition 3.9. We consider A ∈ R k×k max . We set

D(A) = sup u,v∈R k d(Au, Av) .
We call D(A) the projective diameter of A.

It is easy to prove that D(A) is finite if and only if ∀i, j, A ij > ε. A matrix A can be considered as a "linear" (in the (max,+) sense) operator from PR k into PR k . It is a bounded operator if the (decreasing) sequence D(A n ) has a finite limit, i.e. if A is aperiodic (Def. 3.4).

Deterministic Spectral Theory

We recall some results of the deterministic spectral theory in the R max algebra. Theorem 4.2 is due to Cuninghame-Green [START_REF] Cuninghame-Green | Describing industrial processes with interference and approximating their steady-state behaviour[END_REF]. Versions of Theorem 4.4 were proved in [START_REF] Romanovskiȋ | Optimization and stationary control of discrete deterministic process in dynamic programming[END_REF] [19] and [START_REF] Gondran | Valeurs propres et vecteurs propres dans les dioïdes et leur interprétation en théorie des graphes[END_REF]. Under the form proposed here, the result is from [START_REF] Cohen | A linear system-theoretic view of discreteevent processes and its use for performance evaluation in manufacturing[END_REF]. Theorem 4.5 is due to Cohen, Dubois, Quadrat and Viot [START_REF] Cohen | Analyse du comportement périodique des systèmes de production par la théorie des dioïdes[END_REF] and [START_REF] Cohen | A linear system-theoretic view of discreteevent processes and its use for performance evaluation in manufacturing[END_REF]. A complete treatment of the spectral theory can be found in [START_REF] Baccelli | Synchronization and Linearity[END_REF].

We want to find non trivial solutions to the equation:

A ⊗ x = λ ⊗ x ,
where A ∈ R k×k is an irreducible matrix, x is a column vector (the "eigenvector") and λ is a real constant (the "eigenvalue"). We define also periodic regimes for the eigenvalue problem.

A periodic regime of period d is a set of vectors {x 1 , . . . ,

x d } of R k verifying Ax i = λx i+1 , i = 1, . . . , d -1 and Ax d = λx 1 . Definition 4.1. For each path ζ = {t 1 , t 2 , • • • , t j , t j+1 = t 1 }
, we define its average weight by:

p(ζ) = a t 1 t j ⊗ • • • ⊗ a t 3 t 2 ⊗ a t 2 t 1 /j ,
(here the division is the conventional one).

Theorem 4.2. There is a unique (non ε) eigenvalue, λ. It satisfies the relation

λ = max ζ p(ζ) ,
where ζ covers all the circuits of (the graph of ) A. We call also λ the Lyapunov exponent or the cycle time of A.

There might be several eigenvectors. A linear combination (in R max ) of eigenvectors is an eigenvector. An eigenvector has all its coordinates different from ε (due to the irreducibility assumption).

Definition 4.3. For a matrix A, we define:

Critical circuit A circuit ζ of A is said to be critical if its average weight is maximal, i.e. if p(ζ) = λ.
Critical graph It consists of the nodes and arcs of A belonging to the critical circuit(s).

For a general graph, we define :

Cyclicity The cyclicity of a strongly connected graph is the greatest common divisor of the lengths of all the circuits. The cyclicity of a connected graph is the least common multiple of the cyclicities of its maximal strongly connected subgraphs (s.c.s.).

We normalize a matrix by subtracting (in the conventional algebra) the eigenvalue to all the coordinates. The eigenvalue of a normalized matrix is e. For a normalized matrix A of size k, we define:

A + = A ⊕ A 2 ⊕ • • • ⊕ A k .
We check that A + ⊕ A k+1 = A + .

Theorem 4.4. Let A be a normalized matrix.

a. Critical columns A + .i , i belonging to the critical graph, are eigenvectors.

b. For i, j belonging to the critical graph, π(A + .i ) and π(A + .j ) are different iff they belong to two different s.c.s. of the critical graph. c. Every eigenvector of A writes as a linear combination (in R max ) of critical columns A + .i .

In R max , every irreducible matrix is cyclic in the sense of the following theorem.

Theorem 4.5. For an irreducible matrix A of size k and whose eigenvalue is λ, there exists integers d and M such that:

∀m M, A m+d = λ ⊗d ⊗ A m , (6)
furthermore the smallest d verifying the property is equal to the cyclicity of the critical graph of A. We call it the cyclicity of A.

A cyclicity greater than one will provide periodic regimes of period greater than one for the eigenvalue problem.

Proposition 4.6. An irreducible matrix has a unique eigenvector and no periodic regimes of period greater than one for the eigenvalue problem, if and only if its critical graph has a unique s.c.s. and its cyclicity is one. Such a matrix will be called a scs1-cyc1 matrix.

The proof follows from Theorems 4.4 and 4.5.

Definition 4.7 (rank). By analogy with classical linear algebra, we define the "rank" of a matrix A as the number of additively independent columns (resp. lines) of A. More precisely, let A .i denote the i-th column of A. Matrix A is of rank r if there exists J ⊂ {1, . . . , k} such that |J | = r and ∀i = j ∈ J , π(A .i ) = π(A .j ) and ∀i ∈ J , ∃α j , j ∈ J , such that

π(A .i ) = π[ j∈J α j ⊗ A .j ] .
Let A be a rank 1 matrix. Then A is a scs1-cyc1 matrix and verifies A 2 = λ ⊗ A (λ is the eigenvalue of A). The other way round, let A be a scs1-cyc1 matrix and M be defined as in Equation [START_REF] Bambos | On closed ring queueing networks[END_REF]. One can check that A M is a matrix of rank 1.

Example 4.8. [Cyclic Jackson Network 1] Let us consider a basic Cyclic Jackson Network as presented in Section 2.2. We suppose that the service times are deterministic, i.e σ j (n) ≡ σ j . We suppose also that the number of customers, k, is equal to the number of queues. Then we can consider the (max,+) matrix associated with the network, see Equation [START_REF] Baccelli | On a class of stochastic recursive equations arising in queueing theory[END_REF]. The graph associated with this matrix is constituted by the circuit (1, 2, . . . , k, 1) and the recycling loops (1, 1) to (k, k). Let us define I = {i | σ i = max j σ j }. There are two possible cases.

• If the cardinal |I| < k, then the critical graph of the matrix consists of the nodes i ∈ I and the arcs (i, i), i ∈ I. It implies that the matrix is scs|I|-cyc1.

• If |I| = k then the graph and the critical graph of the matrix coincide. It implies that the matrix is scs1-cyc1.

We conclude that the matrix is scs1-cyc1 if and only if |I| = 1 or k.

Borovkov's Renovating Events Theory

Borovkov's theory deals with the problem of regeneration in so-called "Stochastic Recursive Sequences". For a complete treatment, the reader is referred to Borovkov [START_REF] Borovkov | Asymptotic Methods in Queueing Theory[END_REF], Borovkov & Foss [START_REF] Borovkov | Stochastically recursive sequences and their generalizations[END_REF][START_REF] Borovkov | Two ergodicity criteria for stochastically recursive sequences[END_REF] or Brandt, Franken & Lisek [START_REF] Brandt | Stationary Stochastic Models[END_REF]. Let (Ω, F, P ) be a probability space. Let θ be a measurable map from (Ω, F) into itself such that P is θ-invariant and θ-ergodic. Let (E, E) and (G, G) be two Polish spaces (complete, separable metric spaces) equipped with their respective Borel σ-algebra.

Definition 5.1. We call Stochastic Recursive Sequence (SRS), a sequence {x(n)} of E-valued random variables defined by

x(n + 1) = f (x(n), a(n)), n 0, x(0) = x 0 ,
where {a(n)} is an exogenous sequence of G-valued random variables, stationary with respect to the shift θ. The function f is a measurable function from E × G into E. The vector x 0 ∈ E is the initial condition. In order to stress the value of the initial condition, we will sometimes denote the SRS by {x(n, x 0 )}. We talk of an i.i.d. SRS when the sequence a(n) is i.i.d. (an i.i.d. SRS is a Markov Chain and the converse is true).

Definition 5.2. We consider a SRS {x(n)}. We denote by F l the σ-algebra

F l = σ{a(n), n ∈ {-∞, . . . , l -1}}.
The sequence of events {A(n) ∈ F n+m , n ∈ N} is said to be a renovating sequence of length m and of associated function φ :

G m → E if: ∃n 0 , ∀n n 0 , x(n + m) = φ (a(n), a(n + 1), • • • , a(n + m -1)) on A(n) .
A sequence {A(n), n ∈ N} of renovating events of same length and associated function is said to be stationary if

A(n) = A(0) • θ n = θ -n A(0).
We need the following notions of convergence:

Definition 5.3. We say that there is coupling convergence in finite time (or, merely, coupling) of a sequence {X n } to a stationary sequence {Y • θ n } if

P (X n+l = Y • θ n+l , ∀l 0) n→+∞ -→ 1 .
It is easy to show that this notion of coupling convergence implies total variation convergence

(X n → Y in total variation if sup A∈F |P (X n ∈ A) -P (Y ∈ A)| n→+∞ -→ 0).
Definition 5.4. We say that there is strong coupling convergence in finite time (or, merely, strong coupling) of a sequence {X n } to a stationary sequence {Y • θ n } if:

ν = min n 0 | X n+l • θ -(n+l) = Y, ∀l 0 is a.s finite .
Remark Strong coupling implies coupling but the converse is not true.

Theorem 5.5 (Borovkov's renovating events). We consider a SRS {x(n)} defined by:

x(n + 1) = f (x(n), a(n)), n 0, x(0) = x 0 .
If the random process {x(n), n ∈ N} admits a stationary sequence of renovating events {A(n)} such that P (A(0)) > 0, then there exists a finite random variable Z such that:

Z • θ = f (Z, a(0)) ,
and x(n) converges with strong coupling in finite time to Z • θ n .

In the previous theorem, we have considered a SRS defined with a unique initial condition, x 0 .

In the rest of the paper, we will be interested in having results that hold uniformly over the initial conditions. We will then use the following generalization of Borovkov's theorem.

Theorem 5.6. We consider a subset V of E (V = E is in particular possible). We suppose that there exists a stationary sequence of events {A(n)} verifying P (A(0)) > 0 and which is renovating for the SRS {x(n, x 0 )}, ∀x 0 ∈ V . Then, for all (possibly random) initial condition x(0) such that P (x(0) ∈ V ) = 1, the sequence {x(n)} converges with strong coupling to a unique stationary regime.

Theorem 5.7 (converse of Th. 5.5 and 5.6). The conditions of Theorem 5.5 are necessary and sufficient for strong coupling convergence. Let V be a compact subset of E. The conditions of Theorem 5.6 are necessary and sufficient for strong coupling convergence uniformly over initial conditions in V .

Next theorem was proved by Anantharam and Konstantopoulos in [START_REF] Anantharam | Stationary solutions of stochastic recursions describing discrete event systems[END_REF].

Theorem 5.8. Let (Ω, F, P ) be a probability space. We assume that (Ω, F) is a Polish space equipped with its Borel σ-algebra. We consider a SRS "x(n + 1) = f (x(n), a(n))" defined on E. Suppose that, for some x 0 ∈ E, the sequence {x(n, x 0 )} is tight2 on E. Then there is a stationary distribution for the SRS.

The stationary distribution is defined on Ω × E with an Ω marginal equal to P . It provides only a weak stationary regime (wsr) for the SRS, see [START_REF] Anantharam | Stationary solutions of stochastic recursions describing discrete event systems[END_REF] or [START_REF] Brandt | Stationary Stochastic Models[END_REF] for details. All we need to know about wsr is that stationary regimes are wsr. Hence, the uniqueness of stationary regimes implies the uniqueness of wsr.

It is proved in [START_REF] Borovkov | Stochastically recursive sequences and their generalizations[END_REF], that for an i.i.d. SRS (i.e. Markov chain), the conditions of Th. 5.5 are equivalent to the ones ensuring Harris ergodicity. In Harris' framework, the conditions are on the state space. In Borovkov's framework, the conditions are on the exogenous driving sequence. This second approach is better suited for our problem. On the one hand, a direct analysis on the state space appears to be almost inextricable. On the other hand, the renovating events will take a very convenient form because a product of matrices is still a matrix (see Theorems 6.8, 6.10).

η-coupling Coupling and strong coupling, introduced above, are related to total variation convergence. We define now the notion of η-coupling. It is related to weak convergence.

Definition 5.9 (η-coupling). We consider a metric space (E, d). We consider two sequences {X n } n∈N and {Y n } n∈N defined on E. We say that there is η-coupling3 of these two sequences if for each η > 0, one can find versions of {X n } and {Y n } defined on a common probability space and an a.s. finite random time N such that

n N =⇒ d(X n , Y n ) η .
The following proposition is shown in Asmussen [START_REF] Asmussen | On coupling and weak convergence to stationarity[END_REF].

Proposition 5.10. We consider a sequence {X n } n∈N and a stationary sequence {Y • θ n } n∈N defined on the metric space E. Let µ be the invariant distribution of Y . If there is η-coupling of the two sequences, then {X n } converges weakly to µ.

Presentation of the Results

Let us consider a probability space (Ω, F, P, θ). The probability P is stationary and ergodic with respect to the shift θ. We are interested in systems of the type:

x(n + 1) = A(n) ⊗ x(n), n ∈ N x(0) = x 0
where x(n) and A(n) (∀n) are finite, respectively R k max and R k×k max -valued, random variables. We are sometimes going to use the notation x(n, x 0 ) to emphasize the value of the initial condition. We will consider models where the sequences {A(n), n ∈ N} are respectively i.i.d or stationary and ergodic (i.e A(n + 1) = A(n) • θ).

We recall that we have defined first and second order limits in Section 1, Equations ( 3) and (4). We are going to recall results on first and second order limits for open systems and first order limits for closed systems before completing the picture by solving the problem of second order limits for closed systems. [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF] and Cohen [START_REF] Cohen | Subadditivity, generalized product of random matrices and operations research[END_REF] For x ∈ R k and A ∈ R k×k max , we use the notation

Results from Baccelli

x ∞ = k i=1 x i and A ∞ = k i,j=1 A ij .
6.1.1 First order limits for closed systems Theorem 6.1 (Cohen [17]). Let {A(n)} be a stationary and ergodic sequence of matrices. We suppose that the matrix A(0) has a fixed structure (see Definition 3.5), is irreducible and

verifies P (A ij (0) = ε) = 1 or ε < E(A ij (0)) < +∞.
There exists a constant λ ∈ R such that, for all initial condition x 0 and for all i ∈ {1, . . . , k}:

lim n x i (n, x 0 ) n = lim n E x i (n, x 0 ) n = λ, P -a.s.
The constant λ is called the Lyapunov exponent of the stochastic matrix A(0).

The basic idea is to use the inequality A ⊗ B ∞ A ∞ ⊗ B ∞ in order to apply Kingman subadditive ergodic theorem. This definition of a Lyapunov exponent is coherent with the one of Theorem 4.2. Indeed, by Theorem 4.5, for every irreducible and deterministic matrix A, there exists d and M such that ∀m M, A m+d = λ d ⊗ A m , where λ is the eigenvalue of A. It implies that ∀x 0 ∈ R k max , lim n A n x 0 /n = λ.

First order limits for open systems

We suppose that matrices A(n) have a fixed structure. We decompose the graph of A(0) into its maximal strongly connected subgraphs (s.c.s.). If we replace each s.c.s. by one node, we obtain an associated reduced graph which is acyclic. We associate with each node ũ of the reduced graph a constant λ ũ which is the Lyapunov exponent of the corresponding s.c.s. in isolation, see Theorem 6.1. We denote by • ũ the set of predecessors of ũ (including ũ) in the reduced graph.

We have : Theorem 6.2 (Baccelli [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF]). Let {A(n)} be a stationary and ergodic sequence of matrices. We suppose that A(0) has a fixed structure. We suppose also that P (A ij (0) = ε) = 1 or ε < E(A ij (0)) < +∞, ∀i, j. Let us consider i ∈ {1, . . . , k}, i belongs to the s.c.s. ũ.

lim n x i (n, x 0 ) n = lim n E x i (n, x 0 ) n = ṽ∈ • ũ λ ṽ , P -a.s. .
Intuitively, the dynamic of the system is imposed by the s.c.s. having the smallest throughput (largest cycle time λ).

Second order limits for open systems

Matrices A(n) have a fixed structure. In order to simplify the presentation of the results, let us assume that the structure consists of two s.c.s. The general case is completely similar. Up to a permutation of the coordinates, we have :

A(n) = Ũ (n) ε B(n) Ã(n) .
The block Ũ is a square matrix of size I × I, irreducible. It is interpreted as the input of our system. The block à is a square matrix of size (k -I) × (k -I) , irreducible. The block B is the matrix of the communications between the sources ( Ũ ) and ( Ã). We suppose that the block Ũ in isolation has a unique stationary regime. We have the following theorem.

Theorem 6.3 (Baccelli [3]). Let u and a be the Lyapunov exponents of Ũ and à respectively (see Theorem 6.1). If a < u, there is a unique stationary regime for the SRS π(x(n)), regardless of the initial condition. Convergence to the stationary regime occurs with strong coupling. If a > u, then the differences of the form

x j (n, x 0 )x i (n, x 0 ), i = 1, . . . , I, j = I + 1, . . . , k , tend to +∞, Pa.s., for all finite initial condition.

If u > a, the sources which are slower impose their pace. If u < a, everything happens asymptotically as if à were in isolation.

Remark

In the previous theorem, we need the assumption that Ũ in isolation has a unique stationary regime. Knowing if Ũ has a unique stationary regime is precisely the problem which is going to be addressed in the following. Then, to determine if there is a unique stationary regime for π(x(n)), we have to use the results of Section 6.4 (applied to Ũ ) together with the comparison of Lyapunov exponents (of Ũ and Ã).

Remark In the results above (Theorems 6.1, 6.2 and 6.3), the assumption that a matrix M (0) has a fixed structure and is irreducible can be weakened and replaced by:

lim n P ( M (n) ⊗ M (n -1) ⊗ • • • ⊗ M (0) irred. ) = 1, ε < E(M ij (0) | M ij (0) = ε) < +∞ .

Preliminary results

From now on, we concentrate on second order limits in the closed (i.e A(n) is P -a.s. irreducible) case. The limits are expected to be random variables. We are interested in determining whether the limiting distribution is unique. Furthermore, we want to investigate the type of convergence to the limit.

We recall that π is the canonical projection R k π -→ PR k (Def. 3.6). It is clear that the recursive equation x(n + 1) = A(n)x(n) defines a SRS (Def. 5.1). It implies that π(x(n)) is also a SRS. Indeed, let us consider x(n) and x ′ (n) such that π(x(n)) = π(x ′ (n)). We define x(n + 1) = A(n)x(n) and x ′ (n + 1) = A(n)x ′ (n). It is straightforward that π(x(n + 1)) = π(x ′ (n + 1)). We write with some abuse of notation that {π(x(n)), n ∈ N} verifies the recursive equation "π(x(n + 1)) = πA(n)π(x(n))"4 . Lemma 6.4 (Reising & al [START_REF] Resing | Asymptotic behavior of random discrete event systems[END_REF]). For i ∈ {1, . . . , k}, we define

z i (n) = x i (n) -x i (n -1). We have z i (n) = F i ( A(n -1), π(x(n -1)) )
, where F i is an absolutely continuous function.

The sequence {A(n)} being stationary by hypothesis, it implies the following corollary. Corollary 6.5. A sufficient condition for (z 1 (n), . . . , z k (n)) ′ to converge weakly (resp. in total variation) to a unique invariant distribution, uniformly over initial conditions in PR k , is that π(x(n)) has the same property. This sufficient condition is not necessary as demonstrated by the following deterministic example.

Example 6.6. Let us consider

A = e -1 -1 e .

We have A + = A, so u 1 = (e, -1) ′ and u 2 = (-1, e) ′ are eigenvectors of A. The set

{u λ = λ ⊗ u 1 ⊕ (1 -λ) ⊗ u 2 , λ ∈ [0, 1]} ,
is the set of eigenvectors of A, see Theorem 4.4. There is a continuum of stationary regimes for π(x(n)). For example, it is easy to check that for an initial condition u λ , λ ∈ [0, 1], we have:

x 1 (n, u λ ) -x 2 (n, u λ ) = 2× λ -1 .
But on the other hand, we have a unique stationary regime for z i (n). As a direct consequence of the equality

A 2 = A, we have z 1 (n) = z 2 (n) = e, ∀n 2.
We can also easily build stochastic counter-examples of the same kind.

Remark The variables π(x(n)) depend only on the sequence {π(A(n))}. Therefore, all the results on π(x(n)) would still be true under the weaker assumption that only the sequence {π(A(n))} is stationary and ergodic. But, on the other hand, the variables z i (n) depend on the sequence {A(n)} and not only on {π(A(n))}. Corollary 6.5 would not be true under the assumption that {π(A(n))} is stationary and ergodic.

In the rest of the paper, we investigate the existence of a stationary regime for the SRS π(x(n)), i.e. the existence of a finite r.v.

Z : Ω → PR k such that 5 Z • θ = π A(0)π -1 (Z) .
We are interested by conditions ensuring the uniqueness of the stationary regime and the convergence of π(x(n, x 0 )) toward it, for all x 0 ∈ R k . In such cases, we say that the model is stable. Two types of convergence will appear, convergence with η-coupling and convergence with coupling. They imply, respectively, weak convergence and total variation convergence as recalled in §5.

Stability of discrete models

Let {A l , l ∈ L or l ∈ N}, be a finite or countable collection of irreducible matrices of size k × k. We suppose that there exists a discrete probability law {p l } such that A(n, ω) = A l with probability p l > 0.

Definition 6.7 (pattern, 1).

A matrix à is called a pattern of the random sequence {A(n), n ∈ N} if:

1. ∃N | Ã = A u N-1 ⊗ • • • ⊗ A u 0 with u 0 , . . . , u N -1 ∈ L (or N).

P (A(N

-1) ⊗ • • • ⊗ A(0) = Ã) > 0.
If the sequence {A(n)} is i.i.d. then the second condition is always verified.

Theorem 6.8. The sequence of matrices {A(n)} is i.i.d. If there exists a pattern of {A(n)} whose critical graph has a unique s.c.s. and whose cyclicity is 1 (scs1-cyc1 matrix), then {π(x(n))} converges with strong coupling to a unique stationary regime. It implies total variation convergence of {π(x(n))} to its stationary distribution.

Proof. Let C = A u N-1 ⊗ • • • ⊗ A u 0 be a scs1-cyc1 pattern. We have, using the cyclicity 1 assumption (Th. 4.5),

∃M | ∀m M, C m+1 = λC m ,
where λ is the Lyapunov exponent of C. We conclude that for all initial condition y,

C M +1 y = C(C M y) = λ ⊗ C M y. It means that C M y is an eigenvector of C.
By the assumption on the critical graph of C, there is a unique eigenvector (up to a constant) denoted y 0 (Th. 4.4). We have C M ⊗ y = µ(y) ⊗ y 0 , µ(y) ∈ R, or equivalently π(C M y) = π(y 0 ). We define

B i = {ω| A(i + M N -1, ω) ⊗ • • • ⊗ A(i + 1, ω) ⊗ A(i, ω) = C M } .
From the i.i.d. assumption, it follows that P(B i ) > 0. On B i , and for all initial condition y, we have:

x(i + M N ) = C M ⊗ x(i) = µ(x(i)) ⊗ y 0 =⇒ π(x(i + M N )) = π(y 0 ) .
We check that the sequence B i is compatible with the shift, i.e. B i = B 0 • θ i . We conclude that B i is a stationary renovating event sequence for the SRS π(x(n)). We apply Borovkov's Theorem (version 5.6 for the set V = PR k , as we have obtained a sequence of renovating events independent of the initial condition) and the uniqueness of the stationary regime follows. ⊓ ⊔ Example 6.9. [Cyclic Jackson Network 2] We consider a basic Cyclic Jackson Network with k queues and k customers. Such a network can be represented by the (max,+) matrix given in §2.2, Equation ( 5). We suppose that the sequence of service times {(σ 1 (n), . . . , σ k (n)) , n ∈ N} is i.i.d. However the random variables σ 1 (n), . . . , σ k (n) need not be independent for a given n.

We suppose also that the service times have a discrete support, i.e. can only take a countable number of values. We are in the framework of Theorem 6.8. We conclude that a sufficient condition of stability is to find a scs1-cyc1 matrix among the (max,+) matrices corresponding to this network. As a direct application of the result stated in Example 4.8, we obtain that a condition of stability is:

P (∃i | σ i (n) > σ j (n), ∀j = i) > 0 or P (σ 1 (n) = σ 2 (n) = • • • = σ k (n)) > 0 .
We now give a version of Theorem 6.8 in the stationary and ergodic case.

Theorem 6.10. The sequence {A(n)} is stationary and ergodic. We suppose that there exists a finite pattern

C = A u N-1 ⊗ • • • ⊗ A u 0 which is scs1-cyc1 and of rank 1 (see Def. 4.7). We suppose that B = {ω | A(N -1)A(N -2) • • • A(1)A(0) = C} is of strictly positive probability.
Then {π(x(n))} converges with strong coupling to a unique stationary regime.

Proof. The proof resembles the one of Theorem 6.8. As C is of rank 1, we have (see Def. 4.7):

C 2 = λ ⊗ C
, where λ is the Lyapunov exponent of C. We conclude that:

∀y ∈ R k , C 2 y = C(Cy) = λCy .
It implies that Cy is an eigenvector of C. As matrix C is scs1, it has a unique eigenvector y 0 , up to a constant. On

B i = B • θ i , we have π(x(i + N )) = π (Cx(i))
= π(y 0 ) .

We check that the sequence B i is compatible with the shift and we apply Borovkov's Theorem 5.6.

⊓ ⊔

Remark If the dependence between matrices is markovian, a sufficient condition to get P (B) > 0 is that p(A u i , A u i+1 ) > 0, ∀i = 1, . . . , N -1, where p(., .) is the markovian transition kernel.

Remark The conditions of this theorem are, of course, weaker than the i.i.d. assumption of Theorem 6.8. However we made an additional assumption, namely that the pattern C is of rank 1. This assumption cannot be relaxed, as shown by the counter-example 6.11.

Example 6.11. Let Ω = {ω 1 , ω 2 } be the probability space, P={ 1 2 , 1 2 } the probability law, and θ the stationary and ergodic shift defined by: θ(ω 1 ) = ω 2 and θ(ω 2 ) = ω 1 . We consider

A = 1 -η e e 1 ,
n > [1/η], P (∃N | A(N - 1) • • • A(0) = A n ) = P (∃N | A(N -1) • • • A(0) = B n ) = 0.
Hence the conditions of Theorem 6.10 are not verified. In fact, there is a continuum of possible periodic limits. Consider

x 0 = (a, b) ′ with -1 + η < a -b < 1 -η. Then the limit regime of π(x(n)) has a state space which is either {π(a, b) ′ , π(a + η, b) ′ } (with probability 1 2 ), or {π(a, b) ′ , π(a, b + η) ′ } (with probability 1 2 ).

Stability of general models

In this section, we consider a general model where the coordinates of our matrices have a support which can be discrete, absolutely continuous with respect to Lebesgue measure or a mixture of these two cases.

We need the following definition, extending the notion of pattern we have been using for finite models. Let M be a deterministic matrix and η > 0. We denote by B(M, η) the open ball of center M and of radius η for the supremum norm of R k×k . We have

N ∈ B(M, η) iff ∀i, j, N ij ∈]M ij -η, M ij + η[ Definition 6
.12 (pattern, 2). Let A be a random matrix. We say that à is a pattern of A if à is a deterministic matrix verifying ∀η > 0, P A ∈ B( Ã, η) > 0 .

Equivalently, we can say that à belongs to the support of the random matrix A. It includes the cases where à is an accumulation point (discrete case) or a boundary point (continuous case) of the support. Definition 6.13 (pattern, 3). Let {A(n), n ∈ N} be a sequence of random matrices. We say that the deterministic matrix à is a pattern of the sequence

{A(n)} if ∃N s.t. ∀η > 0, P A(N -1) ⊗ • • • ⊗ A(0) ∈ B( Ã, η) > 0 .
Equivalently, we can say that à is a pattern (Def. 6.12) of the random matrix

A(N -1) ⊗ • • • ⊗ A(0). We say that à is an asymptotic pattern of {A(n)} if ∀η > 0, ∃N η s.t. P π(A(N η -1) ⊗ • • • ⊗ A(0)) ∈ π(B( Ã, η)) > 0 .
Remark This definition is coherent with the one given in Definition 6.7 for a discrete model. Note that, for convenience reasons, asymptotic patterns are defined in the projective space PR k×k .

Theorem 6.14. The matrices A(n) are i.i.d. (resp. stationary and ergodic). We suppose that there exists a matrix C which is a pattern of {A(n)} (see Def. 6.13) and which is scs1-cyc1 (resp. of rank 1). Then the SRS {π(x(n))} has a unique stationary regime {Z • θ n }.

The convergence occurs with η-coupling. It implies weak convergence of π(x(n)) to its unique stationary distribution.

Proof. We prove directly Theorem 7.4, a stronger version of the result. It is done in Appendix, §9.3. ⊓ ⊔ Theorem 6.15. The sequence of matrices {A(n)} is i.i.d. or stationary and ergodic. We assume that there exists a set C of matrices such that :

1. ∀C ∈ C, C is a matrix of rank 1.

2. ∀C ∈ C, C is a pattern of {A(n)}.

∃N |

P (A(N -1) • • • A(0) ∈ C) > 0.
Theorem 7.1. The sequence of matrices {A(n)} is i.i.d. or stationary and ergodic. When there is a unique stationary regime, convergence to this regime occurs with strong coupling. A necessary and sufficient condition for the model to have a unique stationary regime is that there exists a matrix C verifying 1. C is a matrix of rank 1 (Def. 4.7).

2. C is a pattern of {A(n)} (Def. 6.13).

Proof. It is given in Appendix, §9.2. ⊓ ⊔ Theorem 7.1 is not true in general when the matrices A l , l ∈ L, belong to R k×k max , see the following counter-example.

Example 7.2. We consider the matrices

A = e -1 -1 -η , B = -η ′ -1 -1 e ,
where 0 < η, η ′ ≪ 1 and η, η ′ are not co-rational, i.e. η/η ′ ∈ Q.

Let u = (u 1 , u 2 ) ′ ∈ R 2 , we set ψ(u) = u 2 -u 1 .
We identify PR 2 and R using the function ψ •π -1 . The matrices A and B are scs1-cyc1. Their respective and unique eigenvectors are ψ(e 1 ) = -1 and ψ(e 2 ) = 1. For a vector u = (u 1 , u 2 ) ′ such that ψ(u) ∈ [-1, 1], we have

ψ(Au) = max(ψ(u) -η, -1), ψ(Bu) = min(ψ(u) + η ′ , 1) . (8)
We consider a Markov chain defined on the set ψ -1 [-1, 1] ⊂ R 2 . The transition probabilities are

• For u such that ψ(u) ∈] -1 + η, 1 -η ′ [, p(u, Au) = 1/2, p(u, Bu) = 1/2.
• For u such that ψ(u) ∈ [1η ′ , 1], p(u, Au) = 1.

• For u such that ψ(u) ∈ [-1, -1 + η], p(u, Bu) = 1.

The behaviour of the Markov chain is illustrated in Figure 2.

-1 1

B B A -1 + η 1 -η ′ PR 2 ≃ R A Figure 2: Markov chain ψ(X(n)) on R.
Let X(n) be a realization of the Markov chain. It is easy to check that this Markov chain is aperiodic. Under the assumption η/η ′ ∈ Q, one can prove using classical arguments that the set {ψ(X(n)), n ∈ N} is P -a.s. dense in [-1, 1]. It implies that the Markov chain is ν-irreducible where ν is the Lebesgue measure on ψ -1 [-1, 1]. Hence there exists a unique stationary distribution Q for the Markov chain. It verifies Q(A) > 0 for all event A such that ν(A) > 0. For a complete treatment of Markov chains on continuous state spaces, see Meyn & Tweedie [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF].

Let us consider a stationary realization X(n) of the Markov chain (i.e. ∀n, P {X(n) ∈ .} = Q(.)).

We define

A(n, ω) = A if X(n + 1, ω) = AX(n, ω), B if X(n + 1, ω) = BX(n, ω).
As X(n) is stationary, it follows that {A(n)} is a stationary and ergodic sequence.

Let us consider the stationary-ergodic finite model "x(n + 1) = A(n)x(n)" and x(0) = x 0 ∈ R 2 . Note that {x(n)} is not a Markov chain anymore.

Let us consider a pattern

C = A n-1 • • • A 0 of {A(n)}, i.e. P {A(n-1) • • • A(0) = A n-1 • • • A 0 } > 0.
Let x 0 , . . . , x n be a corresponding path for the Markov chain X(n), i.e.

x 0 ∈ ψ -1 ] -1, 1[, x 1 = A 0 x 0 , . . . , x n = A n-1 x n-1 ,
and

P {X(n) = x n , . . . , X(1) = x 1 | X(0) = x 0 } > 0 .
Let us denote by c and c the minimal distances between x p , p n and the extremal points of ψ -1 [-1, 1]. We have c = min p n (ψ(x p ) + 1) and c = min p n (1ψ(x p )). It follows from ( 8) that

ψ(Im(A n-1 • • • A 0 )) = [x n -c, x n + c] , (9)
where Im(A) = {Au, u ∈ R k }. From the definition of the Markov chain X(n), it follows that c > 0, c > 0. We conclude that A n-1 • • • A 0 is not a rank 1 matrix. There exists no finite pattern of rank 1 for {A(n)}.

On the other hand, let us prove that there exists asymptotic patterns of rank 1 for {A(n)}.

We define c(n) = min p n (ψ(X(p)) + 1) and c(n) = min p n (1ψ(X(p))). As {ψ(X(n)), n ∈ N} is dense in [-1, 1], we obtain that c(n) → 0 and c(n) → 0. Using [START_REF] Borovkov | Limit theorems for queueing networks[END_REF], we obtain that D(A(n) • • • A(0)) → 0, P -a.s. We conclude following the lines of Theorem 7.4, §9.3. There is a unique stationary regime for the model. For an arbitrary initial condition, we have η-coupling (weak convergence) with this stationary regime.

To summarize, we have exhibited a finite model with a unique stationary regime and no coupling convergence. This type of behaviour is closely related to the non-finiteness of the projective semigroup π < A, B >, see Def. 9.2.

General models

We consider a general model of type "x(n+1) = A(n)x(n)". There exists a set C of matrices such that :

1. ∀C ∈ C, C is a matrix of rank 1.

2. ∀C ∈ C, C is a pattern of {A(n)}.

∃N |

P (A(N -1) • • • A(0) ∈ C) > 0.
We can say equivalently that we must have patterns of rank 1 but with strictly positive probability.

Proof. We have already proved the sufficient part (Th. 6.15). We prove the necessary part of the theorem in Appendix, §9.4.

⊓ ⊔

Convergence with η-coupling appears as a limiting case of coupling in finite time. In a discrete model, we will have only η-coupling when the set C of scs1-cyc1 patterns is non empty but is of probability 0. It means that the scs1-cyc1 patterns are only accumulation points of the support. In a general model, we will have only η-coupling when the scs1-cyc1 patterns are isolated points of the support (which implies that they are boundary points of the support).

Example 7.6. To illustrate the previous remark, let us continue the analysis of Example 7.3.

There is only one matrix (in the projective space PR k×k max ) verifying the first two conditions of Theorem 7.5. It is the matrix

π(C) = π
e e e e , But condition 3. of Th. 7.5 is not verified as ∀N, P (π(A(N -1)

• • • A(0)) = π(C)) = 0.
Let us consider a slightly modified sequence of matrices { Ã(n)} where the diagonal elements are two random variables U (n) and U ′ (n) defined on [0, 1] and such that

P U (n) > U ′ (n) > 0 or P U ′ (n) > U (n) > 0 .
Now, we have scs1-cyc1 patterns with strictly positive probability and there is coupling in finite time with the unique stationary regime.

i (n) > σ j (n), ∀j = i} ∪ {σ 1 (n) = • • • = σ k (n)}) > 0 .
is necessary and sufficient for strong coupling convergence to a unique stationary regime. For i.i.d. Cyclic Jackson Networks, the sufficient condition P ({∃i | σ i (n) > σ j (n), ∀j = i}) > 0 was obtained in [START_REF] Kaspi | Regenerative closed queueing networks[END_REF]. The method of proof was completely different, see the remarks at the end of Section 2.2.

Without Irreducibility

We have supposed from the beginning that the matrices {A(n)} were irreducibles. The relaxation of the irreducibility assumption is very important in terms of modeling power. It enables us to consider, for example, the task graphs with random precedences introduced in Section 2.1.

The irreducibility assumption is used in Prop. 3.8. But the only point we need to prove this proposition is: "∀u ∈ R k , Au ∈ R k ", i.e. if u has only non-ε coordinates then Au has the same property. So the only assumption we need on the matrices {A(n)} is:

I ∀i, P {∃j s.t. A ij (0) > ε} = 1 .
The irreducibility is also essential for the R max spectral theory of Section 4. A reducible matrix A ∈ R k×k may have several eigenvalues. Definition 4.3 and Theorem 4.4 have to be reinterpreted by replacing the unique eigenvalue by the maximal eigenvalue. Theorem 4.5 is not true anymore. But as far as the direct theorems (6.8, 6.10, 6.14 and 6.15) are concerned, we use results of the R max spectral theory only for the pattern C whose existence is critical for the proofs. These theorems are still valid if we state that condition I is verified and that there exists a pattern C which is scs1-cyc1 and irreducible.

Dropping the irreducibility assumption does not influence the converse results. More precisely, the proofs of Theorems 7.1, 7.4 and 7.5 are still valid. Only two conditions need to be verified:

I ∀i, P {∃j s.t. A ij (0) > ε} = 1 . II ∃n | P {D( A(n) • • • A(0) ) < +∞} > 0.
Irreducibility Pa.s. is not necessary to ensure that these conditions hold. We conclude that we can state our converse results under the previous two minimal assumptions.

A counter-example shows that without condition II, the uniqueness of the stationary regime does not imply the existence of a rank 1 pattern.

Example 8.1. Let Ω = {ω 1 , ω 2 } be the probability space, P = { 1 2 , 1 2 } the probability law, and θ the ergodic shift defined by: θ(ω 1 ) = ω 2 and θ(ω 2 ) = ω 1 . We consider

It is immediate that Π is included in the convex hull of these vectors, i.e.

Π ⊂ {π(α 1 ⊗ c 1 ⊕ α 2 ⊗ c 2 • • • ⊕ α k ⊗ c k ), α i ∈ R} .
In the (max,+) algebra, we have the following property, for all

A ∈ R k×k max , u, v ∈ R k max , A(u⊕v) = Au ⊕ Av. It implies ∀x ∈ Π, π(Ax) ∈ {π(α 1 ⊗ Ac 1 ⊕ • • • ⊕ α k ⊗ Ac k ), α i ∈ R} . (11)
We fix η > 0. Using the η-coupling assumption, we have that the random variable N ′ (ω) is Pa.s. finite, where N ′ is defined by:

N ′ = inf{n | d(x(n, c i ), x(n, c j )) η, ∀i, j} .
As both N and N ′ are Pa.s. finite, we have

∀δ > 0, ∃L, L ′ : P {N L, N ′ L ′ } 1 -δ .
As a direct consequence of (11), we have on the event {N L, N ′ L ′ } :

D A(L ′ -1) ⊗ • • • ⊗ A(0)A(-1) ⊗ • • • ⊗ A(-L) η .
We deduce, using the stationarity of {A(n)}, that

P D(A(-1) • • • A(-L -L ′ )) η = P D(A(L ′ -1) • • • A(-L)) η 1 -δ .
It implies that the random variables D(A(-1) • • • A(-n)) converge in probability to 0. But as

D(A(-1) • • • A(-n))
is pathwise decreasing, the convergence occurs also Pa.s.

We have in particular, for all u

∈ R k , d ( A(-1) • • • A(-n)u, A(-1) • • • A(-n -p)u) -→ 0, P - a.s. It implies that {π(A(-1) • • • A(-n)u
)} is a Cauchy sequence which converges. The limit does not depend on u. We denote it by Z. We have :

Z • θ = lim n π(A(0)A(-1) • • • A(-n)u) = πA(0) lim n π(A(-1) • • • A(-n)u) = π(A(0))Z .
The sequence {Z • θ n } is a stationary regime. Let us prove it is the unique one. We want to prove that

∀x 0 ∈ R k , d (x(n, x 0 ), Z • θ n ) n→+∞ -→ 0 , P -a.s. ( 12
)
As Z is P -a.s. finite, for all η > 0, there exists a compact K ∈ PR k such that P {Z ∈ K} > 1-η. We proceed as above (Equation ( 10)) in order to define vectors c i , i = 1, . . . , k such that

K ⊂ {π(α 1 ⊗ c 1 ⊕ α 2 ⊗ c 2 • • • ⊕ α k ⊗ c k ), α i ∈ R} .
We have

Z ∈ K ⇒ Z • θ p ∈ {π(α 1 [A(p -1) • • • A(0)c 1 ] ⊕ • • • ⊕ α k [A(p -1) • • • A(0)c k ]), α i ∈ R} .
Using the η-coupling of trajectories, we also have ∀c i , d(x(n, x 0 ), x(n, c i )) → 0 .

We conclude easily that there is η-coupling of {πx(n, x 0 )} and {Z • θ n }. We can apply Proposition 5.10. There is weak convergence of {πx(n, x 0 )} to the distribution of Z and relation [START_REF] Bougerol | Products of Random Matrices with Applications to Schrödinger Operators[END_REF] establishes the a.s. convergence of {πx(n, x 0 )} to Z • θ n . As a direct consequence, {Z • θ n } is the unique stationary regime. ⊓ ⊔ 9.2 Proof of Theorem 7.1

We are going to prove that the existence of a unique stationary regime implies the existence of a pattern of rank 1 (Def. 6.7). Using Theorem 6.10, the proof will then be complete. We need the following definition.

Definition 9.2. Let us consider A 1 , . . . , A p ∈ R k×k max . We denote by < A 1 , . . . A p >, the semigroup generated by these matrices and by π < A 1 , . . . A p > its projection. We have Let us prove a lemma first.

It implies that the sequence {π(B(n))} belongs to a compact of (PR k×k , |.| P ). Hence there exists a strictly increasing function σ : N → N such that π(B(σ(n))) is converging. Let A ∞ be a representative (in R k×k ) of the limit. By continuity of the projective distance, we have that ∀u, v ∈ R k , d(A ∞ u, A ∞ v) = 0. Therefore A ∞ is a rank 1 matrix.

As the products {π(A(n) • • • A(0))} can only take a finite number of values in compact sets (Proposition 9.3), it implies that the limit matrix A ∞ is attained in finite time. More precisely, there exists N such that ∀n N, π (B(σ(n))) = π(A ∞ ) .

The matrix B(σ(N )) is a rank 1 pattern for {A(n)}. It concludes the proof.

Proof of Theorem 7.4

We first prove the necessary part of the Theorem, i.e. η-coupling with a unique stationary regime implies the existence of an asymptotic pattern. Let Z • θ n be the unique stationary regime. We have for all x 0 , y 0 ∈ R k , d(x(n, x 0 ), Z • θ n ) → 0, d(x(n, y 0 ), Z • θ n ) → 0 ⇒ d(x(n, x 0 ), x(n, y 0 )) → 0 .

We have assumed that ∃N such that P {D(A(N -1) • • • A(0)) < +∞} > 0, see Equation ( 7

As the matrices {B p , p ∈ N} belong to a compact, there exists a subsequence {B σ(p) } which converges to a limit B ∞ . We have (see the proof of Lemma 9.1) that D(A(-1) • • • A(-n)) → 0, Pa.s.. We conclude that B ∞ is a rank 1 matrix.

Proof. Let {Z • θ n } be the unique stationary regime with which the SRS π(x(n)) couples. We consider the event:

Y n = ω | π(x(n, x 0 )) • θ -n ω = Zω, ∀x 0 ∈ R k .
The assumption of coupling in finite time, uniformly over R k , may be written :

P (Y n ) n→+∞ -→ 1 .
Here we implicitly use the assumption that the projective image of A(-1) • • • A(-n) is asymptotically bounded (see Equation ( 7)). Let us consider ω ∈ Y n and p an integer > 0, we have:

π(x(n + p, x 0 )) • θ -n-p ω = π x(n, x(p, x 0 ) • θ -p ) • θ -n ω (20) = Zω (as ω ∈ Y n ) . (21)
The passage from ( 20) to ( 21) uses the fact that coupling occurs uniformly over initial conditions. We have:

Y n = ω | π(x(n + p, x 0 )) • θ -(n+p) ω = Zω, ∀p 0, ∀x 0 ∈ R k , and
P (Y n ) n→+∞ -→ 1 .
This is exactly the definition of strong coupling (Def. 5.4).

⊓ ⊔

We can now use the converse Theorem 5.7. There exists a stationary sequence of events {A • θ n } which is renovating for the SRS {π ( x(n, x 0 ) )}, ∀x 0 ∈ R k , and verifies P (A) > 0. Let m be the common length and Φ the common function of these renovating events. We have, on A: π(x(m)) = Φ (A(m -1), . . . , A(0)) , ∀x 0 ∈ R k .

But we also have:

x(m) = A(m -1) ⊗ • • • ⊗ A(0) ⊗ x 0 , ∀x 0 ∈ R k .
We conclude that, on A, π (A(m -1, ω) ⊗ • • • ⊗ A(0, ω) ⊗ x 0 ) is independent of x 0 . It implies that C = A(m -1, ω) ⊗ • • • ⊗ A(0, ω) is a matrix of rank 1.

B = 1 e e 1

 1 η , e < η < < 1 . {A(n, ω 1 )} = A, B, A, B, . . . {A(n, ω 2 )} = B, A, B, A, . . . . Both matrices A and B are scs1-cyc1 patterns of length 1. But patterns which are scs1-cyc1 and of rank 1 are for example A n or B n for n > [1/η]. We have for any

  ω 1 )} = A, B, A, B, . . . {A(n, ω 2 )} = B, A, B, A, . . . . is invariant by the translation shift. Then by the ergodic Lemma, it is of probability 0 or 1. We have made the assumption that P (A) > 0, hence P (A) = 1. Using the stationarity of the sequence {A(n), n ∈ Z}, we have that ∃N (ω) such thatD ( A(-1) • • • A(-N ) ) < +∞ .Then we can define the projective image of A(-1) • • • A(-N ) which is a bounded subset of PR k and that we denote by Π. The boundedness implies that c = max v∈Π d(e, v) < +∞ , where e = (e, . . . , e) ′ . Let us define the vectors c 1 = (c, e, . . . , e) ′ , c 2 = (e, c, e, . . . , e) ′ , . . . , c k = (e, . . . , e, c) ′ .

< A 1

 1 , . . . , A p > = {(A un • • • A u 2 A u 1 ) , u 1 , . . . , u n ∈ {1, . . . , p}, n ∈ N} , π< A 1 , . . . , A p > = {π(A un • • • A u 2 A u 1 ) , u 1 , . . . , u n ∈ {1, . . . , p}, n ∈ N} .We consider the Euclidean space (PR k×k , |.| P ) as introduced in Definition 3.7. Next proposition was proved by Gaubert in[START_REF] Gaubert | On semigroups of matrices in the (max, +) algebra[END_REF].Proposition 9.3. Let A 1 , . . . , A p ∈ Q k×k max .For all compact set K of (PR k×k , |.| P ), we have π < A 1 , . . . , A p > ∩K is finite.

  ), Section §7. Let K ∈ R be such that P {D(A(-1) • • • A(-N )) < K} > 0. It implies that there existsK ′ such that P {|A(-1) • • • A(-N )| P < K ′ } > 0. Let us denote E 0 = {ω | |A(-1) • • • A(-N )| P < K ′ } .It follows from the stationary-ergodic assumption, that there exists a minimal n 1 > 1 such thatP {E 1 } > 0, E 1 = E 0 ∩ {|A(-n 1 ) • • • A(-n 1 -N + 1)| P < K ′ } .We define in the same way an increasing sequence n p > • • • > n 2 > n 1 and a decreasing sequence of eventsE p ⊂ • • • ⊂ E 2 ⊂ E 1 verifying P {E p } > 0, E p = E p-1 ∩ {|A(-n p ) • • • A(-n p -N + 1)| P < K ′ } .On the event E p , p 1, we have|A(-1) • • • A(-n p -N + 1)| P < 2 × K ′ .The proof is exactly similar to the one proposed in the proof of Theorem 7.1 ( §9.2, Equation (15) and after). Let B(E, K ′ ) denote the open ball of (PR k×k , |.| P ) of center π(E), E ij = e, ∀i, j and of radius K ′ . For all p, we choose a deterministic matrix B p belonging to B(E, K ′ ) and verifyingP {E p } ∩ {πA(-1) • • • A(-n p -N + 1) ∈ B(B p , 1 p )} > 0 .

  i.d. (resp. stationary and ergodic) sequences of service times with general distributions, i.e. to replace ./M/1/∞ servers by ./GI/1/∞ (resp. ./G/1/∞) servers. Finding the minimal assumptions leading to a unique stationary regime for this generalized closed Jackson Network is still an open problem.

  Stability no longer implies coupling in finite time. It was illustrated by Example 7.2. Here is another example, for an i.i.d. model.where U (n) are i.i.d. random variables of uniform distribution over [0, 1]. There is a unique stationary regime for π(x(n)) which is π(e, e) ′ . We denote by d(., .) the projective distance. For an initial condition (y, e) ′ with y 1, we have d (x(n); (e, e) ′ ) = min p n U (p). Thus convergence to π(e, e) ′ occurs only asymptotically. There is no coupling but only η-coupling with the unique stationary regime. The sequence of matrices {A(n)} is i.i.d. or stationary and ergodic. The necessary and sufficient condition for the model to converge with η-coupling to a unique stationary regime is the existence of an asymptotic pattern C of {A(n)} of rank 1 (Def. 6.13). The sequence of matrices {A(n)} is i.i.d. or stationary and ergodic. The necessary and sufficient conditions for the model to converge with coupling to a unique stationary regime are :

	We can show the following results.
	Theorem 7.4. Proof. It is given in Appendix, §9.3.			⊓ ⊔
	Theorem 7.5.
	Example 7.3.
	A(n) =	U (n) e	e U (n)	,

  Example 7.7. [Cyclic Jackson Network 4] We consider the model of Example 6.16. The condition P ({∃i | σ

It is worth mentioning that d(., .) is the Rmax analogue of a distance used in classical Perron-Frobenius theory, which is called the Hilbert's projective metric and is defined by "δ(u, v) = ln (inf{µ/λ | λu v µu})".

Tightness on E means that for any η > 0, there is a compact K of E such that P {x(n, x0) ∈ K} > 1 -η, for all n.

The classical terminology is ε-coupling. We change it to η-coupling to avoid confusions with the notation ε = -∞ of the Rmax algebra.

It would be more rigorous to use different notations π and π for the canonical projections in PR k and PR k×k max respectively. Then we would define more formally π(A(n))π(x(n)) ≡ π(A(n)x(n)).

We will write Z • θ = πA(0)Z with some abuse of notations.

Acknowledgment I would like to thank François Baccelli who introduced me to this problem. F. Baccelli gave me also many ideas and suggestions which appear in this paper. I am also grateful to Serguei Foss, Stéphane Gaubert and Philippe Bougerol for several fruitful talks on the topic. At last, the careful comments of an anonymous referee have greatly helped improving the presentation of this paper.

Then {π(x(n))} converges with strong coupling to a unique stationary regime.

The conditions of Theorem 6.15 are stronger than the ones of Theorem 6.14 as we require the patterns of rank 1 to be of positive probability. On the other hand, we obtain a stronger type of convergence.

Proof. Let us define B = {ω | A(N -1)A(N -2) • • • A(1)A(0) ∈ C} and B i = B • θ i . Using that the matrices C ∈ C are of rank 1, we obtain that, on the event B i , π(x(i + N )) is independent of the value of π(x(i)). It implies that {B n , n ∈ N} is a stationary sequence of renovating events. The result follows.

⊓ ⊔ Remark Theorems 6.8 to 6.15 do not require any aperiodicity (Def. 3.4) assumption on matrices A(n). However, the pattern C whose existence is essential in all of these theorems is aperiodic. The condition "scs1-cyc1" implies aperiodicity.

Example 6. [START_REF] Cohen | A linear system-theoretic view of discreteevent processes and its use for performance evaluation in manufacturing[END_REF]. [Cyclic Jackson Network 3] We consider the same i.i.d. model as in Example 6.9. However, the distributions of the service times are now general. We obtain, by using Theorems 6.14 and 6.15, the stability under the condition:

The support of the random vector (σ 1 (n), . . . , σ k (n)) contains at least one point such that:

If the previous condition occurs with strictly positive probability, we obtain total variation convergence. Otherwise, we obtain weak convergence. Here is a case with only weak convergence. We consider an i.i.d. CJN with three queues and three customers. We assume that σ 1 = σ 2 = 1 and σ 3 is uniformly distributed over [0, 1].

Converse Theorems

We are going to prove converses of Theorems 6.8, 6.10, 6.14 and 6.15. We will consider successively finite and general models of type: "x(n + 1) = A(n)x(n)" where the matrices are of size k × k and are Pa.s irreducible. We will, moreover, always suppose that there exists a pattern whose projective diameter (Def. 3.9) is finite, i.e.

It implies lim n P {D( A(n) • • • A(0) ) < +∞} = 1, see the proof of Lemma 9.1. This condition is very weak. In the i.i.d. case, it is enough that there exists a pattern which is irreducible and aperiodic. We comment further on this condition in Section 8.

Finite models in Q k×k max

We consider a finite model: "x(n + 1) = A(n)x(n)", with A(n) ∈ {A l , l ∈ L = {1, . . . , L}}. We assume that the matrices are irreducible. We assume also that the matrices A l , l ∈ L, belong to Q k×k max , i.e. that their coordinates are rational.

All patterns have an infinite projective diameter. Therefore, condition II is not verified. Nevertheless, there is coupling in finite time with a unique periodic regime. More precisely, there is coupling of π(x(n, u)) to the periodic regime {π(e, e) ′ , π(e, 1) ′ } and coupling occurs for n > 2× [ (u 1u 2 , u 2u 1 )]. We conclude that there is coupling in finite time to a unique stationary regime but no rank 1 pattern. Without condition II, Theorem 7.1 is not true anymore.

Another class of systems where condition II is not verified is the class of open systems studied by Baccelli in [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF]. The results for this type of systems have been recalled in §6.1. In this case also, Theorem 7.1 fails to be true. In such models, there are no patterns which are scs1-cyc1 and irreducible (matrices are non-irreducible with probability 1 !). The good criterion to decide the uniqueness of the stationary regime is the comparison between Lyapunov exponents, see Theorem 6.3. The computation of such exponents involves the whole structure of the stochastic matrices {A(n)}, and not only an extracted pattern.

Condition II is weak and will be verified in most cases. For a discrete i.i.d. model, for example, it is sufficient to have one pattern of finite length Remark For a general model which does not verify condition II, we decompose the model into its maximal sub-models verifying it. Then the complete analysis of the system boils down to an analysis of the sub-models (using the results of Section 6.4) and of their interactions (using Theorem 6.3 and its generalizations, see [START_REF] Baccelli | Ergodic theory of stochastic Petri networks[END_REF]).

9 Appendix

Loynes scheme

Lemma 9.1 is going to be used in several of the forthcoming proofs. Under an assumption of η-coupling of the trajectories, we build a stationary regime using a Loynes' type construction.

Lemma 9.1. We consider a general model "x(n + 1) = A(n)x(n)" (see §6.4). The sequence {A(n)} is stationary and ergodic. We assume that there exists N such that

We assume also that ∀x 0 , y 0 ∈ R k , d(x(n, x 0 ), x(n, y 0 )) → 0, P a.s. (η-coupling of the trajectories). Then there exists a r.v.

is the unique stationary regime of the model.

Proof. We are going to show that the sequence {π(A(-1)

The argument is an analog of the famous backward scheme proposed by Loynes in [START_REF] Loynes | The stability of a queue with non-independent interarrival and service times[END_REF] for G/G/1 queues.

We want to show that

It is easy to see that the event

and A i ∈ Q k×k max , i = 1, . . . , p. We suppose that there is a unique stationary regime. It implies

Equivalently, it implies η-coupling of the trajectories corresponding to different initial conditions.

Proof. We assume that Equation ( 13) is not verified. It implies, using Proposition 3.8, that there exists x 0 , y 0 ∈ R k and c > 0 such that

Let S =< A 1 , . . . , A p > be the projective semigroup generated by the matrices of the model (Def. 9.2). For x ∈ R k , we define S(x) = {π(Ax), A ∈ S}. We have that π < A 1 , . . . , A p > ∩K is finite for all compact K of PR k×k , Proposition 9.3. It implies that S(x) ∩ K is finite for all compact K of PR k . We conclude that S(x) has no accumulation point and verifies S(x) = S(x),

where S(x) is the closure of S(x) in PR k .

We want to apply Theorem 5.8. It is required that the probability space be a Polish space. In order to fulfill this, we consider the canonical probability space consisting of one-sided infinite sequences of matrices {A 1 , . . . A p }, i.e. Ω = {(A u 0 , A u 1 , . . . , A un , . . . ), u i ∈ {1, . . . , p}} .

We recall that we made the assumption [START_REF] Borovkov | Asymptotic Methods in Queueing Theory[END_REF], which implies

It implies that for all η > 0, there exists N ∈ N and K, a compact set of PR k×k such that ∀n N, P {π(A(N )

There exists a compact K ′ (which depends on x) of PR k such that

We conclude that the sequence {π(x(n, x))} is tight in PR k . It implies that it is tight in S(x) = S(x). We can view π(x(n, x)) as a SRS defined on S(x) only. Applying Theorem 5.8, we obtain that, for all x, there exists a stationary distribution Q x defined on Ω × S(x).

Let us consider the initial conditions x 0 and y 0 as defined in [START_REF] Brilman | Synchronisation by resources sharing : a performance analysis[END_REF]. It is a-priori possible to have S(x 0 ) ∩ S(y 0 ) = ∅. As a consequence, one cannot rule out that Q x 0 = Q y 0 . We are going to prove that there exists α ∈ R such that S(x 0 ) ∩ S(αx 0 ⊕ y 0 ) = ∅. It will provide two stationary distributions Q x 0 = Q αx 0 ⊕y 0 , which contradicts the uniqueness of the stationary regime.

We work on the event A, see [START_REF] Brilman | Synchronisation by resources sharing : a performance analysis[END_REF]. We have d(x(n, x 0 ), x(n, y 0 )) > c for all n. Let x, y ∈ R k be two different points. Then there exists an open interval ]λ, λ[ such that

The proof is straightforward. We consider the (random) intervals ]λ(n), λ(n)[ defined as above for the couples of points {x(n, x 0 ), x(n, y 0 )}. For any A ∈ R k×k max , x, y ∈ R k max and λ ∈ R, we have A(λx ⊕ y) = λAx ⊕ Ay. As a consequence, the sequence ]λ(n), λ(n)[ is decreasing. Let λ and λ be the limits of λ(n) and λ(n). On the event A, we have λλ > c (see [START_REF] Brilman | Synchronisation by resources sharing : a performance analysis[END_REF]).

We define the sets

Let x, y, z ∈ PR k be three different points. It is immediate to prove that there exists a unique λ ∈ R such that z = λx ⊕ y. As a consequence, the sets Λ(n) are countable and Λ is countable. It implies that the set ]λ, λ[\Λ is non-empty on A. For all λ ∈]λ, λ[\Λ, we have, by definition of Λ, that S(λx 0 ⊕ y 0 ) ∩ S(x 0 ) = ∅. The conclusion follows.

⊓ ⊔ Remark The proof does not work when matrices A 1 , . . . , A p belong to R k×k max . In this case, it is possible to have S(x) = S(x). In the model detailed in Example 7.2, all the sets S(x) are dense in the interval [-1, 1] (as a classical consequence of the assumption η/η ′ ∈ Q). It implies that S(x) = [-1, 1], ∀x. The stationary distributions Q x are all defined on the same set, Ω × [-1, 1], which prevents the previous proof from working.

We want to prove the existence of a rank 1 pattern of {A(n)} (Def. 6.7). There exists a r.v. N such that A(N ) • • • A(0) ij > ε, ∀i, j (consequence of Equation ( 7)). It follows from the ergodic Lemma, that the set

is infinite, P -a.s. Let σ : N → N be the strictly increasing function such that I = {σ(0), σ(1), . . . }.

We define the subsequence {B(n

can be written under the form B(n

for some indices u, v belonging to the argmax in max ij B(n) ij . We also have ∀i, j, B(n

We consider the Euclidean space (PR k×k , |.| P ) where |.| P is the norm introduced in Definition 3.7. It follows from ( 16) and ( 17) that

We fix η > 0. Let C be such that ∀p > C, we have |B ∞ -B σ(p) | P η/2. For p > max(C, 2/η), we have B(B σ(p) , 1/p) ⊂ B(B ∞ , η). It implies

It means precisely that B ∞ is an asymptotic pattern of {A(n)}, see Definition 6.13.

Let us prove the sufficient part of the theorem. We assume that there exists a deterministic matrix à which is a rank 1 asymptotic pattern of {A(n)}. We want to prove the η-coupling convergence of π(x(n)) to a unique stationary regime. We fix η > 0. Let N η be such that

Using the ergodic Lemma, we have

Let u be the unique eigenvector of the rank 1 matrix à and B(u, η) the ball of center π(u) and radius η in PR k . We have that for all x 0 ∈ R k , {π(x(n, x 0 )) ∈ B(u, η)} ⊂ {πA(n -1)

In particular, it implies that ∀x 0 , y 0 ∈ R k and n large enough, {d(x(n, x 0 ), x(n, y 0 )) < η} ⊂ We deduce that P {d(x(n, x 0 ), x(n, y 0 )) < η} → 1. We conclude by using Lemma 9.1 (the existence of n such that P {D(A(n) • • • A(0) < +∞} > 0 comes from Equation ( 19)).

Proof of Theorem 7.5

We want to prove that the conditions given in Theorem 7.5 are necessary. We suppose that our model couples in finite time with a unique stationary regime, uniformly over initial conditions in R k . Let us prove a lemma first.

Lemma 9.5. If there is a unique stationary regime for π(x(n)), coupling in finite time uniformly over initial conditions in R k implies strong coupling in finite time uniformly over initial conditions in R k .


Fermi-Dirac Distribution

n k = 1 e k B T (k-µ) + 1 Fermi-Dirac distribution

Problems

Occupancy n k and photon m are independent

No fluctuation

Thermal light variance ⇒ No laser Mean values calculable Possible repartitions of electrons for a given energy throughout the paper and we are not interested in the pressure exerted by the electron gas on its boundaries. Thus subscripts introduced in thermodynamics to indicate that the volume is kept constant are omitted. Fig. 2. Ratio of the isolated electronic occupancy ͗N͘ from Eq. ͑6͒ and the Fermi-Dirac ͑FD͒ occupancy ͗N͘FD from Eq. ͑16͒ as a function of the FD occupancy, for the same ͑average͒ energy r taken as a parameter.

N-1,m+1 N+1,m-1 N-1,m N+1,m N,m-1 N,m+1 N,m A A J J E E C C B N=0 Π (N, m) n 2 k (N) = B N=0 Π (N, m+1) α + n 2 k (B -N) Detailed balance at m photons Π (m) = α + C 2 C 1 Π (m + 1)
Ratio between Microcanical and FD occupancies for the same energy r mκ(r ) = -i=1,2,...

(-1) i W [r -i(κ + i/2)]

Canonical ensemble

Hypothesis

We may use canonical occupancy at laser level.

n k (N, q) = Z (N -1, q)q k Z (N, q) [1 -n k (N -1, q)] Canonical Occupancy Z (N, q) = N i=1 q i-1 -q B 1 -q i Partition function n k (N + 1, q) = (1 -q N+1 )q k q N -q B [1 -n k (N, q)]
Recursive occupancy formula

Fig. 1 .

 1 Fig. 1. System containing seven electrons. The energy levels are labeled by , with ϭ0 at the Fermi level ͑defined one-half energy step above the zero-temperature highest energy͒. The column rϭ0 corresponds to zero temperature. The central part of the figure labeled rϭ6 exhibits the 11 ways ͑microstates͒ of incrementing the energy by Eϭ6⑀. The last column gives the number m of electrons occupying some energy level, read off the central part. The average number of electrons ͗N͘ is obtained by dividing m by W(6)ϭ11.

  Photons

  Microscopic: take thermalization into account. Satisfactory from physical point of view but was very slow.Canonical: shortcut in Markov chainSpeed up the program from 1h to 1s with the same physical accuracy.

		Canonical vs. Microscopic difference Dual-mode photon statistics Dual-mode photon statistics Dual-mode photon statistics Bimode photon Fano factor Bimode photon Fano factor Bimode photon Fano factor
				Microscopic Microscopic Microscopic Microscopic Microscopic Microscopic					Canonical Canonical Canonical Canonical Canonical Canonical
	Introduction Laser Introduction Laser Introduction Laser Introduction Laser Introduction Laser Introduction Laser Introduction Laser Introduction Introduction Laser Laser		500 500 500 var(m) <m> <m> <m> var(m) <m> var(m) <m>						500 500 500 var(m) <m> <m> <m> <m> var(m) <m> var(m) <m>
	overview overview overview overview overview overview overview overview overview
	Multimode Multimode Multimode Multimode Multimode Multimode Multimode Multimode Multimode
	Laser Laser models SC laser States Events Statistical specifications FD MC MS Laser Laser models SC laser States Events Statistical specifications Laser Laser models SC laser States Events Statistical specifications Laser Laser models SC laser States Events Statistical specifications Laser Laser models SC laser States Events Statistical specifications Laser Laser models SC laser States Events Statistical specifications Laser Laser models SC laser States Events Statistical specifications Laser Laser Laser models SC laser States Events Statistical specifications Laser models SC laser States Events Statistical specifications FD FD FD FD FD FD FD FD MC MC MC MC MC MC MC MC MS MS MS MS MS MS MS MS	E-E counts	300 400 300 400 300 400 40 60 40 60 40 60 200 200 200				Canonical Microscopic	300 400 300 400 300 400 40 60 40 60 40 60 200 200 200	Double emission events Canonical : Exponential distribution Microscopic : Non exponential distribution
	Canon. Conclusion Canon. Canon. Canon. Canon. Canon. Canon. Canon. Canon. Conclusion Conclusion Conclusion Conclusion Conclusion Conclusion Conclusion Conclusion		0 100 20 20 20 100 100	0.02	0.04	0.06	0.08	0.10	20 20 20 100 100 100	∃ Solution to restrain fast E-E within the Markov framework
					∆t a.u.
			0 0 0 0 0 0						0 0 0 0 0 0	GC mode equality GC mode equality GC mode equality GC mode equality GC mode equality GC mode equality
				Distributions of double emission events 0.61 0.62 0.63 α L 0.61 0.62 0.63 α L 0.61 0.62 0.63 α L 0.60 0.60 0.60 0.60 0.61 0.62 0.63 α L 0.60 0.61 0.62 0.63 α L 0.60 0.61 0.62 0.63 α L		0.60 0.60 0.60 0.60 0.60 0.60	0.61 0.61 0.61 0.61 0.61 0.61	0.62 0.62 0.62 0.62 0.62 0.62	α L α L α L α L α L α L	0.63 0.63 0.63 0.63 0.63 0.63	0.64 0.64 0.64 0.64 0.64 0.64	0.65 0.65 0.65 0.65 0.65 0.65
		Photon statistic for microscopic at pump rate 210 Photon statistic for microscopic at pump rate 210 Photon Fano factor for microscopic at pump rate 210 Photon Fano factor for microscopic at pump rate 210 Photon Fano factor for microscopic at pump rate 210		Photon statistic for canonical at pump rate 210 Photon statistic for canonical at pump rate 300 Photon statistic for canonical at pump rate 100
										16/20 17/20 17/20 17/20 18/20 18/20 18/20 19/20 20/20

<m> Photon statistic for microscopic at pump rate 210 Photon Fano factor for canonical at pump rate 210 Photon Fano factor for canonical at pump rate 300 Photon Fano factor for canonical at pump rate 100 Conclusion Monte-Carlo semiconductor laser simulation based on Markov chain requires electron presence at lasing levels → different statistical ensembles.

Grand-canonical: thermal light → nonsense for a laser.

Micro-canonical: accurately, but considers only closed systems.