Bruno Gaujal

Lazlo Gulyas

Yuri Surdati Mansuri

Eric Thierry Mar

Markov chain analysis of an agent based growth model

Keywords: discrete probabilistic dynamics, face homogeneous Markov chains, Zipf's law dynamique de modèle probabiliste discret, chaîne de Markov, loi de Zipf

In this paper we investigate the asymptotic behavior of a discrete and probabilistic dynamical system which can be described as a growth model where autonomous agents aggregates. The aim of this paper is to give a mathematical analysis of the dynamics. The analysis uses face homogeneous Markov chains and thanks to this study we validate a conjecture set by Laszlo Gulyas and Yuri Mansury concerning a growth model for cities where simulations had shown that the sizes of the cities asymptotically distribute as a Zipf's law. In light of our analysis, we discuss how the emergence of such a Zipf's law could be expected in Gulyas-Mansury' model and in its variants.

Introduction

Many growth phenomena in human or natural multi-site systems exhibit a quite remarkable property, known as Zipf's law, which says that if one calculates the logarithm of the rank and of the site sizes and plot the resulting data in a diagram you will get a remarkable log-linear pattern. In other words, Zipf's law states that the size v of the r'th largest site is a power of its rank: v ∼ Cr -α , where α is a constant greater than or equal to 1 and C a positive constant (in Zipf's original work [START_REF] Zipf | Selected Studies of the Principle of Relative Frequency in Language[END_REF], α was 1).

For instance, it appears that Zipf's law governs many features of the Internet. It has been measured that distributions of the sizes of websites in term of the number of pages they include or the number of links given to other sites follow a Zipf's law [START_REF] Adamic | The small world web[END_REF]. This behavior of the World Wide Web has also been observed on the node degree distribution of the graph underlying the Internet backbone (the physical network) [START_REF] Faloutsos | On power-law relationships of the internet topology[END_REF]. It has lead to the design of growth models in order to explain these features such as growth models by Huberman and Adamic [START_REF] Huberman | Growth dynamics of the world wide web[END_REF] with intuitive and simple assumptions or new random graph growth models by Albert and Barabasi [START_REF] Albert | Statistical mechanics of complex networks[END_REF] which can be seen as agents aggregating with preferential attachments. In economics, some studies suggest that several types of ranking of companies follow Zipf's laws and that some other economic data is also concerned [START_REF] Axtell | Zipf distribution of u.s. firm sizes[END_REF]. Another famous example is the distribution of city populations. After a seminal work of Zipf [START_REF] Zipf | Human Behaviour and the Principle of Least-Effort[END_REF], city sizes growth has been investigated by many authors. One may mention the survey papers by Gabaix [START_REF] Gabaix | Zipf's law and the growth of cities[END_REF][START_REF] Gabaix | Zipf's law for cities: An explanation[END_REF][START_REF] Gabaix | Handbook of Urban and Regional Economics, volume 4, chapter 53 -The Evolution of City Size Distribution[END_REF] which also provide their own explanations of the emergence of Zipf's laws. In another recent paper [START_REF] Tong | Zipf's law for cities: A cross country investigation[END_REF], Soo uses a new data set on 75 countries and shows that the Zipf's law applies (with different parameters according to the countries).

The reader interested by Zipf's laws is referred to the very comprehensive bibliography gathered by Li on http://linkage.rockefeller.edu/wli/zipf/.

However, in many cases, either this remarkable feature remains unexplained or rather complicated models are used to derive it.

During the 2003 Exystence Thematic Institute on Complex Systems [START_REF]discrete and computational aspects of complex systems[END_REF], a simple probabilistic model solely based on individual decisions (autonomous agents) was presented by Laszlo Gulyas. It was introduced as a agent based model to study the growth of city sizes [START_REF] Gulyas | Patterns of firm agglomeration: an autonomous agent based model of city formation[END_REF][START_REF] Mansury | Explaining the zipf-law pattern of city formation: An agent based approach[END_REF]. This model consists in a finite number of sites, which can be seen as the cities, which are ranked by their size, then some agents arrive sequentially and join the cities with the following rule: when arriving, each agent has a random set of eligible cities to join, he then decides to join the largest city proposed to him (and thus increase by one its size). While this process is only based on individual decisions of the agents, Gulyas observed on simulations that the distribution of the city sizes asymptotically behaves like a Zipf's law. He conjectured that its parameter α was equal to 2.

Motivated by this conjecture, our work focusses on the mathematical analysis of a discrete probabilistic dynamics for which Gulyas' model is a special case. As a result, we prove that Gulyas' conjecture is true and that some variants of its model also follow Zipf's law.

As mentioned earlier, a large amount of work has been done to try explaining the recurrent apparition of Zipf's laws when studying city sizes growth. Note that our aim is not to compare all these previous models with Gulyas' new model or to confront its behavior with real data. Our work presents some mathematical tools that prove to be useful in the analysis of such discrete probabilistic models. In particular our proofs are based on a face homogeneous Markov chain description of the dynamics of the system.

In Section 2, we present a general growth model for which we provide an asymptotic analysis. In Section 3, we precisely describe Gulyas' model and we show how it falls into the framework presented before. Applying the preceding asymptotic analysis, we prove that his conjecture about the distribution of sizes is true. We then show that the general growth model apply to variants of Gulyas' model yielding some other Zipf's laws. Section 5 concludes the paper by summing up features that made the analysis work and by presenting some open problems.

2 The growth model and its asymptotic analysis

Presentation of the model

We consider the following growth model in discrete time, where agents sequentially aggregate with several sites. We will consider that these are new inhabitants joining cities. Let C = {c 1 , . . . , c n } be a finite set on n cities. The respective populations in the cities at time t are denoted s 1 (t), . . . , s n (t). The initial population is set to 0 (i.e. s i (0) = 0 for all 1 ≤ i ≤ n).

A set of probabilities {p 1 , . . . , p n } is also given (p

1 + • • • + p n = 1).
At time t, the cities can be ranked by their respective populations, let σ t be one permutation of {1, . . . , n} such that s σt(1) (t) ≥ . . . ≥ s σt(n) (t).

At time t + 1, a new person arrives. This person joins one of the cities with the following probabilities: he joins the i'th largest city c σt(i) with probability p i . In other words, s σt(i) (t + 1) = s σt(i) (t) + 1 with probability p i .

Note that this may induce a new permutation σ t+1 if this arrival has altered the order of the sizes of the cities. Note also that the behavior of this system is not completely well defined when several cities have the same size. Indeed, in such cases, several permutations σ t verify s σt(1) (t) ≥ . . . ≥ s σt(n) (t) and the growth of each city at time t depends on σ t . However it should be clear that if we focuss on the sequence of city sizes s σt(1) (t) ≥ . . . ≥ s σt(n) (t) at time t, it does not depend on the choices of σ 1 , . . . , σ t made at all steps up to t.

Our aim is to study the evolution of the city sizes when t goes to infinity and under the assumption that the sequence p 1 , . . . , p n is strictly positive and decreasing ("people prefer larger cities").

For notation purposes, we also set p n+1 = 0 and we introduce I(x) = max{i | 1/p i ≤ x} that will be used later.

An equivalent model

We consider another representation of the model where each city is replaced by a token and its size is replaced by the position (or slot) of the token over the integer line. If several tokens have the same position, they are piled up on their slot. An example of the old and new representations is given in Figure 2. Note that the selection process for the new model makes one token (the right most among all the picked tokens) move to the right by one slot. The models differ slightly is in the presence of ties. In the new model, if two tokens are in the same position, we decide to select the one on "top", rather than selecting arbitrarily. Similarly when one token moves to the right and joins other tokens at its new position, we choose to insert it at the bottom of the stack. This convention has one main feature which will be used in the following: the relative order between the tokens never changes: the tokens cannot overtake each other. Without loss of generality, we number them in decreasing order: the rightmost token is labeled by 1 and the last token (leftmost) is labeled by n (in a stack of tokens, the labels are decreasing from top to bottom). The position (or value) of token i is denoted by v i .

¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ £ £ £ ¤ ¤ ¤ ¤ ¤ ¤ ¥ ¥ ¥ ¥ ¥ ¥ ¦ ¦ ¦ ¦ ¦ ¦ § § § § § § ¨¨¨¨© 6
It should be obvious that the probabilistic laws of both models are the same. In the remaining, we will rather use the token model where no overtaking can occur because it makes the model easier to describe, however, all the analysis below could have been made using the original model.

Markov chain analysis

The state of the system at step t V (t)

def = (v 1 (t), • • • , v n (t)
) obviously forms a Markov chain with a state space included in N n , since the future evolution only depends on the present state. However, its structure is rather complicated and cannot to be analyzed directly. It should be clear that all the tokens "tend" to move to the right during the evolution. However, to compute the speed of this shift to the right, one needs to look first at second order quantities such as the relative positions of the tokens.

Let g i , 1 ≤ i ≤ n, be the size of the gap between token i and token i + 1. For all

1 ≤ i ≤ n -1, g i = v i -v i+1 and g n = v n
We also denote by h i the smallest label of all tokens piled above token i.

h i = min{j ≤ i, v j = v i } = min{j ≤ i, j k=i g k = 0}.
It should also be obvious that

G = (g 1 , • • • , g n ) is a Markov chain living in N n , where the space can be divised in 2 n zones Z S = {(g 1 , • • • , g n ) ∈ N n |∀i ∈ S, g i = 0}
for all S ⊆ {1, . . . , n} on which the chain is homogeneous. It is a face-homogeneous Markov chain [START_REF] Fayolle | Constructive theory of countable Markov chains[END_REF].

The chains

V = (v 1 , • • • , v n ) and G = (g 1 , • • • , g n ) have
the following transition kernels, respectively denoted by O and P :

At step t + 1, with probability p i , 1 ≤ i ≤ n, (v 1 , • • • , v n ) p i → (v 1 , • • • , v h i + 1, • • • , v n ) (g 1 , • • • , g n ) p i → (g 1 , • • • , g h i -1 -1, g h i + 1, • • • , , g n )
The kernel P of the chain G is displayed in Figure 3 for n = 2.

0 1 2 3 0 1 3 2 g 2 p 1 p 2 g 1 Figure 3: The kernel P for n = 2. Lemma 1. The chain G = (g 1 , • • • , g n ) reaches states where g i > 0 for all 1 ≤ i ≤ n, in finite time, almost surely.
Proof. The proof involves the construction of a new chain

B(t) = (b 1 (t), • • • , b n (t)
) with a bounded state space. As shown below, this new chain will reach states with no null components with probability one and this will imply that G reaches states where g i > 0 for all i, with probability one.

Here is the description of

B(t) = (b 1 , • • • , b n ).
The state space is {0, 1, 2} n . We define a function f i similar to h i but for the chain B.

f i = min{j ≤ i, j k=i b k = 0}.
The transition kernel W of B is the following. For each 1 ≤ i ≤ n,

(b 1 , • • • , b n ) p i → (b 1 , • • • , b f i -1 -1, b f i + 1, • • • , b n ) if f i = i and b f i < 2, (b 1 , • • • , b n ) p i → (b 1 , • • • , b f i -1 -1, b f i , • • • , b n ) if f i = i and b f i = 2, (b 1 , • • • , b n ) p i → (b 1 , • • • , b f i -1 -1, b f i , • • • , b n ) if f i < i. (1)
The finite chain B is displayed in Figure 4. It should be clear that the chain B(t) can reach state (2, • • • , 2) from any starting point, with positive probability. Now, we claim that if in the initial state, B(0) ≤ G(0), component-wise, then at any time t, B(t) ≤ st G(t). This can be shown by induction on t. Assume that B(t) ≤ st G(t), then, by coupling the transitions in both chains and looking at the transition kernel, we will show that B(t + 1) ≤ st G(t + 1).

With probability p i ,

G(t + 1) = (g 1 (t), • • • , g h i -1 (t) -1, g h i (t) + 1, • • • , , g n (t)) B(t + 1) = (b 1 (t), • • • , b f i -1 (t) -1, b f i (t) + 1, • • • , , b n (t)) (2) or (b 1 (t), • • • , b f i -1 (t) -1, b f i (t), • • • , , b n (t)). (3)
First, note that if

f i = h i , then it is straightforward to check that in both cases, if B(t) ≤ st G(t), then B(t + 1) ≤ st G(t + 1).
Since B(t) ≤ st G(t), the only case remaining is f i < h i ≤ i. This implies that b h i -1 (t) = 0 and the transition kernel of B uses case [START_REF] Axtell | Zipf distribution of u.s. firm sizes[END_REF].

Since G(t + 1) remains non-negative by construction, this implies

g h i (t + 1) = g h i (t) + 1 ≥ b h i (t) + 1 ≥ b h i (t + 1)
g h i -1 (t + 1) = g h i -1 (t) -1 ≥ 0 = b h i -1 (t) = b h i -1 (t + 1)
g f i (t + 1) = g f i (t) ≥ b f i (t) = b f i (t + 1) ( if f i < h i -1)
g f i -1 (t + 1) = g f i -1 (t) ≥ b f i -1 (t) ≥ b f i -1 (t + 1).
All the other coordinates are unchanged in G and in B with probability p i . This coupling shows by induction that G(0) ≤ st B(0) implies for all t, G(t) ≤ st B(t). Now,to finish the proof, it is enough to note that since B(t) reaches states where b i > 0 for all i in finite time , then so does G(t) with probability one.

The second step of the analysis of the chain G(t) is to show that there is a drift towards states where the gaps between the tokens grow.

Lemma 2. After some finite time t 0 , the chain G(t) has all its components that remain positive almost surely. Furthermore, for all i,

lim t→∞ g i (t) t = p i -p i+1 .
Proof. Let us first construct a Markov chain S(t) = (s

1 (t), • • • , s n t)) with the following kernel Q. If X has all its components positive , then Q(X, Y ) = P r(G(1) = Y |G(0) = X) = P (X, Y ).
If X has some null coordinates then using Lemma 1, the time

ν(X) = inf{t ≥ 0, s.t.G(t) ≥ (1, • • • , 1), G(0) = X}
is finite almost surely and uniformly bounded. Indeed, for all X with null coordinates, let us denote by

X 2 def = (max(X 1 , 2), • • • (max(X n , 2
)) the state X with all its components truncated at 2. Then

ν(X) ≤ st inf{t ≥ 0, s.t.B(t) ≥ (1, • • • , 1), B(0) = X 2 } ≤ st max Z∈{0,1,2} n inf{t ≥ 0, s.t.B(t) ≥ (1, • • • , 1), B(0) = Z},
which does not depend on X. Then, we set Q(X, Y

) = P r(G(ν) = Y |G(0) = X) = P ν(X) (X, Y ).
Now, let us consider the Lyapounov function min : R n → R, which has bounded average increments for the kernel Q. For all X ∈ N n , if X has some null components

E(min(s 1 (1), • • • , s n (1)) -min(s 1 (0), • • • , s n (0))|(s 1 (0), • • • , s n (0)) = X) ≥ 1,
and if X has no null components

E(min(s 1 (1), • • • , s n (1)) -min(s 1 (0), • • • , s n (0))|(s 1 (0), • • • , s n (0)) = X) ≥ n min i=1 (p i -p i+1 ).
An adequate version of Foster's theorem, (see for example [START_REF] Lamperti | Criteria for the recurrence or transience of stochastic process, I[END_REF]), shows that the chain S is transient and min(s 1 (t), • • • , s n (t)) goes to infinity when t goes to infinity. This implies that there exists a finite time τ such that for all t ≥ τ , s i (t) remains strictly positive for all i. Also, the transition kernels of S and G coincide as soon as X > 0 for all components. Therefore, there exists two finite random variables n 0 and m 0 such that S(t) = G(t + m 0 ) for all t ≥ n 0 . This means that after a finite time t 0 , all the components of G remain positive. This ends the proof the first statement.

As for the second statement, consider (δ j ) j∈N , an i.i.d. sequence such that δ j = 1 with probability p i , δ j = -1 with probability p i+1 , and δ j = 0 with probability 1 -p i+1 -p i .

Since g i (t) is strictly positive almost surely after a finite time t 0 and since for all t ≥ t 0 , g i (t) -g i (t 0 ) = t j=t 0 δ j , the strong law of large numbers yields lim

t→∞ g i (t) t = lim t→∞ t j=1 δ j t = p i -p i+1 a.s.
Using the same method as in the previous lemma, one can use the central limit theorem to get a more precise result on the behavior of g i (t) when t goes to infinity (in distribution).

lim t→∞ |g i (t) -t(p i -p i+1 )| σ 2 √ t = N (0, 1)
Where N (0, 1) is the normal distribution with mean 0 and variance 1 and

σ 2 = p i + p i+1 - p 2 i -p 2 i+1 + 2p i p i+1 .
Let N (t) be the number of tokens whose positions are larger than after t steps.

N (t) = max{i, n j=i g j (t) ≥ }.
Using the central limit theorem,

n j=i g j (t) = p i t + √ tS i (t), (4)
where S i (t) is a centered random variable such that S i (t) converges in distribution to a normal law with zero mean and bounded variance, since the g j (t) are eventually independent variables. If one denotes k := N (t), then by definition of k,

p k t + √ tS k (t) ≥ > p k+1 t + √ tS k+1 (t).
Dividing by t, one gets

p k + S k (t) √ t ≥ t > p k+1 + S k+1 (t) √ t . Now, S k (t)
√ t converges to 0 in probability when t goes to infinity, uniformly in . Finally we get the following result.

Theorem 3. If = p i t + O(t 1/2
) for some i ∈ {0, . . . , N }, then there exists a discrete distribution J whose support is {i -1, i} such that

N (t) ---→ t→∞ J in distribution.
For all the other cases,

N (t) ---→ t→∞ I( /t) in probability, (5)
where I(x) := max{i | p i ≥ x}. Remark 4. Let p i be a Zipf 's law with parameter α. In the above discussion, Equation (4)

shows that when t is large, the sizes of the cities follow a Zipf 's law (with parameter α).

As for Equation (5), it gives the related cumulative distribution function, which is Pareto (with parameter 1/α), i.e. at time t, the probability for city sizes to be larger than is approximatively A -1/α where A = (Ct) 1/α /n since p i ∼ Ci -α and I( /t) = max{i | p i ≥ /t} ∼ ( Ct ) -1/α . The same holds in the first case.

Remark 5. It should also be clear from all the reasoning that initializing the process with city sizes different from zero does not change the asymptotic behavior.

Remark 6. We have used the fact that the sequence (p i ) 1≤i≤n is decreasing to apply Forster's theorem which ensures that from a certain time the gap between city sizes remains strictly positive and then the city of rank i grows almost independently from others at speed p i . If the sequence (p i ) 1≤i≤n is not decreasing, then the asymptotic behavior of the facehomogeneous Markov chain is more difficult to analyze. In general the cities group into clusters, each with the same growing speed. However, given the p i 's, it is an open question to know, for example, if all sites form a single cluster (which is the same as assessing the ergodicity of the corresponding face-homogeneous Markov Chain, see [START_REF] Fayolle | Constructive theory of countable Markov chains[END_REF] for discussions about this).

Application to Gulyas' Model

The next model was proposed by Laszlo Gulyas as an alternative simple model for city growth. Given a set of n cities with initial populations equal to 0, new inhabitants arrive sequentially. When a new inhabitant arrives, a random sample of cities is drawn by first picking at random an integer k between 1 and n (with uniform distribution) and then a random subset of k cities from the n ones (with uniform distribution). The person joins the city with the largest size in the selected subset (if there are several ones, choose one of them arbitrarily) and it increments its size by 1.

Running some simulations, Laszlo Gulyas observed that the distribution of city sizes follows a Zipf's law. He conjectured that its parameter was 2.

We will first show that Gulyas' model is a special case of our general model presented in Section 2, and then apply our analysis to validate Gulyas' conjecture.

The probability to choose one city

With this model, the probability p i that the i'th largest city is chosen and has its size incremented by 1 is equal to:

1 n n-i k=0 n-i k n k+1
, since this city can be chosen only if it belongs to a subset of cardinality k + 1 with 1 ≤ k + 1 ≤ n -i + 1 and once k + 1 is fixed (with probability 1/n), among the n k+1 subsets of cardinality k + 1 only the ones containing the i'th largest city and k cities taken from the n -i smaller ones will lead to its choice (thus it occurs with probability n-i k / n k+1 ).

This sum has a closed and very simple expression that we may compute thanks to formal series. We will use the equality:

k≥0 k + m k x k = 1 (1 -x) m+1
We have:

1 n n-i k=0 n-i k n k+1 = 1 n n-i k=0 (n -i)! n! (k + 1) (n -k -1)! (n -i -k)! = 1 n (n -i)! n! n-i k=0 (k + 1) (n -k -1)! (n -i -k)! = (i -1)! n (n -i)! n! n-i k=0 (k + 1) (n -k -1)! (n -i -k)!(i -1)! = (i -1)! n (n -i)! n! n-i k=0 (k + 1) n -k -1 i -1 ( * )
By denoting [x n ]S(x) the coefficient of x n in the series S(x), the under-braced sum ( * ) is equal to:

n-i+1 k=1 k n -k i -1 = [x n ] k≥1 kx k j≥0 j i -1 x j = [x n ] k≥1 kx k j≥0 j + i -1 i -1 x j+i-1 = [x n ]x i k≥1 kx k-1 1 (1-x) 2 j≥0 j + i -1 j x j 1 (1-x) i = [x n ] x i (1 -x) i+2 = [x n-i ] 1 (1 -x) i+2 = n + 1 i + 1
By reintroducing this value in the initial sum, we obtain for the probability to increment the i th greatest integer:

p i = (i -1)! n (n -i)! n! n + 1 i + 1 = n + 1 n 1 i(i + 1)
.

The sequence of p i 's is strictly positive and decreasing, thus we can apply all the results of Section 2. Since the p i 's follow a Zipf's law with parameter 2 and I(x) = max{i | 1/p i ≤ x} ∼ ((n + 1)x/n) 1/2 , with Remark 4, we prove Gulyas' conjecture: Theorem 7. In Gulyas' growth model, the city sizes asymptotically distribute as a Zipf 's law with parameter 2. More precisely when t tends to infinity, the size of the i'th largest city v i (t) verifies v i (t) ∼ n+1 n t i(i+1) and the number N (t) of cities of size larger or equal to verifies The Zipf's law in the previous model was observed by Gulyas thanks to simulation. For a bench of experimental results, the reader is referred to his work on such models [START_REF] Gulyas | Patterns of firm agglomeration: an autonomous agent based model of city formation[END_REF][START_REF] Mansury | Explaining the zipf-law pattern of city formation: An agent based approach[END_REF].

N (t) ∼ n+1 nt 1/2 -1/2 .

Simulation

We have also run some simulations of the system, mainly to estimate the speed of convergence towards the asymptotic regime given above. Section 2 gives some clues about it, but there are several hidden functions or constants which asymptotically disappear but play a role in the convergence speed. For now, a precise quantification of this speed remains a difficult problem. Figure 5 presents an example of the evolution of the system when n = 100 for a sequence of random draws. At time t = 100, 1000, 10000, 100000 and 1000000, we have plotted the distribution of the 100 city sizes, each city having as coordinates its rank and its size. We have also plotted the ideal curves where the points should be located under the asymptotic regime, namely p i t ≈ t/(i * (i + 1)) for the point of rank i. We have used log-scaled axes to put the data all together and to point out the power law behavior.

It can be clearly seen that the output data match the mathematical analysis: it converges towards the ideal curves which correspond to Zipf's laws with parameter 2.

The simulation uses aliasing techniques to select which city will be joined by the new person in constant time [START_REF] Walker | An efficient method for generating random variables with general distributions[END_REF].

Variations on this agent-based model

A natural way to generalize Gulyas' model is to keep the scheme that a newcomer chooses the largest city proposed to him but to change the probability distribution used to pick the random sample of cities. Suppose that we do not wish that the city names play any role in the random sampling, i.e. we do not wish to impose before starting some correlated growth of cities. Then given the n initial cities c 1 , . . . , c n , for all 1 ≤ k ≤ n, all the subsets of {c 1 , . . . , c n } of cardinality k should be picked with the same probability. Then the system is fully defined by a set of probabilities (α k ) 1≤k≤n , 1≤k≤n α k = 1, where α k is the probability to choose the cardinality k and thus the probability to pick a fixed subset of cardinality k is α k / n k . Gulyas' original model is when you take α k = 1/n for all 1 ≤ k ≤ n.

Given the sequence (α k ) 1≤k≤n , this model is clearly equivalent to the general model of Section 2 where the probability p i to join the i'th largest city is equal to:

n-i k=0 α k+1 n-i k n k+1
. Moreover the sequence (p i ) 1≤i≤n is decreasing if and only if there exists 2 ≤ k 0 ≤ n -1 such that α k 0 > 0. In that case, all the analysis of Section 2 applies. We have seen that asymptotically the size of the i'th largest city is approximatively p i t, thus the emergence of a Zipf's law is directly linked to the form of the p i 's and not to intricate combinatorics due to cities overtaking each other.

Note that formal series may be also used here to compute the values of p i . All the reasoning of Section 3.1 applies and following the same steps, it can be checked that :

p i = (i -1)!(n -i)! n! [x n-i ] S (x) (1 -x) i , (6)
where [x n-i ] is the coefficient of x n-i in the series and S(x) = n k=1 α k x k , S (x) being its derivative.

Non uniform choices for city sets

In the following, we analyse several cases. We have seen above that when the probability to choose a set of cities of size k is uniform (α k = 1/n) then, p i = C n /i 2 , where C n is a normalizing constant. This implies that the long-term city sizes follow Zipf's laws with parameter 2. In the following we consider other values for α k . We focus on two cases. First let α k = a n /k (where a n is a normalizing constant). This means that large sets have a smaller chance to be chosen which could actually be the case in real life where choices are often limited.

The other case is α k = b n n k (again b n is a normalizing constant). This corresponds to a uniform distribution over all subsets of cities, without taking the number of cities into account.

Case (α k = a n /k) If α k = a n /k, then by definition, S(x) = a n n k=1 x k /k and S (x) = a n n-1 k=0 x k-1 = a n 1 -x n 1 -x . Therefore, S (x) (1 -x) i = a n 1 -x n (1 -x) i+1 .
Now, looking at the coefficients,

[x n-i ] S (x) (1 -x) i = a n [x n-i ] 1 -x n (1 -x) i+1 = a n [x n-i ] (1 -x n ) n k=1 k + i i x k = a n n i .
Using [START_REF] Gabaix | Zipf's law for cities: An explanation[END_REF],

p i = 1 n i 1 i n i = a n i . (7)
Therefore, p i follows a Zipf's law with parameter 1. Now, using Theorem 3, the city sizes also follow a Zipf's law with parameter 1. Now, looking at the coefficient of x n-i ,

Case

p i = 1 n i 1 i [x n-i ] S (x) (1 -x) i = b n 2 n 2 i .
In this case, p i does not follow a Zipf's law. It is exponentially distributed as well as the sizes of the cities.

These two cases can be seen as extreme. In general, playing with the probability α k changes the laws of the cities sizes drastically. Actually, it enables one to model many different behaviors of the city sizes, from exponential to power laws with almost arbitrary parameters.

Conclusion

The main feature of this paper is to provide an example of face homogeneous Markov chain that may appear in the modelization of real discrete systems. The mathematical analysis matches simulation results and offers a powerful prediction tool. As shown in the preceding section, the general model we have presented includes some agent-based and multi-site growth models where each site has a value (which may be seen as its power of attraction) and where new agents preferentially choose the more attractive sites among the ones proposed to them (and then increase the attraction of the chosen sites).

Several variations of this type of dynamics may appear in models of real systems. For instance, if the sequence (p i ) 1≤i≤n is not decreasing, the behavior of the asymptotic distribution of sizes is more difficult to describe and the mathematical analysis is much more difficult. Even worse, it might be undecidable to predict some features of the dynamics [START_REF] Gamarnik | On deciding stability of constrained homogeneous random walks and queueing systems[END_REF].

Of course, many other features could be added to the general model described in this paper. One may cite: some synchronicity of the arrivals, some departures from the sites also depending on the site ranks, non-uniform distribution when drawing subsets of the same cardinal e.g. to take into account the geography of the sites and the fact that newcomers will choose their destination among close sites which leads to a correlated evolution of their sizes ... Some of these variations may drastically change the long run behavior and thus the study of each of them constitutes a new interesting challenge.

• The ratio between the sizes of any two cities converges. In particular, if σ i (t) is the i-th largest city at time t, , the ratio with the next largest city exists and can be written

lim t→∞ s σ i (t) (t) s σ i+1 (t) (t) = σ i (∞) σ i+1 (∞)
.

It ranges from 0 to ∞ when σ(i) and σ(i+1) vary from 0 to 1. Its distribution is rather complex and depends on order statistics (same with the laws of σ i (∞) and σ i+1 (∞)) ). Its derivation is not reported here.

• The ratios may depend on which cities are considered It can be shown that the ratios between the sizes of consecuitve cities are smaller for large cities than for small cities (again, this is not reported here). Thefore, the system does not follow a Zip'f law with a constant parameter.

• Furthermore, this system is not ergodic: the limits are random variables, so that one realisation of this system has a very different limit behavior than a second one. Moreover, since i is uniformly distributed, the range of possible behaviors is maximal and equiprobable.

• the non-ergodicity of the system makes it very tricky to control or predict. A simulation will not provide any information on its long term behavior, and the ratio between the cities (Zipf's paramters) can only be computed afterwards. For example, in the simple case of two cities, one very long simulation may provide an asymptotic ratio of σ 1 (∞) = 1 (meaning that almost everyone settles in the largest city), while another very long simulation, starting with the same state may provide a ratio σ 1 (∞) = 1/2 (meaning that only half the population chooses the largest city). Both simulation are valid and equi-probable (as well as any simulation providing any other ratio between 1 and 1/2).

Proof. Let us first prove that the limit exists. This is based on Martingale theory. Let us consider the random variable X i (t) = s i (t) S(t) . Computing its conditional expectation gives E(X i (t + 1)|s Therefore, X i (t) is a martingale with respect to s 1 (t), . . . , s n (t), bounded by 1. Using the Martingale Convergence Theorem [START_REF] Williams | Probability and Martingales[END_REF], this means that there exists a random variable i such that lim t→∞ X i (t) = i a.s.

Next, we show that i is uniformly distributed in [0, 1]. For that, consider the distribution of s 1 (t), . . . , s n (t). We will show by induction that for all x 1 , . . . , x n ≥ 1 such that x 1 + • • • + x n = n + t, P(s 1 (t) = x 1 , . . . , s n (t) = x n ) = 1/K(t), where K(t) = t+n-1 t . At t = 0, x 1 = • • • = x n = 1 is the only case , K(0) = 1 and P(s 1 (t) = x 1 , . . . , s n (t) = x n ) = 1.

Assume that the equation is true at step t and consider the distribution at step t + 1. If x i ≥ 2 for all i, then P(s 1 (t + 1) = x 1 , . . . , s n (t + 1) = x n ) = i P(s 1 (t) = x 1 , . . . , s i (t) = x i -1, . . . , s n (t) = x n ) s i (t) S(t)

= 1 t+n-1 t i x i -1 t + n = 1 t+n-1 t t + 1 t + n = 1 t+n t+1
.

On the other hand, if x i = 1 for some i, then without any loss of generality, one may assume that x 1 = 1, . . . , x k = 1, x k+1 ≥ 2, . . . , x n ≥ 2. .

To end the proof, if P(s 1 (t) = x 1 , . . . , s n (t) = x n ) = 1/K(t) for all t, this means that the distribution of the sizes of the cities is uniform. Taking the limit when t goes to infinity, the distribution remains uniform so that for each city, s i (t)/S(t) converges in distribution to the uniform distribution over [0, 1]. Combining this with the almost sure limit i shown above implies that i is indeed random and uniform over [0, 1].

Actually, a more general theorem can be shown. Theorem 9. Consider a system following Simon's dynamic with α = 0 and n cities, with initial sizes a 1 , . . . , a n . Then, lim t→∞ s i (t) S(t) = i a.s.

where i is a real random variable, with a Dirichlet distribution.

The proof uses the fact that the probablity to go from state a 1 , . . . , a n to state b 1 , . . . , b n has a nice combinatorial structure.

Finally, unlike Simon's model, our model has an ergodic behavior. The ratios between the city sizes converge to a deterministic value which can be evaluated on all simulations and which can also be computed in many cases using the parameters of the model (n, α 1 , . . . , α n ).

Figure 1 :

 1 Figure 1: Once the set of possible choices is fixed ({C,D,E}), the new inhabitant joins the largest city ({C}).

Figure 2 :

 2 Figure 2: The old representation (size on each token) and the new one representation (each token at a position equal to its size).

2 Figure 4 :

 24 Figure 4: The kernel W when n = 2.

Figure 5 :

 5 Figure 5: The sizes of cities depending on their rank at various times during the process, with log-scaled axes.

  (α k = b n n k ) If α k = b n /k, then by definition, S(x) = b n n k=1 n k x k and S (x) = b n n(1 + x) n-1 .

P(s 1

 1 (t + 1) = x 1 , . . . , s n (t + 1) = x n ) = n i=k+1 P(s 1 (t) = x 1 , . . . , s i (t) = x i -1, . . . , s n (t) = x n ) s i (t) S(t)

  1 (t), . . . , s n (t)) =

		s i (t) S(t)	s i (t) + 1 S(t) + 1	+ 1 -	s i (t) S(t)	s i (t) S(t) + 1
	=	s i (t) S(t)
	= X i (t).

A Comparison with Simon's model

In this section, we discuss the relations between our model of city growth and the well known Simon's model.

Let us first give a brief description of Simon's model, which is also agent-based. Consider a system which current state at time t is made of n cities of respective sizes s 1 (t), • • • , s n (t). Each newcomer decides either to create a new city with a given probability α (this new city is then of size s n+1 (t + 1) = 1), or to settle down in one of the existing cities with probability 1 -α. In the latter case, the new comer chooses its home city according to a Bernouilli process with probability s i (t)/S(t), where S(t) = k s k (t).

Simon has shown using a continuous approximation that this model follows a Zipf's law. More precisely that Es σ(i) (t)/s σ(i+1) (t) has a limit (α), when t goes to infinity. This model differs in many ways from ours. The first main difference is the fact that new cities are created during the evolution of the system. However, even when α = 0 (in which case new cities are never created w.p. 1), both systems behave in a very different way.

In the following we will analyse the discrete long-run behavior of Simon's model when α = 0, , which seems to have gone undetected so far. Theorem 8. Consider a system following Simon's dynamic with α = 0 and n cities, all of initial size 1. Then,

where i is a real random variable, uniformly distributed over (0, 1).

Before proving this theorem, several comments can be made.


Pierre Mounoud

Projective properties of Lorentzian surfaces

We investigate projective properties of Lorentzian surfaces. In particular, we prove that if T is a non flat torus, then the index of its isometry group in its projective group is at most two. We also prove that any topologically finite noncompact surface can be endowed with a metric having a non isometric projective transformation of infinite order.

Introduction

Two (pseudo-)Riemannian metrics are called projectively equivalent if they have the same unparametrized geodesics. Projectively equivalent metrics are related to first integrals of geodesic flows. The link is given by the following theorem, already known from Darboux [7, §608] and rediscovered by Matveev and Topalov [START_REF] Matveev | Trajectory equivalence and corresponding integrals[END_REF] (in Russian, see for example [START_REF] Bryant | A solution of a problem of Sophus Lie: Normal forms of two-dimensional metrics admitting two projective vector fields[END_REF] for a recent reference in English).

Theorem 1.1 (Darboux). Two metrics (possibly of different signature) g and ḡ on a surface M are projectively equivalent if and only if the function I : T M → R defined by

I(v) = det(g) det(ḡ) 2 3 ḡ(v, v)
is an integral of the geodesic flow of g.

Let us remark that the integral I given by Theorem 1.1 is quadratic (along fibers). Reciprocally, any non degenerate quadratic integral provides projectively equivalent metrics. It follows therefore from Clairaut Lemma, see Lemma 2.1, that if g has a Killing field then the set of metrics projectively equivalent to g is at least two dimensional. We start this paper by giving direct, but in our opinion interesting, consequences of this fact. For example, we increase the family of Lorentzian tori without conjugate points given in [START_REF] Ch | Sur les surfaces lorentziennes compactes sans points conjugués[END_REF], see Theorem 2.5. We also study the projective extensions of surfaces all of whose spacelike geodesic are closed that have a Killing field, see section 2.3.

The rest of the article, is devoted to non affine projective transformations of Lorentzian surfaces. We recall that a diffeomorphism of a pseudo-Riemannian manifold is called projective if it sends unparametrized geodesics on unparametrized geodesics. The projective transformations of a pseudo-Riemannian manifold (M, g) form a Lie group that contains the group of affine transformations (given by the diffeomorphisms preserving the connection i.e. the parametrized geodesics) that contains the group of isometries. Note that finding non affine projective transformations is the same as finding projectively equivalent Lorentzian metrics with distinct Levi-Civita connections that are isometric.

The projective groups of compact Riemannian manifolds are now well understood, thanks notably to works from Matveev [START_REF] Matveev | Proof of Projective Lichnerowicz-Obata Conjecture[END_REF] and Zeghib [START_REF] Zeghib | On discrete projective transformation groups of Riemannian manifolds[END_REF]. The conclusion of these studies is a theorem à la Ferrand-Obata (compare with [START_REF] Ferrand | The action of conformal transformations on a Riemannian manifold[END_REF]):

Theorem 1.2 ( [START_REF] Zeghib | On discrete projective transformation groups of Riemannian manifolds[END_REF]). Let (M, g) be a compact Riemannian manifold. If M is not a Riemannian finite quotient of a standard sphere, then the projective group of g is a finite extension of its affine group.

It is natural to wonder if such a result still hold in the pseudo-Riemannian setting. In the Lorentzian case, an important step has been done in this direction by Bolsinov, Matveev and Rosemann: Theorem 1.2]). Let (M, g) be the compact Lorentzian manifold. Then any projective vector field on M is an affine vector field.

Theorem 1.3 ([4,
Using the classification of Lorentzian tori whose geodesic flow admits a quadratic integral given by Matveev in [START_REF] Matveev | Pseudo-Riemannian metrics on closed surfaces whose geodesic flows admit nontrivial integrals quadratic in momenta, and proof of the projective Obata conjecture for twodimensional pseudo-Riemannian metrics[END_REF], we prove in section 3.2 the following extension of Theorem 1.3 for surfaces: Theorem 1.4. Let T be a non flat Lorentzian torus. The index of the isometry group of T into the projective group of T is lower or equal to 2.

In section 3.1, we consider, the case of non compact surfaces. We show that in this case the situation is completely different: non affine projective transformations of infinite order often exist.

Theorem 1.5. Every orientable, non compact surface of finite topological type can be endowed with a Lorentzian metric having an affine group of infinite index in its projective group. Moreover, any metric projectively equivalent to it has the same property.

Projective properties of non compact Lorentzian surfaces have already been studied by several authors. We just mention here [START_REF] Bryant | A solution of a problem of Sophus Lie: Normal forms of two-dimensional metrics admitting two projective vector fields[END_REF] by Bryant, Manno and Matveev and [START_REF] Matveev | Two-dimensional metrics admitting precisely one projective vector field[END_REF] (and its appendix) by Matveev. In these papers, the point of view is local (whereas the point of view of Theorem 1.5 is clearly a global one). The authors provide local normal form of metrics admitting projective vector fields. Hence, they work on the disc and the vector fields obtained are not always complete and therefore do not always generate projective transformations.

I wish to thank C. Boubel for making me realize that my first ideas were heading nowhere.

2 Projectively equivalent metrics given by Killing fields

Clairaut Lemma

Let us recall first the very classical Clairaut Lemma:

Lemma 2.1 (Clairaut Lemma). Let (M, g) be a pseudo-Riemannian manifold and K be a Killing field of (M, g). The function

C : T M → R defined by C(v) = g(K, v) is a (linear) integral of the geodesic flow.
Keeping the notations of Lemma 2.1, we note that if C is bounded then for small t ∈ R, the form J = g + tC 2 is a non degenerate quadratic integral of the geodesic flow of g. It follows from Theorem 1.1 that the metric ḡ defined by:

ḡ(v, v) = det(g) det(J) 2 n-1 J(v) (1)
is projectively equivalent to g.

If dim M = 2, on the neighborhood on any point that is not a saddle point of K there exists coordinates (x, y) and a function f such that the metric reads 2dxdy + f (x)dy 2 (see [START_REF] Ch | Extensions maximales et classification des tores lorentziens munis d'un champ de Killing[END_REF] for an account on surfaces with a Killing field). The metric projectively equivalent to g obtained thanks to Clairaut integral are easily given. Proposition 2.2. Let I ⊂ R be an interval and f : I → R be a smooth function. For any (a, ℓ) ∈ R * × R such that 1 + ℓf does not vanish on I, the metrics g and g a,ℓ defined on I × R respectively by 2dxdy + f (x)dy 2 and aℓ (1+ℓf (x)) 2 dx 2 + 2a 1+ℓf (x) dxdy + af (x) 1+ℓf (x) dy 2 have the same unparametrized geodesics.

Proof. Let a, ℓ, z be numbers such that 1 + ℓ z = 0, and

G a,ℓ (z) = a ℓ (1+ℓ z) 2 1 1+ℓ z 1 1+ℓ z z 1+ℓ z , Q a,ℓ (z) = a 2 (1 + ℓ z) 2 1 z z z 2 .
The proof follows directly from Theorem 1.1 and Lemma 2.3 applied to G 1,0 (f (x)).

Lemma 2.3. Setting b = aβ -3 and m = ℓ + µa, we have:

β(G a,ℓ (z) + µ Q a,ℓ (z)) = det(G a,ℓ (z)) det(G b,m (z)) 2 3
G b,m (z).

(

) 2
Proof. It is a direct computation, on one hand we have:

β(G a,ℓ (z) + µ Q a,ℓ (z)) = βa (1 + ℓ z) 2 ℓ + µa 1 + z(ℓ + µa) 1 + z(ℓ + µa) z(1 + z(ℓ + µa))
and on the other hand, as det(G a,ℓ (z)) = -a 2 (1+ℓz) 3 , we have:

det(G a,ℓ (z)) det(G b,m (z)) 2 3 G b,m (z) = a b 4 3 (1 + mz) 2 (1 + lz) 2 b (1 + mz) 2 m 1 + mz (1 + mz) z(1 + mz)
.

When b = aβ -3 and m = ℓ + µa, these two matrices are clearly equal.

Remark 2.4. The distribution orthogonal for g a,ℓ to the Killing field does not depend on a and ℓ.

It means in particular that these metrics have the same "generic reflections". A generic reflection of a surface with a Killing K is a local isometry fixing a non degenerate geodesic perpendicular to K and sending K to -K (as a consequence it permutes the lightlike foliations). In general it is only defined on the saturation of the geodesic by the flow of K, see [START_REF] Ch | Extensions maximales et classification des tores lorentziens munis d'un champ de Killing[END_REF] for details. The metric g a,ℓ has a lightlike geodesic cutting each lightlike orbit of K (what is called a "ribbon" in [START_REF] Ch | Extensions maximales et classification des tores lorentziens munis d'un champ de Killing[END_REF]), it implies that there exits coordinates (u, v) and a function f aℓ such that it reads 2dudv + f a,ℓ (u)dv 2 . According to [START_REF] Ch | Extensions maximales et classification des tores lorentziens munis d'un champ de Killing[END_REF], if K has a lightlike orbit, each of these metrics extends into a bigger surface denoted E u f a,ℓ which is only determined by the function f a,ℓ . These surfaces are still projectively equivalent (some could even be isometric). We can note that it defines an equivalence relation on the set of real functions.

Surfaces without conjugate points

Theorem 1.1 allows us to extend the family of metrics on R 2 inducing tori without conjugate points.

Theorem 2.5. There exists a 2-dimensional family of Lorentzian tori without conjugate points with pairwise non isometric universal cover.

Proof. The main result of [START_REF] Ch | Sur les surfaces lorentziennes compactes sans points conjugués[END_REF] asserts that the Clifton-Pohl plane: (R 2 {0}, 2

x 2 +y 2 dxdy) has no conjugate points. If follows from the discussion above that there exists a two dimensional space of metrics projectively equivalent metrics the Clifton-Pohl plane. Corollary 5.6 of [START_REF] Piccione | Spectral flow, Maslov index and bifurcation of semiRiemannian geodesics[END_REF] implies that having conjugate points is a projective property of Lorentzian surfaces, therefore none of these metrics has conjugate points. The radial field of R 2 is a Killing field of each of these metrics, therefore they all have a quotient diffeomorphic to the torus.

From ( 1), these metrics reads on R 2 {0}:

g a,ℓ = 1 2abxy + a 2 (x 2 + y 2 ) 2 by 2 dx 2 + 2(bxy + a(x 2 + y 2 ))dxdy + bx 2 dy 2 ,
where a, ℓ are real numbers satisfying -a < ℓ < a. We prove now that g a,ℓ is isometric to g a ′ ,ℓ ′ if and only if |a| = |a ′ | and ℓ = ℓ ′ . As the map (x, y) → (-y, x) sends g a,ℓ on g -a,ℓ , we can assume now that a and a ′ are both positive. For any a, ℓ, we have g a,ℓ (K, K) = 0 if and only if xy = 0. Computing the curvature ρ a,ℓ at these points we find -a 2 ℓ. As a > 0, the norm of K is maximal when x = y and ρ a,ℓ (x, x) = -2a 2 (a + ℓ). If there exits an isometry Φ between g a,ℓ and g a ′ ,ℓ ′ , it preserves the set of lightlike orbits of K and the set where K has maximal norm. Consequently we have:

a 2 ℓ = a ′2 ℓ ′ a 2 (a + ℓ) = a ′2 (a ′ + ℓ ′ )
and therefore (a, ℓ) = (a ′ , ℓ ′ ).

Remark 2.6. Each metric g a,ℓ admits a non isometric projective transformation (which is affine if ℓ = 0): the involution (x, y) → (-y, x). It will follow from Theorem 1.4 that the index of the isometry group of the Lorentzian torus associated to g a,ℓ (which is equal to its affine group when ℓ = 0) in its projective group is 2.

Surfaces all of whose spacelike geodesics are closed

Another interesting family of Lorentzian surfaces is given by the set of surfaces all of whose spacelike geodesics are closed. Some of them have also a Killing field, see [START_REF] Mounoud | On spacelike Zoll surfaces with symmetries[END_REF], therefore there exists metrics projectively equivalent to them. However, if two surfaces all of whose spacelike geodesics are closed are projectively equivalent then they have the same lightlike geodesics. Indeed, the lightlike geodesics are the first non closed geodesics of the surface, see [START_REF] Mounoud | Pseudo-Riemannian geodesic foliations by circles[END_REF]. It means they are in the same conformal class, but two conformal metrics that are projectively equivalent are proportional, see [17, Theorem 1.2.3] for example. Hence, the metrics obtained always have a non closed spacelike geodesic. We easily deduce from this fact that a projective transformation of a Lorentzian surface all of whose spacelike geodesics are closed is an isometry.

Let us see what happens in the constant curvature case ie on the de Sitter surface. Let q be quadratic form on R 3 of signature (2, 1). We recall that de Sitter surface is (isometric to) the quadric S q = q -1 (1) endowed with the metric induced by q. Its geodesics are thus the intersections of the vectorial planes with S q . Let q 1 and q 2 be two non degenerate quadratic forms, the signature of q 1 being (2, 1). The radial field of R 3 defines a projection from a subset of S q 1 to S q 2 . In general it is not everywhere defined, but it always sends geodesics of S q 1 on geodesics of S q 2 . When the lightcone of q 2 is contained in the lightcone of q 1 then it is everywhere defined but not surjective. Otherwise said S q 1 is projectively equivalent to a proper open subset of S q 2 . Hence, the de Sitter surface is projectively equivalent to surfaces that do not have all their spacelike geodesics closed but have an extension having this property. Note that we did not assume q 2 to be of Lorentzian type, it could be positive definite. We wonder if any surface projectively equivalent to (an open subset of) a metric all of whose spacelike geodesic are closed can be extended into a metric all of whose spacelike (or timelike) geodesics are closed.

Following [START_REF] Besse | Manifolds all of whose geodesics are closed[END_REF], we give the following definition (we warn the reader that the expression Tannery surface is sometimes used with a slightly different meaning, without reference to a Killing field). Definition 2.7. An pseudo-Riemannian surface is said to be a Tannery surface, if it is orientable, it has a periodic spacelike 1 Killing field, and if all its spacelike geodesics are closed except possibly the one orthogonal to the Killing field (in the Riemannian case).

Lorentzian Tannery surfaces correspond to the surfaces of elliptic type of [START_REF] Mounoud | On spacelike Zoll surfaces with symmetries[END_REF]. Next theorem says in particular that Tannery surfaces have a maximal projective (Tannery) extension. Proof. According to [3, Theorem 4.13], if g 0 is a Tannery surface then there exists global coordinates (r, θ) ∈]0, π[×S 1 on the surface, p q ∈ Q and an odd function h such that g 0 reads:

g 0 = p q + h(cos r) 2 dr 2 + sin 2 r dθ 2 .
It follows from Theorem 1.1, that for any ℓ ∈ R the metric

-1 (1 + ℓ sin 2 r) 2 p q + h(cos r) 2 dr 2 + sin 2 r(1 + ℓ sin 2 r) dθ 2
is projectively equivalent to g 0 on its domain of definition. If we restrict r to ]π/4, 3π/4[ then sin 2 r ∈]1/2, 1[, we can therefore take ℓ = -2. The metric g 1 obtained is Lorentzian. A unitary geodesic γ of g 0 is lightlike for g 1 if and only if 1 -2 sin 4 (r) θ2 = 0.

Moreover, γ is tangent to the parallels {r = π/4} and {r = 3π/4} if and only if g 0 ( γ, ∂ θ ) = ± sin(π/4). But g 0 ( γ, ∂ θ ) = sin 2 (r) θ so these conditions are the same. It means that all the spacelike geodesics of g 1 are closed.

Let us see now that the metric g 1 can be any Lorentzian Tannery surface. Let x : R → R defined by

x(t) =              arcsin -cosh(t) 2 1 -2 cosh(t) 2 , if t ≤ 0 π -arcsin -cosh(t) 2 1 -2 cosh(t) 2 , if t ≥ 0.
It is smooth and made to satisfy sin 2 (a(t))

1-2 sin 2 (a(t)) = -cosh 2 (t). In the coordinates (t, θ), the metric reads:

- p q + (h • v)(sinh t) 2 dt 2 + cosh 2 (t) dθ 2 ,
where v(s) = s √ 1+2s 2 . The function h • v being also odd and v being a diffeomorphism, by a trivial, but rather long, adaptation of chapter 4 of [START_REF] Besse | Manifolds all of whose geodesics are closed[END_REF] (unfortunately it is not the description used in [START_REF] Mounoud | On spacelike Zoll surfaces with symmetries[END_REF]), we see that any Lorentzian Tannery surface can be written this way.

One half of Theorem 2.8 is easily generalizable to Lorentzian surfaces with an other type of Killing field. Proposition 2.9. Let (M, g) be an orientable Lorentzian surface with a Killing K. If the spacelike geodesics of g are all closed, then (M, g) contains a proper open subset projectively equivalent to a Lorentzian surface all of whose spacelike geodesics are closed.

Proof. It is proven in [15, Propositions 3.2 and 3.4] that sup g(K, K) = +∞ and that any spacelike geodesic is somewhere tangent to K. Let ℓ be a positive number and U = {p ∈ M | g(K, K) < ℓ 2 }, it follows from [15, Proposition 3.4] that U is connected. The boundary of U is given by a union of orbits of K. A unitary geodesic γ is tangent to ∂U if and only if C( γ) = ℓ, where C is the Clairaut integral associated to K. From Proposition 2.2, we deduce that J = g -1 ℓ 2 C 2 is non degenerate on U . We endow U with the metric ḡ = det(g) det(J)

2 3
J, it is projectively equivalent to g. A vector v such that g(v, v) = 1 satisfies J(v) = 0 if and only if C(v) = ℓ. It implies as above that the spacelike geodesics of (U, ḡ) are all closed.

Remark 2.10. Unfortunately, Proposition 2.9 cannot be used to deform one of the spacelike Zoll surfaces given in [START_REF] Mounoud | On spacelike Zoll surfaces with symmetries[END_REF] into a spacelike Zoll surface of a different kind. Indeed, it is not difficult to see that the families of surfaces constructed there are precisely the surfaces having the same Clairaut integral as a metric of positive constant curvature. Therefore the metrics obtained have the same Clairaut integral as a metric obtained by deformation of a metric of positive constant curvature. But these metrics also have positive constant curvature and therefore the metrics obtained by deformation are also of the kind described in [START_REF] Mounoud | On spacelike Zoll surfaces with symmetries[END_REF].

3 Projective transformations 3.1 Non compact surfaces, proof of Theorem 1.5

We start by proving the theorem on R 2 . Let f be 1-periodic non negative and non constant function and let m be its maximum. Let a be a real number greater than 1, let ε be a positive number satisfying ε < (a 3 -1)/ma 3 . Let α : [0, 1] → [1, a] be a smooth function equal to 1 on a neighborhood of 0 and equal to a on a neighborhood of 1. Let λ : [0, 1] → [0, 1] be a smooth function equal to 0 on a neighborhood of 0 and equal to ε on a neighborhood of 1.

For any n ∈ Z, we define two functions A n and Λ n on [n, n + 1] by

A n (x + n) = α(x)a n Λ n (x + n) = λ(x) + ε 1 -a 3n
1 -a 3 α 3 (x) Let A (respectively Λ) be the function such that for any n ∈ Z the restriction of A (resp. Λ) to [n, n + 1] is equal to A n (resp. Λ n ). The functions A and Λ are smooth. We see also that for any x ∈ R, Λ(x) > -1 m and therefore 1 + Λ(x)f (x) > 0. Hence the tensor defined by

g = A 3 (x) Λ(x) (1 + Λ(x)f (x)) 2 dx 2 + 2 1 + Λ(x)f (x) dxdy + f (x) 1 + Λ(x)f (x) dy 2
is a Lorentzian metric on R 2 .

When ℓ = Λ(x), a = A 3 (x) and z = f (x), the matrix G a,ℓ (z) defined in section 2 is equal to G x the matrix of g at any point (x, y) and Q a,ℓ (z) is the matrix of the square of Clairaut integral associated to the Killing field ∂ y .

Let x be a point of [0, 1]. Notice that for any n ∈ Z, A(x

+ n) = a 3n A(x) and Λ(x + n) = Λ(x) + ε 1-a 3n 1-a 3 A 3 (x). Therefore relation (2) reads a -n (G x + ε 1 -a 3n 1 -a 3 C 2 (x)) = det(G x ) det(G x+n ) 2 3 G x+n ,
where C is the Clairaut integral. Thus Theorem 1.1 implies that the map τ : (x, y) → (x + 1, y) is a projective map. It is clearly not an isometry. As g is not flat the only parallel endomorphisms of the tangent bundle are the multiples of the identity (see Fact page 2226 of [START_REF] Boubel | Affine transformations and parallel lightlike vector fields on Lorentzian 3-manifolds[END_REF] for a proof). It means that any metric having the same Levi-Civita connection as g is proportional to g. However, the vector field ∂ x is lightlike when x = 0 but not when x = 1, so τ * g is not proportional to g. Consequently τ is not an affine map. If ḡ is a metric projectively equivalent to g then Theorem 1.1 says that it defines an integral I of the geodesic flow of g. If I were not a linear combination of the energy and Clairaut integral, then by a straightforward adaptation of Lemma 2 of [START_REF] Bryant | A solution of a problem of Sophus Lie: Normal forms of two-dimensional metrics admitting two projective vector fields[END_REF], the dimension of the space of projective vector fields of g would be at least 2. But Theorem 1 of [START_REF] Bryant | A solution of a problem of Sophus Lie: Normal forms of two-dimensional metrics admitting two projective vector fields[END_REF] gives a local normal form for such metrics. It is clearly possible to find a function f such that the metric g is not isometric to one of these metrics (in fact the opposite seems to be impossible). Consequently, τ does not preserve ḡ.

We choose now f such that f -1 (0) = Z. We choose also a maximal geodesic γ 0 in ]0, 1[×R perpendicular to ∂ y . For any n ∈ Z, we denote by γ n the image of γ 0 by τ n . Each curve γ n is perpendicular to ∂ y , see Remark 2.4. Let σ n be the generic reflection associated to γ n , ie the isometry of ]n, n + 1[×R sending ∂ y to -∂ y that fixes γ n . We clearly have

τ • σ n = σ n+1 • τ. ( 3
)
We consider then the Lorentzian manifold X obtained by gluing two copies of (R 2 , g) along each set ]n, n + 1[×R thanks to σ n (the fact that X is a Hausdorff Lorentzian surface is proven in [2, Remarque 3.4]). It follows from (3) that τ induces a well defined projective transformation of X.

Finally, Proposition 4.31 of [START_REF] Ch | Extensions maximales et classification des tores lorentziens munis d'un champ de Killing[END_REF] tells us that any non compact, orientable surface of finite topological type can be realized as a quotient of the universal cover of X. Equation (3) tells us that the lift of τ belongs to the normalizer of the group of isometries used to construct these surfaces (the so-called "generic subgroup" in [START_REF] Ch | Extensions maximales et classification des tores lorentziens munis d'un champ de Killing[END_REF]), therefore τ induces a projective transformation on them.

Remark 3.1. We denote by Γ the group generated by τ and the map (x, y) → (x, y + 1). This group preserves the projective connection of (R 2 , g) (the metric defined in the proof of Theorem 1.5). Hence the torus R 2 /Γ is a torus endowed with a projective connection that is locally but not globally Lorentzian. The function g(∂ y , ∂ y ) tends always to 0 when x goes to -∞, therefore, at least when the geodesic flow of g does not admit another first integral (what is probably always the case), there is no Riemannian metric projectively equivalent to g.

Compact surfaces, proof of Theorem 1.4

The starting point of the proof of Theorem 1.4 is the fact proven by Matveev see [12, Theorem 2] that a compact Lorentzian surface whose geodesic flow admits a quadratic integral either admits a Killing field or is of Liouville type (see below for a definition).

Let us assume that there exits a non flat Lorentzian torus T whose isometry group has an index greater than two in its projective group. Let us assume first that it has also a non trivial Killing field K. Lemma 3.2. Let g be a non flat Lorentzian metric on T that possess a non trivial Killing field K. If the index of Isom(T) in Proj(T) is greater than 2, then g(K, K) is non negative or non positive and vanishes and Φ * K = ±K for any Φ ∈ Proj(T).

Proof. Let Φ be a projective transformation of T. As g is not flat, it follows from [START_REF] Matveev | Pseudo-Riemannian metrics on closed surfaces whose geodesic flows admit nontrivial integrals quadratic in momenta, and proof of the projective Obata conjecture for twodimensional pseudo-Riemannian metrics[END_REF]Theorem 6] that the space of projective vector fields of g is spanned by K. As Φ * K is a projective vector field (its flow is a composition of projective transformations), there exists c ∈ R * such that Φ * K = cK. The flow of K being periodic (see for example [START_REF] Ch | Extensions maximales et classification des tores lorentziens munis d'un champ de Killing[END_REF]), |c| has to be 1.

The metric Φ * g is projectively equivalent to g therefore, by Theorem 1.1,

I = det(g) det(Φ * g) 2/3 Φ * g (4)
is an integral of the geodesic flow of g. As g is not flat and has a non trivial Killing field, it follows from [12, Corollary 1] that I is a linear combination of g and the square of Clairaut integral. Hence denoting by C the Clairaut integral there exist a = 0 and ℓ such that I = a -1/3 (g + ℓC 2 ). As Φ preserves ±K, the minimum (resp. maximum) of the functions g(K, K) and Φ * g(K, K) are equal. We denote them by m -and m + . We denote by vol g the volume form of g. The 1-form i K vol g is closed and induces a nowhere zero closed 1-form on the space of leaves of K (that is diffeomorphic to the circle). We choose a coordinate x on this space such that dx = i K vol g . We call such a coordinate a canonical coordinate associated to (g, K). Note that the coordinate x that appears in Proposition 2.2 induces a canonical coordinate associated to (g, K) on the space of leaves. Let f : S 1 → R be the expression of the function g(K, K) in this coordinate. It follows from Proposition 2.2 that the expression of Φ * g(K, K) in this coordinate is af 1+ℓf , where a and ℓ are the numbers given above (note that x is not a canonical coordinate associated to (Φ * g, K)). It is therefore clear that the set of lightlike orbits of K is preserved by Φ and that

if a > 0, Φ(f -1 (m + )) = f -1 (m + ) and Φ(f -1 (m -)) = f -1 (m -) if a < 0, Φ(f -1 (m + )) = f -1 (m -) and Φ(f -1 (m -)) = f -1 (m + ).
If the index of the isometry group in the projective group is greater than 2, we can choose Φ such that leaves invariant both f -1 (m -) and f -1 (m + ). It implies that m ± = am ± 1+ℓm ± , and as m -= m + (because the metric is not flat), m + m -= 0 and a = 1 + ℓm + or 1 + ℓm -.

Let Φ be a projective transformation of T such that Φ * g = ±g. Let p be a point of T belonging to a lightlike orbit of K and let q be the image of p by Φ. There exist a real function f , an interval I (resp. I ′ ) and a neighborhood of p (resp of q) that is isometric to I × S 1 (resp. I ′ × S 1 ) endowed with the metric 2dxdy + f (x)dy 2 . By Lemma 3.2, we can assume that f is non negative.

According to Proposition 2.2 and Lemma 3.2, there exists ℓ ∈ R such that Φ * g reads on each of these neighborhoods

(1 + ℓm + ) ℓ (1 + ℓf (x)) 2 dx 2 + 2 1 + ℓf (x) dxdy + f (x) 1 + ℓf (x) dy 2 ,
note that ℓ = 0 as Φ * g = ±g.

Let J be the Jacobian of Φ at p in the coordinates (x, y). As the orbit of K through q is also lightlike, it satisfies:

t J 0 1 1 0 J = (1 + ℓm + ) ℓ 1 1 0 therefore, knowing that J 0 1 = ± 0 1 , J = ± 1 + ℓm + 0 ℓ/2 1
. The matrix J -1 being the Jacobian of Φ -1 at q, the matrix of Φ -1 * g at q is given by:

t J -1 0 1 1 0 J -1 = 1 1 + ℓm + -ℓ 1+ℓm + 1 1 0 .
As Φ -1 * g is also projectively equivalent to g, by Proposition 2.2 there exist a and c such that Φ -1 * g reads

a c (1 + cf (x)) 2 dx 2 + 2 1 + cf (x) dxdy + f (x) 1 + cf (x) dy 2 .
Clearly a = 1 1+ℓm + and c = -ℓ 1+ℓm + . Therefore, possibly replacing Φ by Φ -1 , we can assume that ℓ < 0 (as ℓ > 0 implies 1 + ℓm + > 0).

We now compare the values of the integrals of the forms i K vol g and i K vol Φ * g on the space of leaves of K. The transformation Φ preserving K, these values must be the same. If we denote by dx the form i K vol g then the form i

K vol Φ * g is 1+ℓm + (1+ℓf (x)) 3 2
dx, where f is the expression of g(K, K) in the coordinate x.

As ℓ < 0 and f ≥ 0, for any x we have

0 < (1 + ℓf (x)) 3/2 < 1 + ℓf (x) < 1 + ℓm + < 1,
therefore for any x, we have 1+ℓm + (1+ℓf (x))

3 2

> 1 and the two integrals are therefore different. We have a contradiction, therefore T has no Killing field.

If (T, g) has no Killing field then by [START_REF] Matveev | Pseudo-Riemannian metrics on closed surfaces whose geodesic flows admit nontrivial integrals quadratic in momenta, and proof of the projective Obata conjecture for twodimensional pseudo-Riemannian metrics[END_REF]Theorem 2], it is of Liouville type, therefore by [12, Lemma 1] it is projectively equivalent to a Riemannian torus (T, g 0 ) also of Liouville type. We recall that a metric is of Liouville type if there exists two real periodic functions h 1 and h 2 and a lattice Λ such that it is isometric to the quotient of R 2 endowed with the metric (h 1 (x)+h 2 (y))(dx 2 ±dy 2 ) by Λ. In the Riemannian case, the function I 0 : T T → R defined by I 0 (v) = (h 1 (x)+h 2 (y))(h 1 (y)dx 2h 2 (x)dy 2 ) is an integral of the geodesic flow of g 0 , cf. [START_REF] Kiyohara | Compact Liouville surfaces[END_REF] for example (the expression in [START_REF] Matveev | Pseudo-Riemannian metrics on closed surfaces whose geodesic flows admit nontrivial integrals quadratic in momenta, and proof of the projective Obata conjecture for twodimensional pseudo-Riemannian metrics[END_REF] looks different because the integral given there is on the cotangent bundle).

As g 0 is not flat, by [12, Corollary 1], the set of quadratic integrals of the geodesic flow of g 0 is 2-dimensional and therefore spanned by I 0 and the energy. It follows by Theorem 1.1 that the set P 0 of Riemannian metrics projectively equivalent to g 0 is diffeomorphic to {(a, b) ∈ R 2 | ah 2 (y)+b > 0 and -ah 1 (x) + b > 0, ∀(x, y) ∈ R 2 } and therefore to R 2 . Clearly, the identity component of the projective group of g 0 acts trivially on this set.

According to [START_REF] Zeghib | On discrete projective transformation groups of Riemannian manifolds[END_REF], the group of isometry of g 0 has finite index in the group of projective transformations. It implies that the group of projective transformations has a finite number of connected components. Hence the action of the projective group on P 0 reduces to the action of a finite group. As this action is smooth it preserves a Riemannian metric. By the uniformization theorem, this action is conjugated to an action by biholomorphisms on C or the disk. Consequently, the action of the projective group fixes a metric of P 0 , we can assume that it is g 0 . The set of non isometric projective transformations of g is contained in the set of isometries of g 0 that do not preserve I 0 . Theorem 1.4 therefore follows from Lemma 3.3. Lemma 3.3. Let g 0 = (h 1 (x) + h 2 (y))(dx 2 + dy 2 ) be a Riemannian metric of Liouville type on the torus. The group of isometries of g 0 preserving the integral I 0 as an index lower or equal to 2 in the isometry group of g 0 . More precisely, there exists an isometry Φ of g 0 that does not preserve the Liouville integral I 0 if and only if there exist k ∈ R and c ∈ R such that h 1 is 2k-periodic and h 2 (x + k) = h 1 (x) + c or h 2 (-x -k) = h 1 (x) + c. Moreover in that case (after a possible change of origin) Φ(x, y) = (y + k, x + k) or Φ(x, y) = (-y + k, -x -k).

Proof. An isometry Φ of g 0 is a conformal transformation of the flat metric induced by dx 2 + dy 2 and therefore an isometry of it. It implies that the function h 1 + h 2 is Φ invariant. As h 1 + h 2 is not constant, Φ preserves or permutes the lines R∂ x and R∂ y , therefore Φ 0 , the linear part of Φ in the coordinates (x, y), is one of the following: ±Id, ± 0 1 1 0 , ± 0 1 -1 0 , ± 1 0 0 -1 .

If Φ 0 = ± 0 1 1 0 then, by a short computation, Φ preserves h 1 and h 2 (and not just their sum)

and therefore I 0 . It implies that the product of two isometries always preserves I 0 proving the first part of the lemma.

If Φ 0 = 0 1 1 0 then (after a possible change of origin) there exists k ∈ R such Φ(x, y) = (y + k, x + k). We therefore have h 1 (x) + h 2 (y) = h 1 (y + k) + h 2 (x + k), for any x and y. It follows that h 1 (x) -h 2 (x + k) is constant and h 1 (x + 2k) = h 1 (x). Reciprocally, if h 1 (x) -h 2 (x + k) is constant and h 1 (x + 2k) = h 1 (x) then the map (x, y) → (y + k, x + k) is an isometry that does not preserve the Liouville integral.

If Φ 0 = -0 1 1 0 then (after a possible change of origin) there exists k ∈ R such Φ(x, y) = (-y + k, -x -k). We therefore have h 1 (x) + h 2 (y) = h 1 (-y + k) + h 2 (-x -k), for any x and y. It follows that h 1 (x)-h 2 (-x-k) is constant and h 1 (x+2k) = h 1 (x). Reciprocally, if h 1 (x)-h 2 (-x-k) is constant and h 1 (x + 2k) = h 1 (x) then the map (x, y) → (-y + k, -x -k) is an isometry that does not preserve the Liouville integral.

Theorem 2 . 8 .

 28 Any Lorentzian Tannery surface is projectively equivalent to an open subset of a Riemannian Tannery surface. Reciprocally any Riemannian Tannery surface contains a subset projectively equivalent to a Lorentzian Tannery surface.

in the Riemannian case, we call any non zero vector spacelike