Cheng-Hung Chang

Tyll Krüger

Roman Schubert

Serge Troubetzkoy

Quantisations of piecewise affine maps on the torus and their quantum limits

For general quantum systems the semiclassical behaviour of eigenfunctions in relation to the ergodic properties of the underlying classical system is quite difficult to understand. The Wignerfunctions of eigenstates converge weakly to invariant measures of the classical system, the so called quantum limits, and one would like to understand which invariant measures can occur that way, thereby classifying the semiclassical behaviour of eigenfunctions.

We introduce a class of maps on the torus for whose quantisations we can understand the set of quantum limits in great detail. In particular we can construct examples of ergodic maps which have singular ergodic measures as quantum limits, and examples of non-ergodic maps where arbitrary convex combinations of absolutely continuous ergodic measures can occur as quantum limits.

The maps we quantise are obtained by cutting and stacking.

Introduction

The correspondence principle in quantum mechanics states that in the semiclassical limit → 0 classical mechanics emerges and governs quantum mechanical quantities for small de Broglie wavelength. One manifestation of this principle is that the Wignerfunctions of eigenfunctions converge weakly to invariant probability measures on phase space, the so called quantum limits. It is one of the big open problems in the field to classify the set of quantum limits, and it is in general not known which invariant measures can occur as quantum limits. In particular the case that the classical system is ergodic has attracted a lot of attention. In this case the celebrated quantum ergodicity theorem, [START_REF] Šnirel ′ Man | Ergodic properties of eigenfunctions[END_REF][START_REF] Zelditch | Uniform distribution of eigenfunctions on compact hyperbolic surfaces[END_REF][START_REF] De Verdière | Ergodicité et fonctions propres du laplacien[END_REF], states that almost all eigenfunctions have the ergodic Liouville measure as quantum limit, and one would like to know if in fact all eigenfunctions converge to the Liouville measure, i.e., if quantum unique ergodicity holds or if there are exceptions. Possible candidates for exceptions would be quantum limits concentrated on periodic orbits, a phenomenon called strong scaring. Another very interesting case is when the classical system is of mixed type, i.e., the phase space has several invariant components of positive measure, or, if there exist several invariant measures which are continuous relative to Liouville measure. Here the question is to what extend the quantum mechanical system respects the splitting of the classical system into invariant components, i.e.., is a typical quantum limit ergodic, or can every convex combination of invariant measures appear as a quantum limit.

There has been recently considerable progress in some of these questions. For the cat map it was shown that quantum unique ergodicity does not hold, in [START_REF] Faure | Scarred eigenstates for quantum cat maps of minimal periods[END_REF] a sequence of eigenfunctions was constructed whose quantum limit is a convex combination of the Liouville measure and an atomic measure supported on a periodic orbit. It was furthermore shown that the orbit can carry at most 1/2 of the total mass of the measure [START_REF] Bonechi | Controlling strong scarring for quantized ergodic toral automorphisms[END_REF][START_REF] Faure | On the maximal scarring for quantum cat map eigenstates[END_REF]. The eigenvalues of the cat map have large multiplicities and this behaviour depends on the choice of the basis of eigenfunctions, in [START_REF] Kurlberg | Hecke theory and equidistribution for the quantization of linear maps of the torus[END_REF] it was shown that for a so called Hecke basis of eigenfunctions quantum unique ergodicity actually holds.

On compact Riemannian manifolds of negative curvature quantum unique ergodicity was conjectured in [START_REF] Rudnick | The behaviour of eigenstates of arithmetic hyperbolic manifolds[END_REF] and for arithmetic manifolds it was recently proved by Lindenstrauss for Hecke bases of eigenfunctions, [START_REF] Lindenstrauss | Invariant measures and arithmetic quantum unique ergodicity[END_REF]. The non-arithmetic case is still open, but in [START_REF] Anantharaman | Entropy and the localization of eigenfunctions[END_REF][START_REF] Anantharaman | Half-delocalization of eigenfunctions for the laplacian on an anosov manifold[END_REF] the authors succeeded in proving lower bounds for the entropy of quantum limits on manifolds of negative curvature.

In this paper we introduce a class of model systems for which the set of quantum limits can be determined very precisely. This work was motivated by a paper of Marklof and Rudnick were they gave an example of a quantum ergodic map which one can prove to be quantum uniquely ergodic [START_REF] Marklof | Quantum unique ergodicity for parabolic maps[END_REF]. They mention that there are no examples known where a quantum ergodic map is not quantum uniquely ergodic. The purpose of our work was to provide such examples, in fact examples which are quite close in nature to those considered by Marklof and Rudnick. The map they considered was a skew product map of the torus T 2 = R 2 /Z 2 of the form:

F : p q → p + 2q f (q) mod 1
where f (p) is an irrational rotation of the circle R/Z. In this article we consider skew products of the same form for other functions f (p). In particular since a circle rotation is an interval exchange transformation (IET for short) on two intervals, one of the examples we will consider for f (p) are interval exchanges on more intervals.

There are examples known of IETs which are not uniquely ergodic. A consequence of our main result is that if f (p) is an IET which is not uniquely ergodic then each of the invariant measures of F of the form Lebesgue measure cross an invariant measure (with the exception of finitely supported ones) of f (p) is a quantum limit. The plan of the paper is as follows. In Section 2 we give a quick review of quantisation of maps on the torus, and introduce the maps we study and their quantisation. In particular we prove Egorov's theorem for these maps. In Section 3 we turn our attention to quantum limits, we first give a general proof of quantum ergodicity for maps with singularities, and then show that for our particular class of maps the quantum limits can be understood purely in terms of the orbits of discretisations of the classical map. Then, in Section 4, we finally come to our main result. We first review the cutting and stacking construction to obtain maps and then show how it can be combined with discretisations to get a detailed understanding of quantum limits. In Theorem 4 we summarise our main findings. Finally in the last two sections we discuss two examples and give some conclusions.

Quantisation

We give a short summary of the quantisation of maps on the torus, for more details and background we refer to [START_REF]The mathematical aspects of quantum maps[END_REF][START_REF] Bièvre | Quantum chaos: a brief first visit[END_REF]. The Hilbert space: For (p, q) ∈ R 2 we introduce the phase space translation operator T(p, q) := e -i (q p-pq) , where pψ(x) := i ψ ′ (x) and qψ(x) := xψ(x) for ψ ∈ S(R), are the momentum and position operators, respectively. These operators are unitary on L 2 (R) and satisfy for (p, q), (p ′ , q ′ ) ∈ R 2

T(p + p ′ , q + q ′ ) = e -i 2 (qp ′ -pq ′ ) T(p, q)T(p ′ , q ′ ) (1)
and they provide therefore a unitary irreducible representation of the Heisenberg group on L 2 (R).

The state space of the classical map is obtained from R 2 by identifying integer translates which gives the two torus

T 2 = R 2 /Z 2 .
By mimicking this procedure the quantum mechanical state space is defined to be the space of distributions on R which satisfy

T(1, 0)ψ = ψ , T (0, 1)ψ = ψ .
One finds that these two conditions can only be fulfilled (for

ψ = const.) if Planck's constant meets the condition 1 2π = N (2)
where N is a positive integer. The allowed states then turn out to be distributions of the form

ψ(x) = 1 √ N Q∈Z Ψ(Q)δ x - Q N (3) with Ψ(Q) a complex number satisfying Ψ(Q + N ) = Ψ(Q) .
So the Ψ(Q) are functions on Z N = Z/N Z and the space of these functions will be denoted by H N , it is N -dimensional and through the coefficients

Ψ(Q), Q = 0, 1, • • • , N -1 it can be identified with C N . There is a map S N : S(R) → H N defined by S N ψ := n,m∈Z (-1) N nm T(n, m)ψ
which is onto. If we equip H N furthermore with the inner product

ψ, φ N := 1 N Q∈Z N Ψ * (Q)Φ(Q)
then H N is a Hilbert space and S N is an isometry.

Observables: In classical mechanics observables on the torus are given by functions on T 2 , these can be expanded into Fourier series

a = n∈Z 2 ân e(-ω(z, n)) .
where z = (q, p) ∈ T 2 , ω(z, n) = qn 2 -pn 1 and ân := T 2 a(z)e(ω(z, n)) dz denotes the n-th Fourier coefficient. We use here and in the following the notation e(x) = e 2πix and e N (x) := e 2πi N x . These observables can be quantised by replacing e(-ω(z, n)) by the the translation operator

T N (n) := T n 1 N , n 2 N
which acts on H N . This is called Weyl quantisation, to a classical observable a ∈ C ∞ (T 2 ) a corresponding quantum observable is defined by

Op N [a] := n∈Z 2 ân T N (n) , (4)
which is an operator on H N . For example, if a depends only on q then the corresponding operator is just multiplication with a,

Op N [a]ψ(q) = a(q)ψ(q) , (5)
and in terms of the coefficients Ψ(Q) the action of Op N [a] is given by Ψ

(Q) → a Q/N Ψ(Q).
The trace of a Weyl operator can be expressed in terms of the symbol, from (4) follows easily, see [START_REF]The mathematical aspects of quantum maps[END_REF][START_REF] Bièvre | Quantum chaos: a brief first visit[END_REF], that for a ∈ C ∞ (T 2 ) lim

N →∞ 1 N Tr Op N [a] = T 2 a dz . (6)
Quantisation of a map: Let

F : T 2 → T 2
be a volume preserving map. One calls a sequence of unitary operators U N : H N → H N , N ∈ N, a quantisation of the map F if the correspondence principle holds, i.e., if for sufficiently nice functions a one has

U * N Op N [a]U N ∼ Op N [a • F ] , (7)
for N → ∞. If this relation holds, it is often called Egorov's theorem and it means that in the semiclassical limit, i.e., for N → ∞, quantum evolution of observable approaches the classical time evolution.

Let us now turn to the specific class of maps we want to quantise. They are given by

F : p q → p + 2q f (q) mod 1 (8)
where f : [0, 1] → [0, 1] is a piecewise affine map given by a cutting and stacking construction which we will describe in detail in Section 4. For the construction of the quantisation we only need the property that the singularity set S ⊂ [0, 1] is nowhere dense. In order to quantise this map we proceed similar to the construction in [START_REF] Marklof | Quantum unique ergodicity for parabolic maps[END_REF], i.e., use a sequence of approximations to f . Consider the discretized interval

D N := {Q/N ; Q ∈ {0, 1, 2, • • • , N -1}} , (9)
i.e., the support of the Hilbert space elements (3). For each N ∈ N we will call a map f N : D N → D N an approximation of f if it is close to f in a certain sense which we will now explain. Since f is not assumed to be continuous we do not approximate it uniformly in the supremum norm. Let f (S) := {q : ∃q 0 ∈ S such that q = lim q ′ →q 0 f (q ′ )}. We measure the difference between f and an approximation f N only away from the set f (S). Let us call the relevant set

I ε := {q ∈ [0, 1]; dist(q, f (S)) ≥ ε}.
In the construction of f N in Section 4 we will choose a sequence ε N with lim N →∞ ε N = 0. For any fixed ε N the relevant measure for the quality of the approximation will be

δ N := δ N (ε N ) := sup Q/N ∈Iε N |f N (Q/N ) -f (Q/N )| . ( 10
)
Any approximation f N then defines via (8) an approximation F N of F . The quantisation of F is now defined to be the sequence of unitary operators

U N Ψ(Q) = e N -( f -1 N (Q)) 2 Ψ( f -1 N (Q)) , (11)
where fN (Q) := N f N (Q/N ) denotes the map induced by f N on Z N = Z/N Z. This is indeed a unitary operator on H N , with its adjoint given by

U * N Ψ(Q) = e N (Q 2 )Ψ( fN (Q)) .
That this sequence of operators U N is really a quantisation of the map F is the content of the Egorov theorem (7) which we will now prove. In our case we have to be careful at the singularities of the map. The singularities of f and F can be naturally identified, thus without confusion we can denote by S the set of singularities of F as well. By C ∞ S (T 2 ) we denote the space of functions in C ∞ (T 2 ) which vanish in a neighbourhood of F (S). We then find Theorem 1. For any a ∈ C ∞ (T 2 ) we have

U * N Op N [a]U N = Op N [a • F N ] ,
and for any a ∈ C ∞ S (T 2 ) there are constants C(a), ε 0 (a) > 0 such that for ε N < ε 0 (a)

||U * N Op N [a]U N -Op N [a • F ]|| ≤ C(a) δ N .
Proof. The map F and its quantisation U N can be decomposed into a product of two simpler maps and operators. Namely, with

F (1) (p, q) = (p + 2q, q) , and
F (2) (p, q) = (p, f (q))
we have

F = F (2) • F (1) .
These maps can be quantised separately as

U (1) N Ψ(Q) := e N (-Q 2 )Ψ(Q) , and
U
(2)

N Ψ(Q) := Ψ( f -1 N (Q)) ,
where f N denotes a discretisation of f on the Heisenberg lattice. We then have

U N = U (2) N U (1) N ,
and therefore it is sufficient to study the conjugation of an Weyl operator for the two operators separately. In the case of U

(1)

N it is well known that Egorov's theorem is exactly fulfilled

U (1) N * Op N [a]U (1)
N = Op N [a • F (1)
] , see [START_REF] Marklof | Quantum unique ergodicity for parabolic maps[END_REF]. For the study of the second operator we use that

Op N [a]U (2) N -U (2) N Op N [a] = 0
for any observable which is constant in q, therefore we can restrict ourselves in the following to the case that a is constant in p. But then Op N [a] is just multiplication with a, and we obtain

U (2) N * Op N [a]U (2) N ψ(q) = a(f N (q))ψ(q) .
For general observables we therefore obtain

U (2) N * Op N [a]U (2) N = Op N [a • F (2) ] ,
and this proves the first part of the theorem.

For the second part we have to estimate

|| Op N [a • F ] -Op N [a • F N ]||
and this can again be reduced to the case that a depends only on q, and then

Op N [a • F ] -Op N [a • F N ]
is the multiplication operator with a(f N (q)) -a(f (q)).

Since we have for a b ∈ C ∞ (T 2 ) which depends only on q by (5) that

|| Op N [b]|| = sup q∈D N |b(q)| , we obtain || Op N [a • f ] -Op N [a • f N ]|| = sup q∈D N |(a • f N -a • f )(q)| .
But for the right hand side we obtain by using (10) and that a ≡ 0 in a neighbourhood of f S

sup q∈D N |(a • f N -a • f )(q)| ≤ sup Q/N ∈Iε N |(a • f N -a • f )(Q/N )| + sup Q/N ∈I\Iε N |(a • f N -a • f )(Q/N )| ≤ C(a)δ N .
since the second term on the right hand side is 0 if ε N is small enough.

So if we can choose our approximations f N in a way that ε N → 0 and δ N → 0 for N → ∞, then the sequence of unitary operators U N reproduces the classical map F in the semiclassical limit N → ∞, and so the correspondence principle holds. Definition 1. A sequence of operators U N for which δ N (ε N ) and ε N tend to 0 for N → ∞ will be called a proper quantisation of F

The restriction on the support of the classical observables is necessary in order that a•F N and a•F are smooth for N large enough. For a general a the composition a • F is discontinuous which causes problems with the Weyl quantisation. Theorem 1 is not valid without the assumption on the singularities. This is shown by the following counter-example.

Proposition 1. Let s ∈ S, a(q) ∈ C ∞ (T 2 ) depend only on q with a(s) = 0 and let g s (q) := √ N e -(q-s) 2 /N be a Gaussian centred at s and ψ s := S N f s be its projection to H N . Then there exists a constant C such that

lim N →∞ || U * N Op N [a]U N -Op N [a • F ] ψ s || ≥ C|a(s)|||ψ s || Proof. Since a depends only on q the operator U * N Op N [a]U N -Op N [a • F ] is given by multiplication with a • f N -a • f , and we have || U * N Op N [a]U N -Op N [a • F ] ψ s || 2 = 1 N N n=1 (a • f N -a • f )(n/N ) 2 |Ψ s (n/N )| 2 , (12)
where Ψ s (q) = m∈Z g s (q -m). Since f is discontinuous at s there exists an ε > 0 and a C > 0 such that

(a • f N -a • f )(q) ≥ C|a(s)| , for q ∈ [s -ε, s + ε] . (13)
If we use now that g s is exponentially concentrated around q = s, which in particular implies Ψ s (q) = g s (q

) + O(e -c/N ) , for q ∈ [s -ε, s + ε] , (14)
we obtain

1 N N n=1 (a • f N -a • f )(n/N ) 2 |Ψ s (n/N )| 2 = 1 N |n/N -s|≤ε (a • f N -a • f )(n/N ) 2 |Ψ s (n/N )| 2 + 1 N |n/N -s|>ε (a • f N -a • f )(n/N ) 2 |Ψ s (n/N )| 2 ≥ C 2 |a(s)| 2 ||ψ s || 2 + O(e -c/N ) , (15)
where we have furthermore used

||ψ s || 2 = 1 N |n/N -s|≤ε |Ψ s (n/N )| 2 + O(e -c/N ).
We want to close this section with some comments about the underlying motivation for the specific quantisation assumptions on the neighbourhood of the singularities. Classically the singularities act like points with infinite local expansion rate respectively Lyapunov exponent. Therefore any perturbation in a small neighbourhood of the singularity set gives rise to an error which becomes unbounded if the perturbation approaches the singularity set. Since the quantised maps are a specific kind of perturbation it is natural to leave the allowed error big for points close to the singularity set.

Quantum limits and orbits

We will now discuss the implications of the Theorem 1 for the eigenfunctions of the quantised map. We will denote a orthonormal basis of eigenfunctions of

U N by ψ N k , k = 1, • • • , N , U N ψ N k = e N (θ N k )ψ N k ,
where we use the notation e N (x) = e 2πi N x and θ N k are the eigenphases. Each eigenfunction defines a linear map on the algebra of observables

Op N [a] → ψ N k , Op N [a]ψ N k
and the leading term for N → ∞ depends only on the principal symbol σ(a). The limit points of the sequence of all these maps defined by the eigenfunctions define measures on the set of classical observables and are called quantum limits (see, e.g., [START_REF] Marklof | Quantum unique ergodicity for parabolic maps[END_REF]). To put it more explicitly, a measure ν on T 2 is called a quantum limit of the system defined by the U N if there exist a sequence of eigenfunctions {ψ

N j k j } j∈N such that T 2 a(z) dν = lim j→∞ ψ N j k j , Op N [a]ψ N j k j .
One of the major goals in quantum chaos, and more generally in semiclassical analysis, is to determine all quantum limits that can occur and the relative density of the corresponding subsequences of eigenfunctions. We say that a subsequence of eigenfunctions {ψ

N j k j } j∈N has density α ∈ [0, 1] if lim N →∞ #{k j : N j = N } N = α , (16)
provided that the limit exists. Egorov's theorem usually implies that all quantum limits are invariant measures for the classical map. In our case the same is true, but we have to be careful at the singularities. If the set of singularities S is nowhere dense then the space C ∞ S (T 2 ) in Theorem 1 is large enough so that as an immediate consequence we have:

Corollary 1. Let us denote by M inv (F ) the convex set of F-invariant probability measures on T 2 , and by M qlim (U N ) the set of quantum limits µ of

U N with µ(S) = 0, then M qlim (U N ) ⊂ M inv (F )
So we only have to look at invariant measures as candidates for quantum limits. In the simplest case that there is only one invariant probability measure, i.e., that the system is uniquely ergodic, all eigenfunctions must converge to this measure, and we have the so called unique quantum ergodicity. This was the situation in the example of Marklof and Rudnick, [START_REF] Marklof | Quantum unique ergodicity for parabolic maps[END_REF].

We will study now the relationship between properties of quantum limits and the density of subsequences of eigenfunctions converging to them more closely. Our first result gives an upper bound on the density.

Theorem 2. Let U N be a proper quantisation of F and let µ be a quantum limit of U N with support Σ ⊂ T 2 . Then any sequence of eigenfunctions which converge to µ has at most density µ T 2 (Σ), where µ T 2 is the Lebesgue measure on T 2 . Proof. Let a ε ∈ C ∞ (T 2 ), ε ∈ (0, 1], be a sequence satisfying a ε | Σ = 1 and lim ε→0 a ε (z) = 0 for all z ∈ T 2 \Σ, i.e., a sequence approximating the characteristic function of Σ. If F = {ψ N j k j } j∈N is a sequence of eigenfunctions with µ as quantum limit, and

F N := {ψ N j k j ; N j = N }, then lim j→∞ ψ N j k j , Op N [a ε ]ψ N j k j = 1 (17)
and therefore

α(F) = lim N →∞ 1 N ψ∈F N ψ, Op N [a ε ]ψ ≤ lim N →∞ 1 N N k=1 ψ N k , Op N [a ε ]ψ N k = T 2 a ε dµ T 2 , (18)
where we have used (6). We now take the limit ε → 0 and the theorem follows.

Since µ T 2 (S) = 0 it follows in particular that a possible sequence of eigenfunctions converging to a quantum limit concentrated on S must have density 0. This result is as well interesting for non-ergodic maps, because it gives an upper bound on the number of eigenfunctions whose quantum limits are supported on an invariant subset Σ of T 2 by the volume of Σ.

In case that the system is ergodic, we can actually determine the quantum limit of most eigenfunctions.

Theorem 3. Let U N be a proper quantisation of F and assume that µ T 2 is ergodic. Then there exists a subsequence of eigenfunctions of density one which converges to µ T 2 . This is the usual quantum ergodicity result, but our proof differs from the standard one (see e.g. [START_REF]The mathematical aspects of quantum maps[END_REF]) in that we rely on the convexity definition of ergodicity, this is more convenient when dealing with maps with singularities as has been observed in [START_REF] Gérard | Ergodic properties of eigenfunctions for the Dirichlet problem[END_REF]. Recall that µ T 2 is ergodic if it is extremal in the convex set of invariant probability measures, i.e., if

µ T 2 = αµ 1 + (1 -α)µ 2 with µ 2 = µ T 2 then α = 1 and µ 1 = µ T 2 .
Proof. The existence of a subsequence F = {ψ

N j k j } j∈N of density one of eigenfunctions with quantum limit µ T 2 is equivalent to lim N →∞ 1 N N k=1 | ψ N k , Op N [a]ψ N k -a| = 0 , (19)
see [START_REF]The mathematical aspects of quantum maps[END_REF]. We first observe that by ergodicity every subsequence F = {ψ

N j k j } j∈N of positive density satisfies lim N →∞ 1 |F N | ψ∈F N ψ, Op N [a]ψ = a . (20)
To see this we consider the sequence

a N := 1 |F N | ψ∈F N ψ, Op N [a]ψ , (21)
this is a bounded sequence since Op N [a] is bounded, and therefore there exists a convergent subsequence {a N j } j∈N . Now using (6) we have with a convergent subse-quence

a = lim j→∞ 1 N j N j k=1 ψ N j k , Op N [a]ψ N j k = lim j→∞ |F N j | N j 1 |F N j | ψ∈F N j ψ, Op N [a]ψ + lim j→∞ N j -|F N j | N j 1 N j -|F N j | ψ N j k ∈H N j \F N j ψ N j k , Op N [a]ψ N j k = α(F) a dµ 1 + (1 -α(F)) a dµ 2 (22)
where µ 1 and µ 2 are invariant measures defined by

lim j→∞ 1 |F N j | ψ∈F N j ψ, Op N [a]ψ = a dµ 1 (23) lim N j →∞ 1 N j -|F N j | ψ N j k ∈H N j \F N j ψ N j k , Op N [a]ψ N j k = a dµ 2 (24)
These two measures exist by the assumption that the subsequence {a N j } j∈N was convergent, and they are invariant by Theorem 1. But equation ( 22) can be rewritten as

µ = αµ 1 + (1 -α)µ 2 (25)
and if µ is ergodic and α = 0 this is only possible if µ 1 = µ, and this proves that

lim j→∞ a N j = a . (26)
Since this holds for every convergent subsequence of {a N } N ∈N a is the only limit point and (20) follows. Now assume that lim

N →∞ 1 N N k=1 | ψ N k , Op N [a]ψ N k -a| = C > 0 , (27)
then there must either exists a subsequence {k j } j∈N of positive density with

ψ N j k j , Op N [a]ψ N j k j -a ≥ C/2 (28)
or one with ψ

N j k j , Op N [a]ψ N j k j -a ≤ -C/2 . ( 29
)
But the mean value of the sequence ψ

N j k j , Op N [a]ψ N j
k j -a must tend to 0 by (20) and so we have a contradiction if C = 0.

The previous results, Corollary 1, Theorem 2 and Theorem 3, are quite general, they are valid for all quantised maps which satisfy Egorov's theorem. We now turn to a more concrete study of the eigenfunctions for the specific quantum maps (11). Our aim is to show that the quantum limits are determined by the spatial distribution of the periodic orbits of the discretisation of the classical map. The eigenvalue equation U N ψ = e N (θ)ψ leads to the following recursion equation for ψ

Ψ( fN (Q)) = e N (θ -Q 2 )Ψ(Q) . ( 30
)
From this recursion relation we obtain

|Ψ( fN (Q))| 2 = |Ψ(Q)| 2 ,
and this implies that the probability densities in position space defined by the eigenfunctions are invariant under the map f N . In order to determine these densities it is therefore sufficient to determine the spatial distribution of the orbits of f N .

For the further investigation we note that each periodic orbit of f N carries at least one eigenfunction. And we can determine the eigenfunctions and eigenvalues more explicitly, let O be an periodic orbit of period |O| = K, then the recursion relation (30) gives

Ψ(Q) = e N Kθ - K-1 k=0 f k N (Q)] 2 Ψ(Q) .
So if ψ should be an eigenfunction with eigenvalue e N (θ) we get the condition

Kθ - K-1 k=0 f k N (Q)] 2 = N m with m ∈ {0, 1, • • • , K -1}.
This determines the eigenvalues, and then the corresponding eigenfunctions follow from the recursion relation and the normalisation condition. Summarising we get: Proposition 2. Let O be an orbit of period K = |O| of f N , then there exists K orthogonal eigenfunctions of U N with support O. The eigenphases are given by

θ k = S O + N K k with k ∈ {0, 1, • • • , K -1} and
S O := 1 K K-1 k ′ =0 f k ′ N (Q)] 2 ,
and a normalised eigenfunction corresponding to θ k is given by

Ψ k ( f k ′ N (Q 0 )) = e N k ′ θ k - k ′ m=0 [ f m N (Q 0 )] 2 N K 1/2
where Q 0 ∈ O is an arbitrary point on the orbit and

k ′ ∈ {0, 1, • • • , K -1}.
The quantum lattice D N of N points is a disjoint union of all periodic orbits of f N , and on each of these orbits are as many eigenfunctions concentrated as the orbit is long. But that means that the orbits determine the quantum limits and the relative density of the corresponding sequence of eigenfunctions.

To each periodic orbit O we can associate a probability measure on [0, 1]

δ O (q) := 1 |O| Q∈O δ q - Q N ( 31
)
which is invariant under f N .

Corollary 2. Let O (N ) j , j = 1, . . . , J N be the periodic orbits of f N and let δ (N ) j be the corresponding probability measures (31). Assume that there is an invariant measure ν of f and a sequence of periodic orbits {O

(N k ) j k } k∈N such that lim k→∞ δ (N k ) j k = ν , then ν is a quantum limit of U N .
In all our examples the sequence (N k ) contains all natural numbers. Thus there are two possible definitions of the density of a sequence of periodic orbits, G = {O

(N k ) j k } k∈N , α(G) = lim N →∞ |O (N ) j N | N or β(O (N k ) j k ) = lim N →∞ 1 N k:N k =N |O (N k ) j k | ,
whenever the limit exists, which we will call the α-density or β-density of G, respectively. This corollary suggests that the set of quantum limits coincides with the set of limits points of the sequence of orbit measures δ (N ) j and that moreover the relative densities of the convergent subsequences coincide too. This is true if there are no multiplicities in the eigenvalues. If there are eigenvalues of multiplicity larger than one, then the eigenspace can mix the contribution of the different orbits. But even in this case there always exists a choice of a basis of eigenfunctions corresponding to the orbit measures δ (N ) j . Notice that in this case the β-density of the sequence of orbits coincides with the density of the corresponding sequence of eigenfunctions defined in (16).

Cutting and Stacking constructions

Cutting and stacking is a popular method in ergodic theory to construct maps on the interval which are isomorphic models of arbitrary measure preserving dynamical systems. The construction gives a piecewise isometric mapping on the interval with Lebesgue measure as an invariant measure. One can also think of this transformation as a countable interval exchange transformation.

The final mapping will be defined only Lebesgue almost everywhere. None the less we can use this model to study certain other invariant measures which are well behaved with respect to the cutting and stacking construction.

For a very readable introduction into cutting and stacking see the recent book by Shields [START_REF] Shields | The ergodic theory of discrete sample paths[END_REF] or the old book by Friedman [START_REF] Friedman | Introduction to ergodic theory[END_REF]. We will now give a short description of the basic construction scheme and the relevant definitions.

A stack (or column) S is a finite family of enumerated disjoint intervals {I j }

h(S) j=1

where h(S) is called the height of S. The I j are subintervals of [0, 1] of equal length which is called the width of S. There are two possible conventions: either one can take all the intervals to be open, or all the intervals to be closed on the left and open on the right. Which convention we choose is not important for this paper. The intervals I 1 and I h(s) are called the bottom and top of S respectively. We define a transformation f S as follows, if the point x ∈ I j is not on the top of S and not on the boundary of I j then it gets mapped to the point directly above it (see Figure 1(a)). Since I j+1 and I j have equal width, f S is simply the canonical identification map between I j and I j+1 . Interpreted in [0, 1] this means that f S : I j → I j+1 such that x → x + ∂ -I j+1 -∂ -I j where ∂ -denotes the left boundary point of an interval. The construction clearly defines f -1 on all I j except at the bottom.

A stack family S is a finite or countable set of stacks

{S i } = {{I i j } h(S i )
j=1 } such that all I i j are disjoint and ∪I i j = [0, 1]. In this paper we will work only with finite stack families. On S one defines a transformation f S by f S | S i = f S i except on the collection of top intervals.

Given two stacks S i and S j , i = j of equal width one can define a new stack S ′ by stacking S j on S i that is

S ′ = {I ′ k } h(S i )+h(S j ) k=1 I ′ k = I i k for k ≤ h(S i ) and I ′ k = I j k-h(S i ) for k > h(S i ).
Correspondingly one gets a new transformation f S ′ which agrees with f S i and f S j except on I i h(S i ) , where f S i was not defined before. It remains to define the cutting of stacks. A cutting of a stack S = {I k } is a splitting of S into two (or more) disjoint stacks S 1 and S 2 with intervals {I 1 k } and 1(b)). The definition of f {S 1 ,S 2 } is as above. Multiple cutting of S is defined analogously.

{I 2 k } such that I 1 k ∪ I 2 k = I k and ∂ -I 1 k < ∂ + I 1 k = ∂ -I 2 k < ∂ + I 2 k ∀k that is I 1 k is always the left component of the partition of I k into I 1 k and I 2 k (Fig- ure
A stack family S(n) is obtained from a stack family S(n -1) by cutting and stacking, if each S i (n) from S(n) can be obtained by successive cuttings and stackings of stacks from S(n -1). By construction f S(n) is an extension of f S(n-1) . If one has a sequence {S(n)} n≥1 of stack families such that each S(n) is obtained from S(n -1) by cutting and stacking and furthermore This gives an extension f S ′ of the map f S which was not defined on the top of substack two before.

lim n→∞ S k (n)∈S(n) width(S k (n)) = 0, (32)
then lim f S(n) = f is an invertible transformation on [0, 1] defined everywhere except at a set of zero Lebesgue measure. Note that f is always aperiodic.

The "partition" of [0, 1] into the intervals of S(1),1 the starting object of the construction, gives a natural symbolic dynamics for f . The coding is unique for all points whose infinite orbit is defined.

For convenience, we denote the intervals of S(n) by I i j (n) where i indexes the stack and j the interval in the stack. We consider the set

S = ∪ i,n ∂I i h(S i ) (n) ∪ ∪ i,n ∂I i 1 (n) .
The set S = S (f ) is called the singularity set of the map f . It consists of all the points of discontinuity of f and all the points where the map f or f -1 is not defined. The boundary points of the intervals I i j which are not top or bottom intervals are not included in this set, the map is defined and continuous on such points! Furthermore

∪ i,j,n ∂I i j (n) is included in the set ∪ k∈Z f k S, where f k S := {x : f -k x ∈ S}.
Historically the cutting and stacking construction was invented to represent the dynamics with respect to a single invariant measure as a countable interval exchange transformation with the canonical invariant Lebesgue measure. The construction is universal in the sense that for every measurable dynamical system (M, g, µ) one can explicitly give a cutting and stacking representation ([0, 1] , f (g) , µ L ) [START_REF] Shields | The ergodic theory of discrete sample paths[END_REF]. The following proposition which will be needed for the application of Theorem 3 seemed to be unknown. Proposition 3. Let (Σ, σ) be a symbolic dynamical system over a finite or countable alphabet and let µ be a shift-invariant measure such that µ ([w]) > 0 for all cylindersets corresponding to finite words w appearing in Σ. Then the associated cutting and stacking construction yields a representation ([0, 1] , f, µ L ) of (Σ, σ, µ) such that all nonatomic invariant measures ν * on (Σ, σ) have a corresponding isomorphic invariant measure ν on ([0, 1] , f ) with the property ν (S (f )) = 0.

Proof. By definition all finite symbolic words have a representation as orbit segments of the cutting and stacking construction. Furthermore to each nonperiodic symbol sequence corresponds a unique point in [0, 1] whose orbit under f is well defined and commutes with the shift. By the Poincaré recurrence theorem every f -invariant measure is supported on recurrent points. Since the singular orbits, respectively symbolic sequences, are the ones which eventually or asymptotically fall onto the singularity set they do not intersect with the recurrent nonperiodic symbol sequences. Therefore all invariant measures-except the finitely supported ones-give zero measure to the singularity set.

The approximating family

Each f defined by cutting and stacking provides us with a natural approximation family {f S(n) } which we will use now to define the approximation mappings on the rational points D N = { Q N : Q ∈ {0, ....., N -1}} for the quantisation. Let the points in G i,j (n, N ) := D N ∩ I i j (n) be enumerated from left to right and let

K(N, I i j (n)) := ♯G i,j (n, N ) and K(N, S i (n)) := min j {K(N, I i j (n))}. K(N, S i (n)
) is just the smallest number of points from the discretisation in an interval in the stack S i (n). Let Ĝi,j (n, N ) be the set of the first

K(N, S i (n)) points from G i,j (n, N ) denoted by {x j,i e } K(N,S i (n)) e=1 . We define f N,n first on DN,n = ∪ i,j Ĝi,j (n, N ) by setting f N,n x j,i e = x j+1,i e for j < h i (n) := h (S i (n)) . (33)
We call these the internal orbit segments. Clearly f N,n x j,i e -f x j,i e = O 1 N . We call each approximation mapping f N,n on D N whose restriction to DN,n is given by the above construction an ergodic approximation.

Let ĎN,n be the set of points not in DN,n and not in in any of the top intervals The approximate mapping f N on the discrete set D N (denoted by the full dots), and the way the orbit segments from different stacks are concatenated in order to produce a quantum limit of the form α 1 µ 1 + α 2 µ 2 . A singular limit µ * is obtained by a stack of small width and where the orbits are concatenated from top to the bottom. necessarily an invertible map, thus the construction implies that max

I i h(S i ) (n). For x ∈ ĎN,n let f N,n x be the closest point to f x. Note that f N,n is not (x) f(x)
x∈ ĎN,n |f N,n x-f x| = O 1 N . It is clear that for fixed n we have lim N →∞ ♯ DN,n ♯D N = 1. (34)
To complete the definition of f N,n , it remains to define the mapping of the points Ĝi,h i (n) (n, N ) on the tops of the stacks to points Ĝi,1 (n, N ) on the bottoms of the stacks. This will be done in a way to produce periodic orbits which approximately mimic a given invariant measure. Furthermore it remains to link N to a given n to get a good approximation. In essence we have to require that for each fixed stage n construction we have enough discretisation points in each stack. That means that with the increase of N we pass from n to n+1 only as we pass a critical threshold value N n . This can be done already without a precise description of the gluing between top and bottom of the stacks. We need that the approximation family f N is good enough to apply Theorem 1 for a sequence ε N → 0 such that δ N (ε N ) → 0. Note that Theorem 1 does not impose any requirements on the rate of convergence. The basic idea is to keep N large enough compared with n such that all intervals I i j (n) of S(n) contain sufficiently many points from D N . Let b n = min i width(S i (n)). Choose a function n(N ) going to infinity such that lim

N →∞ min i K(N, S i (n(N ))) = ∞ (35) which is equivalent to require 1 bn = o(N ). Furthermore let ε N = max i width(S i (n)).
Equation (32) implies that ε N tends to 0 since n(N ) tends to infinity. With this choice of ε N the points where

f N,n(N ) is not yet defined do not contribute to δ N (ε N ), hence one obtains δ N (ε N ) = O( 1 N
). An approximation family f N := f N,n(N ) for any function n(N ) satisfying the above requirements is called proper.

We say that a measure µ appears as a quantum limit if one can find a proper approximating family f N and associated quantisation U N such that µ is a quantum limit of U N . The notion of quantum limit as well as the notion of density where introduced in Corollary 2.

Theorem 4. a) If µ is an absolute continuous ergodic measure for f , then µ appears as a positive α and β-density quantum limit. Furthermore the quantum limit has full α and β-density if µ is the Lebesgue measure.

b) If µ is a nonatomic, singular ergodic measure for f then µ appears as a quantum limit. Furthermore the quantum limit must have zero α and β-density. c) If µ 1 and µ 2 are two absolute continuous ergodic invariant measures, then α 1 µ 1 + α 2 µ 2 appears as a positive α and β-density quantum limit for any α 1 , α 2 ∈ (0, 1) with α 1 + α 2 = 1.

d) If µ 1 and µ 2 are two ergodic measures at least one of which is singular, then α 1 µ 1 + α 2 µ 2 appears as a quantum limit for any α 1 , α 2 ∈ (0, 1) with α 1 + α 2 = 1. Furthermore the quantum limit must have zero α and β-density.

A transformation f : [0, 1] → [0, 1] is called a finite rank transformation if one can construct it via cutting and stacking such that the number of stacks #{S i (n)} in the n-th stack family S(n) is bounded (independent of n).

Corollary 3. If f is a finite rank transformation, then every ergodic non atomic invariant measure appears as the quantum limit of any proper ergodic approximation family.

Proof. For a finite rank transformation, the singularity sent S is a countable set which has only a finite number of points of density. Thus any non atomic invariant measure can not be support on S.

A point x is called µ-typical if lim m,u→∞ 1 m+u+1 m i=-u δ f i x = µ.
For the proof of the theorem we need the following simple fact whose proof is omitted since it is immediate from the definition of weak convergence of measures. Proposition 4. Fix x be µ-typical. Let j(m, u) and ε(j) be functions such that j(m, u) → ∞ for m, u → ∞ and ε(j) → 0 as j → ∞. Let {y (j) k } k∈Z be a family of sequences with the property that y

(j) k -f k x ≤ ε(j) for -u ≤ k ≤ m. Then lim m,u→∞ 1 m+u+1 m k=-u δ y (j(m,u)) k = µ.
Proof of Theorem 3. For the proof of part b) and part a) for α-density we complete the definition of f N as follows. Let us complete each internal orbit segment into a periodic orbit by setting f N x h i (n),i e = x 1,i e (compare (33)) . For points x ∈ D N \ S where the map is not yet defined we have the freedom to map x anywhere, for preciseness define f N x to be the closest point to f x. By Corollary 2 it is enough to show that there is a sequence of periodic orbits on D N whose point-mass average converge for N → ∞ to the considered measure µ.

Fix an arbitrary enumeration O be the set of periodic orbits in ∪j Ĝi,j (n(N ), N ). They stay width(S i (n(N )))-close to the orbit segment {f -m 0 x, ....., f n 0 x} where m 0 and n 0 are the smallest and largest iterates such that f k S i (n(N )) x is still defined. Note that

m 0 + n 0 + 1 = h (S i (n(N ))). Since width(S i (n(N ))) ≤ ε N := max i width(S i (n(N ))) → 0 for N → ∞
we can apply the above proposition to the family O (N ) j l (with fixed l) to conclude that µ is a quantum limit, however this construction has not yet proved the positive density.

To prove the positive density we need a quantified version of the above. Let J q be the set of subintervals of [0, 1] with boundary points of the form p q . For the convergence of a sequence to a measure it is clearly enough to check the characteristic function averages with respect to the elements of ∪ q J q . A stack S i (n(N )) is called ε -q-good with respect to µ if for all x ∈ I i 1 (n(N )) (i.e. x in the base of the stack S i (n(N )))

µ(J) -ε ≤ 1 h i (n(N )) h i (n) i=0 1 J f | i S i (n(N )) (x) ≤ µ(J) + ε f or ∀J ∈ J q
where h i (n(N )) denotes the height of the stack S i (n(N )). Denote the family of such stacks by G(n, q, ε, µ) and by G(n, q, ε, µ) the set of points contained in G(n, q, ε, µ).

Clearly one has ∀q, ε > 0

lim n→∞ µ L (G(n, q, ε, µ)) = µ L (x : x is µ -typical) and lim n→∞ µ(G(n, q, ε, µ)) = 1.
Let ε(n) be a sufficiently slowly decreasing function and q(ε(n)) be a sufficiently slowly increasing function that

lim n→∞ µ L (G(n, q(n), ε(n), µ)) = µ L (x : x is µ -typical) . (36)
With this new notation we are ready to prove that µ is a quantum limit of positive density. We define the sequence O (N ) j l of periodic orbits which give rise to the desired quantum limit as follows. For fixed N the set G(n(N ), q(n(N )), ε(n(N )), µ)∩ DN,n(N) consists of a collection of points of periodic orbits. The sequence O (N ) j l consists of the set of these orbits. The positive β-density and full density in the case of Lebesgue measure then follows from Equation (36) and 34. To prove part a) for the α-density we only need to modify the map f N on top of the stacks. This will be done such that the collection of periodic orbits O (N ) j l becomes just one periodic orbit for each N. This completes the proof of part a).

To prove part b) observe that due to Proposition 3 we have µ (S) = 0 and hence can apply Theorem 1 to get an invariant measure out of the quantum-limit. One considers the set of µ-typical points. From the proof of Proposition 3 follows that the orbit of every µ-typical point does not intersect nor converge to the singularity set S. Since the set of µ-typical points has zero Lebesgue measure we can obtain only a zero density quantum limit just as in the proof of part a).

To prove part c) one has to modify the construction of the approximating mapping f N in the following way. Instead of making f N periodic within each stack S i (n(N )) we want to connect two stacks say S i (n(N )) and S j (n(N )) where the orbit segments in the i-th stack respectively j-th stack are approximately typical for µ 1 respectively µ 2 to get an average of µ 1 and µ 2 .

For

l = 1, 2 let A l (N ) := {i : S i (n(N )) ∈ G(n (N ) , q(n), ε(n), µ l ). On DN,n(N) de- fine f N as before by f N x j,i e = x j+1,i e for j < h i (n). For x j,i e ∈ G(n (N ) , q(n), ε(n), µ l )∩ DN,n(N) one has for ∀J ∈ J q(n) µ l (J) -ε(n) + O( 1 N ) ≤ 1 h i (n) 0≤k≤h i (n)-1 1 J f k N x 1,i e ≤ µ l (J) + ε(n) + O( 1 N ).
Let θ l := µ L (x : x is µ l -typical) and note that

#A 1 (N ) ∩ DN,n(N) #A 2 (N ) ∩ DN,n(N) → θ 1 θ 2 for N → ∞.
Thus by gluing all the orbit segments of f N in the sets A 1 (N )∩ DN,n(N) and A 2 (N )∩ DN,n(N) in such a way that they form one periodic orbit we obtain a family of periodic orbits with quantum limit α 1 µ 1 + α 2 µ 2 where α l = θ l θ 1 +θ 2 . The α and β-densities are just θ 1 + θ 2 .

It is easy to construct in the same spirit approximation families f N for any values α 1 and α 2 = 1-α 1 . Suppose first that α 1 < θ 1 θ 1 +θ 2 and hence α 2 > θ 2 θ 1 +θ 2 . Take in each stack S i (n(N )) with i ∈ A 1 (N ) approximately α 1 α 2 of the internal orbit segments. The function n (N ) is sufficiently slowly growing (35) to ensure that there are enough points in the discretisation set DN,n(N) we can guarantee the convergence to α 1 α 2 . Gluing these segments together with all the internal orbit segments of A 2 (N ) yields a single periodic orbit O (N ) . The family of periodic orbits O (N )

N defines a quantum limit for the measure α 1 µ 1 + α 2 µ 2 with α and β-density α 1 α 2 θ 1 + θ 2 > 0. The case α 1 > θ 1 θ 1 +θ 2 is analogous.

The proof of d) follows immediately by combining the arguments from parts b) and c).

Examples

Interval Exchange Maps

Consider a permutation π of {1, 2, . . . , n} and a vector

v = (v 1 , . . . , v n ) such that v i > 0 for all i and n i=1 v i = 1. Let u 0 = 0, u i = v 1 + • • • + v i and ∆ i = (u i-1 , u i ). The interval exchange transformation T = T π, v , T : [0, 1] → [0, 1]
is the map that is an isometry of each interval ∆ i which rearranges these intervals according to the permutation π.

The Lebesgue measure is always an invariant measure for an IET. A typical IET is uniquely ergodic, however there exist minimal, non uniquely ergodic IETs. The first example of a minimal, non uniquely ergodic IET was given by Keynes and Newton [START_REF] Keynes | A "minimal", non-uniquely ergodic interval exchange transformation[END_REF] and Keane [START_REF] Keane | Non-ergodic interval exchange transformations[END_REF]. The number of ergodic invariant measures for a minimal IET on m intervals is bounded by the ⌊m/2⌋. [START_REF] Katok | Invariant measures of flows on orientable surfaces[END_REF][START_REF] Veech | Interval exchange transformations[END_REF]. The set of invariant measures always includes absolutely continuous measures but can also include singular measures. It is known that an interval exchange transformation on m intervals is at most of rank m, in particular it is a finite rank transformation (see for example [START_REF] Ferenczi | Systems of finite rank[END_REF]). In fact the typical IET is of rank 1 [Vee84], although we will not use this fact. Thus we can apply Corollary 3 to conclude:

1. any uniquely ergodic IET is quantum uniquely ergodic, 2. any minimal, non uniquely ergodic IET is not quantum uniquely ergodic, 3. any absolutely continuous invariant measures appear as a positive density quantum limit, 4. any singular ergodic invariant measure appears as a zero density quantum limit.

The full shift

Another example of a cutting and stacking transformation f B that has µ L as an ergodic invariant measure and admits further singular measures µ such that µ(S) = 0 is given by the full shift. Take any cutting and stacking model of the full two-sided shift on two symbols with Bernoulli-measure p 0 = p 1 = 1 2 (for details of such models we refer to the book [START_REF] Shields | The ergodic theory of discrete sample paths[END_REF]). Note that although the full shift has many periodic orbits the cutting and stacking model has none. We remark that one could introduce some periodic orbits at the boundaries of the subintervals but they would all sit or fall at singularity points and hence do not appear as quantum limits, in other words there are no scars in quantised cutting and stacking skew product mappings. By Proposition 3 all other invariant measures of the full shift have no support on the singularity set. Hence we can apply Theorem 3. It is interesting to note, that the fractal-dimensions (box or Hausdorff dimension) of the singularity set are rather large and that the upper and lower dimensions do not coincide. A straightforward counting argument shows for instance that the upper and lower box dimensions are in the open interval 1 2 , 1 .

Comments and Conclusions

We have shown in this paper that for a rather general class of dynamical systems on the torus the variety of different invariant measures can be recovered as quantum limits of the corresponding proper families of quantised maps. The quantisation scheme here used is based on the one introduced by [START_REF] Marklof | Quantum unique ergodicity for parabolic maps[END_REF]. For a discussion of alternative quantisation procedures and a critical comparison we refer to the recent work [START_REF] Zelditch | Quantum maps and automorphisms, The breadth of symplectic and Poisson geometry[END_REF].

One of the main features in our systems is the presence of singularities. In the quantisation procedure this provides enough freedom to obtain eigenfunctions reflecting the typical orbit structure with respect to any non atomic ergodic measure. It is an interesting question whether our results are still valid in case the classical dynamical system has no singularities. We conjecture that similar statements can be obtained. For this it seems natural to replace the top-bottom gluing scheme in the interval exchange approximating family by cutting and "crossover-concatenation" of touching period orbits.

Concerning the quantisation of flows one might hope that a good understanding of the associated quantised Poincare maps can guide one to a deeper understanding of concrete features of eigenfunctions and spectrum. An natural class of examples to study this questions are polygonal billiards. In the case of rational polygons the associated Poincare maps for the directional flow are interval exchange transformations which can be quantised similar to the quantisation used in this paper. It would be interesting to compare the results obtained that way with the semiclassical properties of the direct flow quantisation via the billiard Hamiltonian.

Figure 1 :

 1 Figure1: The cutting and stacking construction. In a given stack (a) the mapping f S is defined, except at the top interval. In (b) the stack is cut into three substacks, and in (c) the third substack is stacked onto the second one. This gives an extension f S ′ of the map f S which was not defined on the top of substack two before.

  Figure 2: The approximate mapping f N on the discrete set D N (denoted

  orbits on D N . For µ absolute continuous the set of points x ∈ [0, 1] with lim m,u→∞ 1 m+u+1 m i=-u δ f i x = µ has positive Lebesgue measure and in the case µ is the Lebesgue measure it has full measure. Let x be µ-typical and consider for each N the stack S i (n(N )) in which x is placed, where n(N ) is a function satisfying the requirements of Equation (35). Let O (N ) j l

Here we can ignore the boundary points of the interval since f is not defined on them.


Vincent Miele
email: vincent.miele@univ-lyon1.fr

Catherine Matias

Revealing the hidden structure of dynamic ecological networks

Keywords: Dynamic networks, Network clustering, Stochastic block model, Animal contact network, Trophic network

1. Recent technological advances and long-term data studies provide interaction data that can be modelled through dynamic networks, i.e a sequence of different snapshots of an evolving ecological network. Most often time is the parameter along which these networks evolve but any other one-dimensional gradient (temperature, altitude, depth, humidity, . . . ) could be considered.

2. Here we propose a statistical tool to analyse the underlying structure of these networks and follow its evolution dynamics (either in time or any other one-dimensional factor). It consists in extracting the main features of these networks and summarise them into a high-level view.

3. We analyse a dynamic animal contact network and a seasonal food web and in both cases we show that our approach allows for the identification of a backbone organisation as well as interesting temporal variations at the individual level.

4. Our method, implemented into the R package dynsbm, can handle the largest ecological datasets and is a versatile and promising tool for ecologists that study dynamic interactions.

Introduction

Networks are widely used in ecology as they provide a powerful tool for modelling the complex interplay between ecological entities. Depending on the context, those entities can be different species or different individuals while their interplay may be as diverse as trophic, competitive, cooperative relations or even contacts measured through physical proximity. Studying these networks can help answering important ecological questions about for e.g. the structure of these interactions and their robustness to external factors.

As Newman and Leicht pointed out some 10 years ago "much of the current research on networks [...] aimed at answering the question: how can we tell what a network looks like, when we can't actually look at it? ". One first answer has been to develop and use descriptive statistics and network measures such as connectance or centrality (see Rayfield et al., 2011, for a comprehensive list). This approach considers any ecological network as a whole, assuming that the network is homogeneous. A next proposal has been to go beyond descriptive statistics and consider network clustering, i.e. grouping entities according to their common properties. This technique allows answering fundamental questions about any underlying network structure: is there any peculiar non-random mixing of entities that would be a sign for a structural organisation [START_REF] Kéfi | How structured is the entangled bank? The surprisingly simple organization of multiplex ecological networks leads to increased persistence and resilience[END_REF]? is there, for instance, compartmentalisation [START_REF] Montoya | Functional group diversity increases with modularity in complex food webs[END_REF], hierarchical organisation [START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF] or nestedness [START_REF] Bascompte | The nested assembly of plant-animal mutualistic networks[END_REF]?

Nowadays, recent technological advances (sensors, GPS technology, . . . ) and long-term data studies have given rise to an avalanche of temporal data that need to be appropriately modelled. Data acquired over time can be aggregated within relevant time intervals (days, seasons, years, . . . ) and consequently produce snapshots of a same ecological network at different time steps. With these new data, one can potentially address new ecological questions which might not be tackled through the analysis of the static network where data is aggregated over the full recording time. In the same way, snapshots of an ecological network along any one-dimensional factor (such as temperature, altitude, depth, humidity, . . . ) may help analyse the evolution of the network structure along this gradient [START_REF] Stegen | Evolving ecological networks and the emergence of biodiversity patterns across temperature gradients[END_REF]. However addressing those new questions requires the development of new methodological tools. Up to now, very few proposals have been made to handle what we call here "dynamic networks", namely any sequence of snapshots of a same ecological network along a one-dimensional parameter (that we most often call time).

The two fundamental questions we will focus on here are the following: are there any relevant statistical patterns in the dynamic network? If so, how does this structure vary with time (or along the sequence)? In this article, we answer these two key points and argue that this is a first stone for further understanding and predicting processes on dynamic ecological networks such as event spreading (infection or extinction, for instance).

We thus propose a statistical modelling approach to address the lack of tools to analyse dynamic ecological networks [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. Our approach mainly consists in extending one of the techniques dedicated to find structural patterns in static networks, now focusing on their dynamics. We thus first present our methodological proposal, stating the key concepts and introducing the vocabulary required for handling dynamic networks. It is important to stress that our model is suitable for integrating arrivals and departures of entities (corresponding either to species invasion/extinction or birth/arrival/death/departure of individuals) through the possible presence or absence of nodes at the different time steps. Moreover, it can also deal with quantitative edges (and is not restricted to binary interactions) which are often available in datasets.

We first illustrate our approach through the analysis of the dynamic contact network in a colony of ants [START_REF] Mersch | Tracking individuals shows spatial fidelity is a key regulator of ant social organization[END_REF]. Contact networks represent a relevant proxy to study animal sociality [START_REF] Wey | Social network analysis of animal behaviour: a promising tool for the study of sociality[END_REF]. In the literature, these networks may be built from field observations of association between animals (e.g. giraffes in [START_REF] Carter | Fission-fusion dynamics in wild giraffes may be driven by kinship, spatial overlap and individual social preferences[END_REF], trapping data (e.g. field voles in [START_REF] Davis | Spatial analyses of wildlife contact networks[END_REF] and more recently and predominantly from sensorsbased measurements (e.g. song birds in [START_REF] Farine | Interspecific social networks promote information transmission in wild songbirds[END_REF]. These data are now available for large time periods, ranging from days to years of observations for instance. It is therefore possible to investigate the (in-)stability of the social structure [START_REF] Pinter-Wollman | The dynamics of animal social networks: analytical, conceptual, and theoretical advances[END_REF] and potentially question the impact of other time-related factors (seasonal changes, response to stresses such as draught, arrival/departure of a peculiar individual, . . . ).

We lastly present the study of a seasonal trophic network (or food web; [START_REF] Woodward | Quantification and resolution of a complex, size-structured food web[END_REF]. The structure of trophic relations has been intensively studied in the network framework (see [START_REF] Thompson | Food webs: reconciling the structure and function of biodiversity[END_REF], for a review on food webs). Nowadays, following the seminal work of [START_REF] Baird | The seasonal dynamics of the Chesapeake bay ecosystem[END_REF], new datasets allow for monitoring the variation of this structure along temporal gradients (seasons or years), spatial gradients (latitudinal or longitudinal for instance; [START_REF] Kortsch | Climate change alters the structure of arctic marine food webs due to poleward shifts of boreal generalists[END_REF] or qualitative gradients (increasing habitat modification; [START_REF] Tylianakis | Habitat modification alters the structure of tropical host-parasitoid food webs[END_REF]. We will restrict here to dynamic trophic networks corresponding to different temporal snapshots of a food web. In this context, studying such structural variation (or on the contrary, structural stability) can be appropriate to analyse the system's response to major changes (species extinctions, environmental perturbations, climate change, etc). An underlying issue is whether there is resilience of this structural organisation, or rather cyclic dynamics with a return to an equilibrium state.

Materials and methods

From static to dynamic networks

An ecological network is composed of nodes that correspond to any ecological entitiese.g. species, individuals or communities; while edges (or links) characterise presence/absence of an interaction between any two entities and may be valued in some cases. For instance, values may be the frequencies of contacts between two individuals (de [START_REF] De Silva | The dynamics of social networks among female Asian elephants[END_REF] or the number of field observations of interactions between two species. When this network is unique and covers an entire time period, it is called a static network. While many empirical data were aggregated over a whole period of observation recording, it is important to realise that such aggregation could lead to an incorrect understanding of the network structure due to the smoothing aggregation process (cf. Figure 1). An approach to study the temporal dynamics of a set of interactions is the discrete time snapshots approach (see [START_REF] Blonder | Temporal dynamics and network analysis[END_REF], for a complete perspective). It consists in aggregating data over specific time frames (days, months, seasons, years or any relevant frame regarding the ecological system of interest) and to obtain what Blonder et al. call time-aggregated dynamic networks. In the following, we use the term dynamic networks and while we refer to time as being the parameter that drives the evolution, we recall that this could be any other relevant one-dimensional factor.

Formally, we assume T time steps, a number N t of nodes at each time step t, a total number of nodes N (with N N 1 +• • •+N t ) and edges record the presence (possibly valued) or absence of an interaction between any two pair of nodes at each time step. Note that our set-up is different from the one corresponding to the observation of the full interaction flow, namely when data consists in the complete knowledge of edges appearance and disappearance along a continuous gradient. Indeed in our case data is still aggregated over some time intervals or corresponds to a sequence of networks which are specific to a set of discrete values of a one-dimensional factor. When considering continuous time interaction flow data, the object of interest (the flow) is called a temporal network [START_REF] Holme | Modern temporal network theory: a colloquium[END_REF] and this setup will not be explored in this article.

Lastly, it is important to mention that the time frame selection may be an issue in cases where choosing the resolution for the time aggregation is not driven by the ecological question. Indeed, in many cases, the choice of the time frame is expert-based: for instance the dataset from [START_REF] Reynolds | Raccoon contact networks predict seasonal susceptibility to rabies outbreaks and limitations of vaccination[END_REF] consists in T = 52 days of observation including the breeding season, but it is possible to restrict to T = 3 networks (before, during and after the breeding season) to study the network variations due to reproduction period. This choice might not be harmless and for instance [START_REF] Blonder | Temporal dynamics and network analysis[END_REF] showed that the degree distribution in networks can be sensitive to the time frame selection (see also [START_REF] Kun | Locally boosted graph aggregation for community detection[END_REF], for a statistical perspective). It is out of the scope of the present work to explore this frame selection problem.

Stochastic block models (SBM)

In the field of network analysis, one of the most exciting research problems of the last decade has been the network clustering question. Moving beyond descriptive statistics, the goal here is to propose algorithms to extract a high-level view of complex networks, i.e. zooming out the network. Network clustering consists in grouping nodes based on their common characteristics. It often rhymes with finding modules or communities (or compartments; see [START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF]. A module is a set of nodes with much more edges between these nodes than with the others. An important drawback of module-based approaches appears when, quoting Newman and Leicht, we ask: "could there be interesting and relevant structural features of networks that we have failed to find simply because we haven't thought to measure the right thing? ". In other words, is it relevant to search for modular structure in a network that can be structured in any other ways? Following this objection, methods based on statistical inference arose which rely on the principle of grouping nodes that have similar interaction patterns (e.g. hubs, modules, peripheral nodes; see Figure 2) without any a priori knowledge. This is the purpose of a general class of models called stochastic block models (SBM). Network with a complex structure including different patterns: a module (yellow), a hub (red) and a set of peripheral nodes (green). Both networks can be modelled by a SBM with different parameters that capture the structural organisation, with two and three groups respectively. Besides, it is not clear what would be the results of a modular detection algorithm on the second network (as it is not modular).

SBM have been developed for analysing complex networks [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF][START_REF] Newman | Mixture models and exploratory analysis in networks[END_REF][START_REF] Daudin | A mixture model for random graphs[END_REF][START_REF] Guimerà | Missing and spurious interactions and the reconstruction of complex networks[END_REF][START_REF] Matias | Modeling heterogeneity in random graphs through latent space models: a selective review[END_REF] and more recently used to decipher the structure of ecological networks such as hosts-parasites [START_REF] Mariadassou | Uncovering latent structure in valued graphs: a variational approach[END_REF], food webs [START_REF] Picard | Deciphering the connectivity structure of biological networks using MixNet[END_REF][START_REF] Baskerville | Spatial guilds in the Serengeti food web revealed by a Bayesian group model[END_REF] and multi-interactions network [START_REF] Kéfi | How structured is the entangled bank? The surprisingly simple organization of multiplex ecological networks leads to increased persistence and resilience[END_REF].

Let {Y ij } 1≤i,j≤N be the set of edges between any possible couple of nodes (i, j). The model may be defined either for directed or undirected networks and may allow the presence of selfinteractions (edges Y ii ). The principle of SBM is the following: we assume that the ecological entities (nodes) can be gathered into Q groups based on their common interaction properties. Therefore the distribution of Y ij is specified conditionally on the group memberships such that

Y ij ∼ f (Θ ql ) | i ∈ group q, j ∈ group l
where f is any probability distribution parametrised by Θ (called interaction parameter).

The group memberships are unknown, as well as the interaction parameters. An EM-like algorithm (Expectation-Maximisation, see [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] allows for jointly estimating memberships and parameters [START_REF] Daudin | A mixture model for random graphs[END_REF]. The statistical procedure finally displays a high-level view of the network: what kind of interaction patterns are present (through the interaction parameters Θ ql ) and which nodes participate in those patterns (through the group memberships).

A key advantage of SBM is the possibility to plug any probability distribution f in order to fit any kind of interactions. For instance, one can use a Bernoulli distribution for binary interactions [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF], a Gaussian for frequencies or a Poisson for a number of interactions (weighted interactions; [START_REF] Mariadassou | Uncovering latent structure in valued graphs: a variational approach[END_REF], a multinomial for finitely many values or even a multivariate distribution for multivariate edges. One can also use the combination of any of those distributions with a Dirac mass at 0 so as to obtain a 0-inflated distribution that accounts for sparsity in the network (in the valued case, not all interactions necessarily exist; [START_REF] Ambroise | New consistent and asymptotically normal parameter estimates for random-graph mixture models[END_REF]. Relying on a probabilistic framework allows for modelling some randomness and variability in the observations and consequently provides robustness to possible errors or missing data.

Dynamic stochastic block models (dynSBM)

How can we analyse dynamic ecological networks to extract structural information? At the time of writing, only a few alternatives based on descriptive statistics [START_REF] Holme | Modern temporal network theory: a colloquium[END_REF] or on evolving modules [START_REF] Mucha | Community structure in time-dependent, multiscale, and multiplex networks[END_REF] have been considered. Following the above mentioned objections, we claim that a model-based clustering approach could be relevant and we recently proposed to extend the SBM approach to dynamic networks introducing dynamic stochastic block models (dynSBM, [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF].

To develop such an extension, an important question to ask is what could be the meaning of zooming out an object that can change with time? Our answer is to capture the dynamics of a high-level view of the network. This means tracking the evolution of the group behaviours (i.e. the interaction parameters) as well as the nodes group memberships with time. Technically, we rely on a collection of SBM for modelling the different snapshots at each time step combined with N (the number of nodes) independent and identically distributed Markov chains that capture the evolution of a node group through time. Thus at any time step t, our estimate of the group of a node i depends on the SBM estimated for the network at time t and on the group of this node at time step t -1. The model is now characterised through

Y t ij ∼ f (Θ t ql )
|i ∈ group q, j ∈ group l at time t P(i ∈ group q at time t | i ∈ group q at time t -1) = Π qq where Π is the (common) transition matrix of the N different group memberships Markov chains. Reconstructing the different SBM and the common Markov chain parameters has to be done jointly.

As demonstrated in [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF], without adding some constraints this model is not identifiable because of a possible label switching phenomenon between the time steps (which is not the usual label switching encountered in any model with latent groups). To illustrate this phenomenon on a toy example, let us consider a dynamic network where the same static network is observed at two different time steps. We assume that this network is a star (one central node -called the hub -is connected to all others -called peripheral nodes). A SBM is fitted with 2 groups. Supposing group 1 is the hub and group 2 is composed of the peripheral nodes at t = 1, two alternative scenarios are possible at t = 2: either group 1 is still the hub and group 2 is still composed of the peripheral nodes, or the reverse (group 1 is the peripheral nodes and group 2 the hub). Both scenarios are equivalent clusterings when considered at each time step separately (because clusters are defined up to a permutation of their labels only), but globally (meaning when considering the 2 time steps jointly) these clusterings are different and may be fitted by dynSBM with different parameter values that induce different interpretations. Indeed, in the first scenario there are no group-switches while in the second all the nodes undergo a group change. Therefore, we need to add constraints to the dynSBM parameters to remove this ambiguity. In [START_REF] Yang | Detecting communities and their evolutions dynamic social networks -a Bayesian approach[END_REF], the authors proposed to constrain the group memberships to be constant with time (see alternative 1 in Figure 3). But this is clearly not suited to ecological networks where entities can evolve and move from one role to another; e.g. fission-fusion societies [START_REF] Rubenstein | Similar but Different: Dynamic Social Network Analysis Highlights Fundamental Differences between the Fission-Fusion Societies of Two Equid Species, the Onager and Grevy's Zebra[END_REF], modification of animal behaviour between breeding and non-breeding seasons [START_REF] Reynolds | Raccoon contact networks predict seasonal susceptibility to rabies outbreaks and limitations of vaccination[END_REF], response to climate change [START_REF] Kortsch | Climate change alters the structure of arctic marine food webs due to poleward shifts of boreal generalists[END_REF]. We consequently choose to allow the group memberships to vary freely with time but constrain some of the interaction parameters to be stable in time (see alternative 2 in Figure 3). To be more specific, we assume that the intra-group interactions are constant (namely Θ t qq does not depend on t, for any group q). However the interaction parameters between different groups may vary with time (namely Θ t ql depends on t, for any pair of groups (q, l) with q = l). This restriction is sufficient to characterise the groups and solve the identifiability issue [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. Note that this will also be helpful to interpret dynSBM results as the groups are now unambiguously defined and meaningful: each group has its specific connectivity behaviour, with constant intra-group and possibly varying extra-group connectivities. In the meantime, the memberships are free and may change, such that any entity can change its behaviour.

We stress that our approach is different from a naïve one that would separately cluster each network and use an ad-hoc procedure to resolve the label switching problem between two time steps. Indeed, our hidden Markov chain modelling induces dependencies between the networks at different time steps. As a consequence, the clusters recovered for one specific network use information about the others. We also mention that the (maximal) number of groups is fixed with time though some groups might be empty at some time steps. This number may be selected either through a statistical model selection criterion called ICL or relying on heuristic procedures [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF].

To summarise, the dynSBM approach allows for exploring the following questions: i) Is there any high-level structure in the network, i.e. does dynSBM find more than a single group of nodes? ii) Does this network structure vary with time, i.e. are the nodes group memberships evolving with time? iii) What are the group switches trends and frequencies, i.e. what are the values of the underlying Markov chain parameters? iv) Are there any stable or unstable individuals, i.e. are there peculiar group memberships trajectories? Datasets Ants interaction networks. Colonies of the ant Camponotus fellah were followed with a tracking system that monitored the individual positions over days of observations and dynamic social interactions were deduced from physical proximity [START_REF] Mersch | Tracking individuals shows spatial fidelity is a key regulator of ant social organization[END_REF].

Broadstone Stream seasonal food webs. This dataset concerns the aquatic macroinvertebrate community of Broadstone Stream in south-east England [START_REF] Schmid-Araya | The importance of meiofauna in food webs: Evidence from an acid stream[END_REF][START_REF] Woodward | Quantification and resolution of a complex, size-structured food web[END_REF]. Six seasonal connectance food webs were recorded, one every two months from May 1996 to April 1997. We restricted here to simple presence/absence Figure 3: Adding constraints to avoid label switching between time steps. At time t = 2 there are two equivalent alternative scenarios. In scenario 1 group memberships stay constant in time and interaction parameters change with time; in scenario 2 group memberships change between the two time steps but the intra-group interactions stay constant. In this latter case, the green group is characterised by being a 'peripheral nodes' group, the yellow group is a community and the red group is a hub. information on species (nodes) and binary feeding links (edges) and did not consider quantitative data.

Results

In a previous work, we proposed a statistical analysis of two animals contacts networks (sparrows and onagers respectively, see [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. Here, we aim at focusing at a more ecological perspective and specifically explore issues raised by dynSBM on two ecological datasets: a dynamic contact network of ants [START_REF] Mersch | Tracking individuals shows spatial fidelity is a key regulator of ant social organization[END_REF] and a seasonal food web [START_REF] Woodward | Quantification and resolution of a complex, size-structured food web[END_REF].

Dynamic animal contact networks

The data corresponds to a colony of N = 152 Camponotus fellah ants observed during T = 10 days. Edges of the resulting dynamic network are weighted by the number of interactions between each pair of ants and the network is thus undirected with no self-interactions. After examination of the weights distribution, we chose to bin those weights into M = 3 categories corresponding to low, medium and high interaction intensity. We consequently fitted a dynSBM with a multinomial distribution f (in fact as many multinomials as the number of group pairs {q, l}). We selected Q = 3 groups with the heuristic "elbow" method (see [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF].

We first focus on the overall structure of the dynamic network by observing the inter/intragroups interaction properties, as shown in the different cells of Figure 4. Note that the global Q × Q matrix shown here is symmetric as the network under consideration is undirected. The first key concept here is the sparsity level, i.e the amount of edges that are present over all the possible relations (without considering edge values). We clearly see that intra-group interactions are very frequent, in particular in groups 1 and 2 where almost any pairs of ants of these groups are in contact (Figure 4, large blue areas in diagonal plots). This pattern is stable in time (10 days, x-axis in Figure 4). The most interesting trend about inter-groups interaction concerns group 3 which contains ants that interact with those of group 1 but much less with those of group 2 (Figure 4, smaller blue areas in cells (2,3) or (3,2) than in cells (1,3) or (3,1)). These properties are the key factor determining the group boundaries (i.e the memberships) as the other inter-group interactions remain frequent. The next key notion is the intensity level that focuses on the values of present edges and reflects the point to which ants of two given groups are more likely to be in contact with low to high intensity. Interestingly, when two ants of group 2 are in contact (edge is present), they are likely to be in contact with a high intensity/frequency value (Figure 4, large dark blue area in cell (2,2)). On the contrary, even if some contacts exist between these ants and those of group 3, which is already unusual, these sporadic contacts exhibit low intensity (Figure 4, larger light blue area in cell (2,3) or (3,2)). Here, the remaining intra/inter-group intensity levels do not reveal any other interesting pattern (equal proportion of intensity categories). With all these observations, we deduce that group 2 is a so-called module (highly intra-group connected ants) relatively disconnected from group 3 and that group 1 gathers ants "at the interface" i.e. interacting with partners from any of the three groups.

We now investigate whether there are some interesting trends in the turnover of group membership. In other words, we wonder whether all ants have the same propensity to move from one group to another one. We first observe that the global group turnover (i.e. the amount of group switches) is low: 46% of ants never switch group. Moreover, there are no group switches between groups 2 and 3 (Figure 5, no fluxes between these groups over time). This observation, along with the low level of interactions between group 2 and 3 that we discussed before, suggests the existence of a "barrier" between these groups that could be a consequence of space positioning. Indeed, [START_REF] Mersch | Tracking individuals shows spatial fidelity is a key regulator of ant social organization[END_REF] showed that ants were distributed over three social groups (obtained by analysing each daily static network and combining those analyses) with different interaction patterns and that there existed some spatial segregation of the groups. We thus propose to compare our groups obtained with dynSBM (which are evolving with time) and the social groups of [START_REF] Mersch | Tracking individuals shows spatial fidelity is a key regulator of ant social organization[END_REF] (which are fixed with time). Focusing on the ants that stay in the same group at least 8 days over 10 (111 ants over N = 152), we note a quasi-perfect match between Mersch et al. groups and our groups (see Table 1). The modular group 2 corresponds to the foragers of Mersch et al., while the other groups 1 and 3 correspond to the cleaners and nurses respectively. Besides retrieving this functional group, we provide another relevant information: it is now possible to study ants playing different social roles over time, i.e. those that experience group switches at certain time points. Indeed, our dynSBM groups allow to pinpoint these q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 Figure 4: Interaction properties between groups on the ants dataset. Interaction presence and intensity between nodes in any of the Q = 3 groups to the others are represented in a global Q × Q matrix; each cell contains T = 10 time points on the x-axis corresponding to the different time steps. Each square represents four areas: the white area is the proportion of absent edges among all possible pairs of interactions; low to dark blue areas correspond to the proportion of edges (among present ones) with low to high intensity value, respectively. Plot obtained with the connectivity.plot function of the dynsbm package.

interesting individuals that modified their behaviour over time and that can be of peculiar interest for specialists of the Camponotus fellah system.

Broadstone Stream seasonal food webs

The number of sampled species of this aquatic macro-invertebrate community varies seasonally (up to N = 26 in total including 10 predators) as well as the number of directed links. This dataset forms a dynamic trophic network with T = 6 snapshots (May, August, October, December 1996, February, April 1997) and we are interested in addressing the temporal variation in the web structure. Five species were not sampled each month but this situation where nodes are present/absent over time is supported by our model (see supplementary Between two days (t = 1, . . . , 10 on the x-axis), each line is a flux that represents the switch of one or more ants from a group to another group (q = 1, . . . , Q represented on the y-axis).

D1-1 D1-2 D1-3 D2-1 D2-2 D2-3 D3-1 D3-2 D3-3 D4-1 D4-2 D4-3 D5-1 D5-2 D5-3 D6-1 D6-2 D6-3 D7-1 D7-2 D7-3 D8-1 D8-2 D8-3 D9-1 D9-2 D9-3 D10-1 D10-2 D10-3
Here, Dt -q denotes group q from day t. The thickness of each line is proportional to the corresponding counts. Plot obtained with the alluvial.plot function of the dynsbm package.

material of Matias and Miele, 2016, for details). It is also important to mention that selfinteractions (cannibalism) exist for 6 out of 10 predator species. Again, our model allows for this behaviour which might distinguish predators among them and structure the network.

We then fitted a dynSBM with Bernoulli distributions f and we selected Q = 4 groups with the ICL criterion (see [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. The inter/intra-groups interaction properties shown in Figure 6 are not symmetric as we consider directed networks. Therefore, for each pair of groups {q, l} their interaction characteristics are twofold: how often species of group q eat those of group l, and the reverse. As such, group 4 is composed of omnivorous species that eat species of any other groups, but are only eaten by species of their own group (this includes cannibalism). Group 3 has the same properties than group 4 with a significant difference: species of group 3 do not eat those of group 4. We conclude that species from group 3 occupy intermediate positions in food chains whereas those of group 4 are top predators. Indeed, group 4 is mainly composed of the three largest species (the top predators Cordulegaster boltonii, Sialis fuliginosa and Plectrocnemia conspersa) whereas group 3 contains mostly three small species (the larvae of the tanypod midges Macropelopia nebulosa, Trissopelopia longimana and Zavrelimyia barbatipes). Group 2 overall gathers preys that are mostly eaten by predators of groups 3 and 4. Species from group 1 are "hidden" species: they do not eat much, and are not much eaten either. This group is obtained by our statistical procedure due to fewer feeding links from/to other species compared to species of the other groups. It is not coherent from a taxonomic point of view as it gathers a mixture of predators with little activity and secondary preys that we call peripheral species. Lastly, our model can deal with the overall decrease in the number of links after October (see Figure 6, blue area decreasing with time in boxes) which is partly due to the fact that tanypods become less predatory and more detritivorous after autumn [START_REF] Woodward | Quantification and resolution of a complex, size-structured food web[END_REF]. q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 q q q q q q 0 1 Figure 6: Interaction properties between groups on the food web dataset. Same as Figure 4 for Q = 4 groups and T = 6 time steps. In this case, only interaction presence is shown (blue area) as we consider binary edges. Moreover, the Q × Q matrix shows directed interaction from group q (lines) to group l (columns). Now, we explore whether species positions in the food chains (namely being top or intermediate predators, peripheral species or common preys) evolve or stay constant across the seasons. We do not expect a low-level prey to become a top predator, but the group boundaries may change due to seasonal diet variations; for instance, Macropelopia nebulosa eats lots of Nemurella pictetii in August but not in April as this species becomes too large; see [START_REF] Woodward | Quantification and resolution of a complex, size-structured food web[END_REF]. Figure 7 shows that group memberships remain stable before winter, but some changes are observed between October and December. In particular, the tanypod species Macropelopia nebulosa belongs to group 4 and changes to group 3 in winter. Indeed, in summer and autumn only, this species diet is similar to the one of the other members of group 4 (the three competitive top predators, that eat each other) while also being their prey. Still between October and December, the stonefly Siphonoperla torrentium becomes an active predator (with only 1 prey in October and 5 in December) and moves from group 1 to group 3. The prey Prodiamesa olivacea becomes commonly eaten during winter and is consequently integrated into group 2 during this period (and moves back to group 1 in April). Lastly, we observe that Brillia modesta changes from group 2 to group 1 between December and February: this species becomes the exclusive prey of the top predators during this period whereas it is a common prey during the other months. 5 for months t = 1, . . . , 6. Here, M t -q denotes group q on month t (for 1 ≤ q ≤ 4) and special group 0 gathers absent entities at each time step.

Discussion

The inclusion of time in network analysis has been a recent challenge that requires ad-hoc modelling approaches. The success of these approaches has to be measured by their ability to extract substantial additional information that would not be caught by a traditional static network analysis. To this aim, we propose the use of our new dynamic stochastic block model to decipher temporal trends or temporary patterns in dynamic ecological networks.

On the ants interaction network dataset, while the overall group behaviour trends are captured by our model, different individual behaviours are also highlighted. This way, our results can be interpreted at different scales. On the food web dataset, our model underlines a clear trophic organization but also seasonal differences in the prey assemblage. These results require further investigation by experts, but it is interesting to note that our approach can play a key role in extracting unexpected patterns.

Our model is grounded on a rigorous statistical method [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF] and is implemented in an efficient R/C++ package that can handle hundreds to thousands of nodes. It is henceforth one of the very first tools for ecologists facing the recent availability of time-ordered datasets or that would like to explore the evolution of ecological networks with respect to a one-dimensional factor.

Figure 1 :

 1 Figure 1: Same data a) modelled by a two time steps dynamic network (left and centre) or b) aggregated over the whole time period into a static network (on the right). The structure in the static case does not reflect the complexity of the network structure which clearly varies with time. Indeed, edges present at t = 1 and t = 2 are disjoint.

Figure 2

 2 Figure2: a) Network with a clear modular structure with two modules (green and yellow) b) Network with a complex structure including different patterns: a module (yellow), a hub (red) and a set of peripheral nodes (green). Both networks can be modelled by a SBM with different parameters that capture the structural organisation, with two and three groups respectively. Besides, it is not clear what would be the results of a modular detection algorithm on the second network (as it is not modular).

Figure 5 :

 5 Figure 5: Alluvial plot showing the dynamics of the group memberships on the ants dataset.Between two days (t = 1, . . . , 10 on the x-axis), each line is a flux that represents the switch of one or more ants from a group to another group (q = 1, . . . , Q represented on the y-axis).Here, Dt -q denotes group q from day t. The thickness of each line is proportional to the corresponding counts. Plot obtained with the alluvial.plot function of the dynsbm package.

Figure 7 :

 7 Figure7: Alluvial plot showing the dynamics of the group memberships on the food web dataset. Same as Figure5for months t = 1, . . . , 6. Here, M t -q denotes group q on month t (for 1 ≤ q ≤ 4) and special group 0 gathers absent entities at each time step.

Table 1 :

 1 Contingency table between Mersch et al. functional groups and our dynSBM groups (restricted to 75% of ants staying at least 8 over 10 days in the same group).

		cleaners foragers nurses
	dynsSBM group 1	29	1	4
	dynsSBM group 2	2	29	0
	dynsSBM group 3	0	0	42

Acknowledgements

The authors would like to thank Sonia Kéfi for helping us in finding the food web dataset and for her comments on earlier versions of this work.

Data and code accessibility

The ants dataset is available at http://datadryad.org/resource/doi:10.5061/dryad. 8d8h7. The food web dataset is available in plain text in Table 2 of [START_REF] Woodward | Quantification and resolution of a complex, size-structured food web[END_REF]. The R software package dynsbm is available at http://lbbe.univ-lyon1.fr/dynsbm.