A dominance tree approach to systems of cities Thomas Louail * CNRS, Géographie-cités, Campus Condorcet, 93322 Aubervilliers Cedex, France Marc Barthelemy † Institut de Physique Théorique, CEA, CNRS-URA 2306, F-91191, Gif-sur-Yvette, France and CAMS (CNRS/EHESS) 54 Avenue de Raspail, 75006 Paris, France

Characterizing the spatial organization of urban systems is a challenge which points to the more general problem of describing marked point processes in spatial statistics. We propose a nonparametric method that goes beyond standard tools of point pattern analysis and which is based on a mapping between the points and a 'dominance tree', constructed from a recursive analysis of their Voronoi tessellation. Using toy models, we show that the height of a node in this tree encodes both its mark and the structure of its neighborhood, reflecting its importance in the system. We use historical population data in France (1876-2018) and the US (1880-2010) and show that the method highlights multiscale urban dynamics experienced by these countries. These include nonmonotonous city trajectories in the US, as revealed by the evolution of their height in the tree. We show that the height of a city in the tree is less sensitive to different statistical definitions of cities than its rank in the urban hierarchy. The method also captures the attraction basins of cities at successive scales, and while in both countries these basin sizes become more homogeneous at larger scales, they are also more heterogeneous in France than in the US. Finally, we introduce a simple graphical representation -the height clock -that monitors the evolution of the role of each city in its country.

Understanding the organization of urban systems has always been a central challenge in geography and economics [START_REF] Helsley | Matching and agglomeration economies in a system of cities[END_REF][START_REF] Sanders | SIMPOP: a multiagent system for the study of urbanism[END_REF][START_REF] Rozenblat | Firm linkages, innovation and the evolution of urban systems[END_REF]. On one hand the population sizes of cities have been extensively discussed since Zipf's work [START_REF] Zipf | Human behavior and the principle of least effort: An introduction to human ecology[END_REF][START_REF] Gabaix | Zipf's law for cities: an explanation[END_REF][START_REF] Verbavatz | The growth equation of cities[END_REF], and have been shown to follow a broad distribution: there is a hierarchy of cities characterized by many small cities, a few medium cities and a very few large ones, whose population sizes are much larger than the rest. While this distribution has been fully characterized and discussed [START_REF] Verbavatz | The growth equation of cities[END_REF][START_REF] Cottineau | MetaZipf. A dynamic metaanalysis of city size distributions[END_REF], the spatial distribution of cities has also been an important subject of debate [START_REF] Christaller | Die zentralen Orte in Suddeutschland: Eine okonomisch-geographische Untersuchung uber die Gesetzmafligkeit der Verbreitung der Siedlungen mit stddtischen Funktionen PhD thesis[END_REF][START_REF] Ullman | A theory of location for cities[END_REF][START_REF] Haggett | Locational analysis in human geography[END_REF][START_REF] González-Val | The Spatial Distribution of US Cities[END_REF].

More generally, this problem of characterizing the spatial organization of urban systems appears often in spatial statistics [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF][START_REF] Jacobsen | Point process theory and applications: marked point and piecewise deterministic processes[END_REF] where the points (cities in our case) have a position and are described by (at least) one quantity such as the population. In spatial statistics this general setting is referred to as a 'marked point process' [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF], where each point x i of the process carries extra information called a mark m i . The mark can be a random or a categorical variable, or any other additional information about the points. This problem is not only important in geography with the study of human settlements, but is also relevant for many fields ranging from ecology (positions of plants of different species), to epidemiology (locations of infected individuals), material science (positions of defects) or astronomy that is interested in the location of stars and galaxies. An urban system is an example of such marked point process, where the points are the cities and the mark is their population size. Space is obviously important in this system, and has to be considered jointly with population: if for example two cities are close to each other, it makes a difference if they have similar or very different population sizes. Standard tools developed for the analysis of point processes usually consist in measuring spatial autocorrelations [START_REF] Moran | Notes on continuous stochastic phenomena[END_REF], or testing the null hypothesis of complete spatial randomness with the K or L statistics which summarize the deviations from a uniform (Poisson) distribution [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF][START_REF] Jacobsen | Point process theory and applications: marked point and piecewise deterministic processes[END_REF]. These tools were then extended to marked point processes and describe deviations from known cases such as the Poisson process, or intensity and moments measures [START_REF] Haenggi | Stochastic geometry for wireless networks[END_REF]. A tool that would go beyond these measures and provide a more precise characterization of these processes would then be extremely useful for a wealth of different problems.

To further characterize marked point processes it is possible to start from geometrical structures constructed on top of the point pattern. These 'secondary structures' [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF] comprise in particular tessellations and networks [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF]. The Voronoi tessellation is one of the most relevant structure in computational geometry, and is of major importance in the resolution of many problems, in particular in location science [START_REF]Location Science[END_REF]. Networks can also be constructed over a set of points, and useful tools include proximity graphs (such as the random geometric graph) or excluded volume graphs (such as the Gabriel graph). Measures on these secondary structures can then characterize the point process itself, and constructing a spatial network on top of the point process enables to import all the networks knowledge into spatial statistics. For example, a recent approach uses first-passage times of random walks on networks constructed over a set of points in order to quantify correlations in complex systems [START_REF] Bassolas | First-passage times to quantify and compare structural correlations and heterogeneity in complex systems[END_REF]. However, approaches based on secondary structures are in general used to estimate deviations from uniformity, or the importance of some correlations. In contrast, we want to construct a tool that helps us understand the multiscale structure of a system of marked points, and enable us to compare it to other systems. Such a tool for characterizing the organization of spatial hierarchies should thus encode both the spatial information and the population (in more general terms both the positions and the marks of the points). The purpose of this tool is not to describe the deviation from a uniform distribution, but to compare two systems (such as urban systems in different countries for example) and to understand the temporal evolution of these systems from a non-local, high-level perspective. We will illustrate the construction of such a tool by considering the case of urban systems.

Cities are not scattered at random but follow some logic based on geographical constraints, economical considerations and historical path dependency, and central place theory [START_REF] Christaller | Die zentralen Orte in Suddeutschland: Eine okonomisch-geographische Untersuchung uber die Gesetzmafligkeit der Verbreitung der Siedlungen mit stddtischen Funktionen PhD thesis[END_REF] has seeked to explain the spatial distribution of cities of different sizes based on the idea that settlements function as 'central places' that provide services to surrounding areas. The result of consumers' preferences is then a system of centers of various sizes, forming different levels of a hierarchy. A consequence obtained by Christaller [START_REF] Christaller | Die zentralen Orte in Suddeutschland: Eine okonomisch-geographische Untersuchung uber die Gesetzmafligkeit der Verbreitung der Siedlungen mit stddtischen Funktionen PhD thesis[END_REF] is that the most efficient pattern to serve areas without overlap is a triangular or hexagonal arrangement of settlements. Although this idea has been very successful and inspirational to many scientists, few works have tried to validate it quantitatively. For example, cities have been studied from the correlation point of view [START_REF] Glass | Uniform distribution of objects in a homogeneous field: Cities on a plain[END_REF][START_REF] Hernando | Space-time correlations in urban sprawl[END_REF], but an important contribution to this problem is due to Okabe and Sadahiro [START_REF] Okabe | An illusion of spatial hierarchy: spatial hierarchy in a random configuration[END_REF] who showed that random (uniform) arrangements of cities could explain Christaller's findings. To reach this result they had to define a quantitative tool that captures the spatial dominance of a city on another. They used a tree representation of the dominance relation between marked points, and they used it to characterize both a Poisson point pattern and a real-world case study -market places in Nishinomiya, near Osaka. They analyzed this tree in the perspective of testing Christaller's ideas and in particular measured the ratio K between the number of nodes in two consecutive hierarchical levels. However they did not go further in the characterization of the importance of each point in the system, and in the following we will adapt and extend this method in order to construct a general tool able to characterize a system of marked points. First we will describe the dominance tree introduced in [START_REF] Okabe | An illusion of spatial hierarchy: spatial hierarchy in a random configuration[END_REF], and that will constitute the starting point of our analysis. We will then introduce different measures to extract information from this tree. In particular, we will show that the height of a node in this tree encodes both its mark and space-related information. For systems of cities, this means that the height of city in the dominance tree encodes both its population rank and its location in space. We will first illustrate this method on toy models, and then on empirical data, to compare the evolution of the French and the US urban systems over the 20th century.

Constructing the dominance tree

The method we propose can be applied to any marked point process, and we will illustrate it on the case of cities in a country. We assume that a given country has N cities and each city i is characterized by its location x i in some coordinate system, and its population P i (t) which can vary over time. In order to characterize this system we have to define a data structure that encodes both the spatial information (the location of cities) and their importance (their population size).

The first step is to construct the dominance tree proposed by Okabe and Sadahiro [START_REF] Okabe | An illusion of spatial hierarchy: spatial hierarchy in a random configuration[END_REF]. The idea is to recursively construct Voronoi tesselations over the set of nodes (see for example the book [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF] and references therein). The Voronoi cell V i of a node i is the set of points that are closer to i than to any other node,

V i = { x | d(x, x i ) < d(x, x j ), ∀j = i}.
We show in Fig. 1(A) an example of a Voronoi tesselation computed for nodes distributed in the plane. Starting from this Voronoi tessellation, we identify local maxima or local 'centers': a point i is a local maximum if its population is larger than those of the neighboring Voronoi cells. In other words, for any point i we define with the help of the Voronoi tesselation the set of neighbors Γ(i) that are the nodes whose Voronoi cell is adjacent to V i . A city i is then a local center if its population P i is larger than the populations of its neighbors: P i > P j ∀ j ∈ Γ(i). In the example of Fig. 1(A), we thus have 3 local maxima labeled a, l, and m. In a second step, we keep these local centers and construct a new Voronoi tessellation over them (Fig. 1(B)). Cities that do not appear anymore (i.e. that were not local centers at the previous step) belong now in the new Voronoi cell of a local center, and are therefore 'dominated' by this city and belong to its 'attraction basin' (for example, nodes b, c, d, f , g, belong to the attraction basin of node a). Here again we determine local maxima, i.e. cells that do not have any neighbor whose population is larger than their own. We repeat this procedure recursively until only one city is left (Fig. 1(C)) which, by construction, is the city with the largest population. We therefore understand that a city will 'survive' many iterations if its population is large but also if it is well located. A large city located very close to an even larger one will quickly be absorbed by this larger city. It is interesting to note that this process bear some ressemblance to the coarse-graining obtained by means of real-space renormalization [START_REF] Kadanoff | Statistical physics: statics, dynamics and renormalization[END_REF]. The dominance tree -shown in Fig. 1(D) for the simple case we just described -is then some sort of bookkeeping of the changes of scale and the hierarchical organization of the marked points: each node has children that correspond to other nodes which belong to its attraction basin. We show on Fig. 1(E-H) the Voronoi tesselations that correspond to the successive steps of this process applied to US cities. In this case the root is New York City (successive Voronoi cells of New York and Los Angeles are colored) and the depth of the tree, that is the number of iterations to reach the root, is H = 6.

We thus represent this process by the dominance tree (Fig. 1(D)) where we keep at each level the remaining points and where the links denote the spatial dominance relation. At height zero (the leaves of the tree), we then have all cities of the system that are locally dominated by a larger neighbor city. After one iteration of the process, we have the local centers, etc. until we reach the root of the tree which is the city with the largest population. Okabe and Sadahiro proposed this construction [START_REF] Okabe | An illusion of spatial hierarchy: spatial hierarchy in a random configuration[END_REF] and used it to prove that the most important quantitative statements in Christaller's central place theory were also observed for a random spatial Poisson point process. More precisely, for both the random Poisson point patterns and their real world case study -marketplaces in Nishinomiya, Japan -they measured Christaller's K indicator, which is the ratio between the number of local centers at consecutive hierarchical levels. They also calculated the statistics of the number of edges of the Voronoi cells as a function of the number of points in the point pattern, and for different families of point patterns.

Consequently their characterization of the structure of the dominance tree was limited, and mainly focused on geometrical aspects.

. We can then represent a system of cities by a dominance tree, where the height of a city is the largest iteration before it is absorbed by a larger city. Once we have constructed the dominance tree, each city i has its height h i . We recall here that, by definition, the height of a node in a tree is the number of edges between this node and the furthest leaf going down in the tree. In other words, it is the length of the longest path (i.e. its number of edges) from the node to the deepest leaf (in contrast, the depth of a node in a tree is the distance from a node to the root). The height of a leaf is then h = 0, while the root has the largest height (denoted by H in the following). We will show that, due to the statistical properties of Voronoi tesselations built from spatial Poisson point processes, H ≈ log 6 (N ) where N is the number of points in the system, and H is in general of order 5 or 6 in our empirical analysis of urban systems. We introduce n f (h) which is the number of cities such that their final height is h. In contrast, the total number of cities at a given height h is denoted by n i (h) and we have the following relation

n i (h) = h ≥h n f (h ).
Once we have constructed the dominance tree, we have to characterize it. There is a large number of possible measures, but we will focus here on the height of a city, and we will discuss it for toy models and for empirical data. A city i at time t is then characterized by its population P i (t) (or equivalently by its rank r i (t) when population are sorted by decreasing order), and its height h i (t) in the dominance tree. This height characterizes the role of the city in the hierarchical organization of the urban system. While in the following we analyze urban systems at the national scale, the method could be applied at different scales as well, such as the regional or continental scale.

Toy models

In order to define a toy model we have to specify both the point distribution and the population distribution. Different models are possible and we will mainly consider the following variants. First, we will consider that populations are distributed according to the power law ρ(P ) ∼ 1/P α with α = 2. This case corresponds to the classical result of Zipf and even if recent data show that this exponent can fluctuate considerably [START_REF] Verbavatz | The growth equation of cities[END_REF][START_REF] Cottineau | MetaZipf. A dynamic metaanalysis of city size distributions[END_REF], it won't affect our discussion here, as the relevant quantity is not the population size itself, but its corresponding rank in the hierarchy: indeed, what matters for the construction of the dominance tree is to know if a city is larger than its neighbors. In the following, we will thus indifferently discuss population or rank for characterizing a given city.

For the point distribution we will consider a uniform distribution (spatial Poisson process) of the x i in the square [-1, +1] 2 . We thus have two lists: a list of population sizes P i and a list of positions x i (i = 1, . . . , N ). In order to define a model we have to specify how to match these lists, and the resulting correlations encode in a simple way the spatial organization of the system. To simplify the description we will assume that the populations are sorted in decreasing order P 1 > P 2 > • • • > P N and that locations are sorted in increasing order according to their distance to the center (0, 0):

||x 1 || < ||x 2 || < • • • < ||x N ||.
We will consider the following three cases:

• The 'deterministic model' : we associate the largest population to the closest point to the center (0, 0): P 1 ↔ x 1 and follow the order P i ↔ x i . This model mimics in some way a system with a central organization around the main city.

• The 'random model': we associate the largest population to the closest point to the center (0, 0):

P 1 ↔ x 1
, and then associate randomly the rest of the P i s to the rest of x i s. In this case, we have a central large city but there is no specific organization around it.

• Finally, the 'tunable model' is less deterministic. We first assign P 1 to x 1 as before, but for the remaining N -1 cities we proceed as follows. For P 2 we choose at random a position among {x 2 , x 3 , . . . , x N } according to the probability p k = exp(-d(x 1 , x k )/L)/Z where d(x 1 , x k ) is the distance between x 1 and x k , Z = k=2,...,n exp(-d(x 1 , x k )/L) is the normalization and L is a positive parameter that can be interpreted as a polarization or interaction distance. We then proceed in a similar way for P 3 with the remaining x i left, etc. until the list of populations We first construct the Voronoi tessellation (shown with dotted lines) for these points and we observe that there are three local maxima: a, l, m, in red and blue. These are the points which do not have any neighbor whose mark is larger than their own. (B) From step (A) we keep the local maxima only, we construct the Voronoi tessellation over this set, and determine the local maxima (which is here the node a). (C) At the end of the process, we are left with the node with the largest mark.

(D) The entire process can be described by the 'dominance tree': leaves correspond to the points which are not local maxima, and the children of a given node are the points that belong to its attraction basin. For example, as can be seen in (B), nodes k, j and e belong to the attraction basin of node l, which is a local maximum at level/height h=1. (E-H) Successive Voronoi tesselations obtained when applying this 'decimation' process for US cities marked with their 2010 population sizes [START_REF] Census | Spatial History Project. United States Historical City Populations, 1790-2010, Center for Spatial and Textual Analysis[END_REF] (we didn't show the case h = 0 which corresponds to the Voronoi tessellation of all US cities and which is unreadable). After H = 6 iterations, we obtain a single city. By construction, the node with the largest mark will always be the root of the dominance tree, but for the other nodes their height in this tree will not only depend on their mark but also on their location.

is exhausted. This tunable model is then able to interpolate between the deterministic case with L 1 and the random case with L 1.

These toy models allow us to explore the effect of various parameters and to test various aspects of our characterization of spatial hierarchies. We show in Fig. 2(A) the Voronoi tessellation of 10 3 points in the square. The color code represents here the rank of cities according to their population (the darker, the larger the associated population). On Fig. 2 we have represented the 'tunable' model that interpolates from the deterministic case (on the left) with a clear 'monocentric' pattern with population decreasing with the distance to the center, to the random model (last figure on the right) where there are no correlations between rank and distance to the center. Before turning to these specific cases, we first discuss some general properties about the depth of the tree, and the number of nodes at a certain level. First, by construction, the root of the dominance tree is the city with the largest population. An important quantity is the depth H of the tree that will measure the number of different hierarchical levels in the system. In order to estimate H, we first evaluate the number of cities n i (h) that are left at the height h. We denote by z h the average number of neighbors in the Voronoi tessellation at level h, and the number of cities at level h is then

n i (h) = N z 1 z 2 . . . z h N e -h ln z (1)
where the brackets • denote the average over all levels. This expression implies that the number of cities and the level h are related as follows

ln n i = a -bh (2)
where a = ln N and b = ln z ≡ ln z eff . The depth of the tree H corresponds to the height of its root, at which there is only city left (n i = 1). Hence we obtain

H = ln N ln z (3)
In the 'random' toy model, cities are uniformly distributed and we can roughly expect z eff ≈ 6 [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF] and with N = 10 4 cities we obtain a tree of height H ≈ 5.

In contrast, for the deterministic model, there is essentially a single local maximum leading to a small depth (see Supplementary Figures S2 andS3)

We consider now the relation between the rank and the height of a city in the dominance tree, and the result is shown for the random model in Fig. 2 fixed locations (in other words, location is the quenched disorder, while the mark is the annealed disorder). This plot shows that on average, smaller ranks lead indeed to larger heights (and eventually the largest city with rank r = 1 has the largest height h = H). However, we also observe on this plot that this relation is not univocal, as there are large fluctuations around the average behavior.

In other words, the rank (or population) of a city does not determine its height in the tree. This shows that the height encodes simultaneously the rank of a city in the urban hierarchy and its location. In particular, we see on Fig. 2(B) that cities with ranks r ≥ 2 can have different heights depending on the disorder realization.

For each value of L for the tunable model, using the general relation Eq. 2, from the measure of n i (h) versus h, we compute the average number z eff = exp(b) of neighbors of the Voronoi tessellations obtained in the decimation process (see Fig. 2(C)). For small L values, this model reproduces the deterministic case where cities are distributed around the largest one in decreasing order.

In this case in many regions there are no local maxima and the effective number of neighbors z eff is very large. As L increases we tend towards the random model where cities are distributed randomly in the plane. In this case we expect a regular Voronoi tessellation constructed over a Poisson point process, with z eff ≈ 6 [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF].

We have shown on these simple toy models that the height of a city in the tree does not depend on its rank alone, but also on its location and its neighborhood. In order to illustrate further the effect of the spatial arrangement of cities on their heights, we consider the following simple situation. We locate the largest city 1 (of rank r = 1) at the center of the plane and we assign to a point at distance d(1, 2) the second largest city 2. We then locate at random all the other cities and measure the height h 2 of city 2 as a function of the distance d(1, 2). We also record the maximum and minimum height reached by 2 in the random configurations. The resulting plot is shown in Fig. 3. We see in this plot that the height h(2) of city 2 depends strongly on the spatial arrangement of cities in the space between it and the larger city (1): if there is a smaller city in between, the city 2 is a local maximum and has a larger height (Fig. 3(A)). In the opposite case, city 2 will be absorbed by city 1 at the next iteration of the decimation process. We test this numerically and the result is shown in Fig. 3(B). First, as expected we observe that the height h( 2) is on average an increasing function of the distance to the largest city. At small distance, the city 2 has a large chance to belong to the attraction basin of city 1 and will be quickly 'absorbed' by it, leading to a small height. In contrast, as the distance increases it will take more steps before the city 2 is in the attraction basin of city 1, leading to a larger height. In addition, for a given distance d(1, 2) we observe important fluctuations of h(2). This clearly shows that the height encodes the spatial organization of the neighborhood of a city, and in particular of the space between the city and the closest larger one. The height is then related to the number of local maxima in this intermediate space.

These toy models demonstrate the relevance of the height of a city in the dominant tree for monitoring its importance in the system. We will now apply these ideas on empirical data.

Empirical studies: US and French urban systems

The previous section helped us to understand the main properties of the dominance tree and the factors that determine the height of a marked point. We will now use this tool to measure the spatial organization of the French and the US urban systems, and their evolution over the last 130 years (see the Data section at the end for details).

Height fluctuations

We first focus on the height of cities as discussed above. We show on Fig. 4(A) and (B) the heights, at different dates, of the 20 cities that are currently the largest ones in the US and in France.

The root of the tree is by construction the largest city of the system, which in both countries has remained the same over the entire period -New York City in the US, and Paris in France. We observe small height variations for French large cities, a fact that was already observed in [START_REF] Batty | Rank clocks[END_REF] for the ranks of old European cities in the urban hierarchy, but with notable exceptions such as Lille. In the US it is interesting to note that even if the population of San Francisco grew steadily since its beginning, its height in the US tree has decreased, signalling a changing environment and the higher growth of neighboring cities.

We observe that in both countries there are different cities which at some point of their history reached the height H -1 (H is the root): in the US Los Angeles, Chicago, or San Francisco; and for France, Marseille, Lyon, Nantes, Bordeaux or Lille. We also observe that the 20 municipalities that are currently the largest ones in France have experienced much less height fluctuations than the 20 largest US cities, some of the latter spanning 3 or 4 different levels in the tree. The height thus captures ascending and descending city trajectories in the US that were not observed in France during the last 130 years. Like other countries in Europe, France has experienced a long-standing urbanization, with a slow and regular evolution, and a city like Paris was already an international hub in the Middle-Age. In contrast, the US is a relatively new country, which experienced successive waves of urban development in the last centuries [START_REF] Bretagnolle | The Organization of Urban Systems[END_REF]. The trajectory of a given city in the dominance tree highlights how the role and the local importance of this city has evolved in time. It provides a richer information that the rank alone [START_REF] Batty | Rank clocks[END_REF], because it encapsulates information not only about the position of the city in the urban hierarchy, but also about its importance among its neighbors.

Basin sizes and statistical definition of cities

Basin sizes: Gini and Zipf

Beyond single cities trajectories, the dominance tree also contains information about the spatial organization of the urban hierarchy at larger scales. At each step of the tree construction process, a local maximum -say city i -contains a number of cities in its Voronoi cell V (i). We denote by φ i the total population of this attraction basin: φ i = j∈V (i) P j (the population of the 'seed' city P i is not taken into account in φ i ). At each step of the iterative construction process we thus have a collection of values of φ, one for each Voronoi cell. We compute the Gini coefficient G of these values (see for example [START_REF] Dixon | [END_REF]): a large value of G (i.e. close to one) indicates that there are a few very large basins of attraction and all the others are small, while in contrast a small G (close to zero) indicates that most attraction basins have roughly the same population size. On Fig. 4(C) we plot G for France and the US at different tree levels h, that correspond to different spatial scales, and for different time periods. There are two remarkable features on this plot. First, in both countries the Gini is decreasing as we approach the root, which means that even if we start from a very heterogeneous situation at small scales (at the urban agglomeration level), at a large scale the population sizes of the attraction basins become comparable. Second, we observe an important difference between the two countries. In France the system has evolved towards a situation where the attraction basins are becoming more heterogeneous (at all spatial scales), while the data show the opposite in the US. Surprisingly enough, it seems that in France, the urban system has not evolved towards a uniform distribution of important basins, but that inequalities in population sizes have increased in the recent periods, at all spatial scales. Furthermore the fact that Zipf's plots of φ for different levels (see Supplementary Figure S1 and the corresponding section in the SI) are all well described by power laws, confirms that the hierarchy of population sizes is preserved at different scales. The power law fits give exponent values displaying a decreas- ing trend, indicating less fluctuations when h increases, in agreement with the results obtained with the Gini coefficient.

Height sensitivity to different statistical definitions of cities

Defining consistent city boundaries is a central issue in urban studies, especially when the goal is to characterize and compare cities [START_REF] Cottineau | MetaZipf. A dynamic metaanalysis of city size distributions[END_REF][START_REF] Pumain | Multilevel comparison of large urban systems[END_REF][START_REF] Cottineau | Diverse cities or the systematic paradox of Urban Scaling Laws[END_REF]. In general, administrative boundaries such as municipalities fail to capture meaningful borders, and many statistical definitions of cities have been proposed. Some are based on morphological aspects (built-up area), some other on functional ones (journey-to-work commuting flows, see Supplementary text in SI for details). While different criteria lead to different boundaries, key variables that quantify the importance of a city in the system should not depend on small fluctuations at finer scales. For example, for cities that are regional, national or international hubs, key variables such as their (population) rank in the urban hierarchy or their height in the dominance tree should not vary too much for different city definitions. We will test this robustness for the french case. In 2021, there were about 2, 400 urban units (morphological definition) in Metropolitan France, that gathered approx. 7, 500 municipalities representing a total population of 50M people. Following the functional definition, there were 682 urban areas, gathering 26, 000 municipalities for a total population of 60M people. For each definition we determine the rank of cities in the population hierarchy, and then attribute their rank to all the municipalities that compose the city according to this definition. Since Paris' urban unit is the most populated in France, all the municipalities that are part of Paris urban unit are attributed rank 1, all the municipalities included in Lyon's urban unit are attributed rank 2, etc. We do the same for the urban area definition: we determine the population hierarchy of urban areas, and attribute their rank to all the municipalities that compose the urban area. The 7k+ municipalities that are simultaneously part of an urban area and of an urban unit thus have two population ranks r i and r i , one according to the morphological definition, the other according to the functional definition. In order to compare these lists of ranks, we compute the Kendall correlation coefficient τ (see for example [START_REF] Abdi | The Kendall rank correlation coefficient[END_REF]). This quantity τ is equal to one when the order of the cities is strictly identical in the two lists, and -1 when the lists have an opposite order. Here, we obtain for these lists a value τ r = 0.28. We then compute the dominance tree for France according to the two different city definitions, and determine the height of each city in both trees (we allocate this height to all the municipalities that compose the city according to this definition). For the 7k+ municipalities that are part of a urban unit and also part of a urban area, we end up with two lists of heights, and we compute the Kendall τ coefficient between these two lists. We obtain τ h = 0.30, and we thus have τ h > τ r , which shows that the height of a city in the dominance tree is in fact slightly less sensitive to different city definitions than its rank in the population hierarchy.

Relation between city size and height in the tree

We also investigate the relation between the population size of a city and its height in the dominance tree. In order to quantity this relation, we compute two lists in which cities are ranked in decreasing order of (i) population sizes and (ii) heights. We compare these lists using Kendall's τ correlation coefficient. For both countries we calculate this quantity for each year available in the data (see Supplementary Figure S4), and we also calculate it for the toy models previously discussed. In particular, we obtain for the 'random' toy model the average value τ ≈ 0.43 (represented on Supplementary Figure S4 by the dashed line). For the determinisc, monocentric toy model we obtain a value close to 0, and for the tunable model, τ interpolates between these two values for different L. Overall for the US and France, we observe relatively small values of Kendall's τ in these systems, confirming that the height is not determined by population alone, and that the location of a city has its own importance and effect. Supplementary Figures S5 and S6 complement this observation, and show that there are many pairs of cities such as P i > P j and h i < h j .

We can understand these different results with the following considerations. When large cities are close to one another, some are dominated by even larger ones, and consequently their height in the tree is small. This leads to a small value of τ as observed for the deterministic model. In contrast, when the population size and the location are uncorrelated as in the random model, large cities can reach higher levels in the tree before being 'absorbed' by an even larger city. In this case, there is a higher correlation between rank and height, leading to a larger value of τ . We observe that τ is smaller for France than for the US, suggesting that the French system of cities displays a stronger concentration of population within large urban agglomerations. These are composed of large municipalities located close to one another, in the vicinity of a historical large municipality that is the core of the agglomeration. In addition, we observe that τ is decreasing in time for both countries, signalling a decrease of the correlation between the height and the rank in the urban hierarchy. This can be due to two different effects: smaller cities becoming more important (large h) due to their strategical location, or in contrast cities with small h whose population size increased because they are located near an important city from which they depend.

Height clocks

A city i of population P i (t) has thus a certain height h i (t) in the dominance tree which characterizes the importance of its role at a regional level. In contrast, the leaves of the tree (with height h = 0) do not dominate any other city and depend directly on a more important city in their vicinity. In order to characterize the evolution of the urban system and how different cities can see their role evolving in this system, we study the temporal variation of their height. A first approach was shown in Fig. 4(A,B), but for a better visualization we propose a 'height clock', in the same spirit as the rank clock introduced by Batty in [START_REF] Batty | Rank clocks[END_REF]. On a polar plot, the radius is equal to the normalized height max(h) -h i (t) and the angle is proportional to time. During a time range [t 1 , t 2 ], the total height variation for a city i is then given by ∆ i = We illustrate here different classes of possible dynamics with selected cases: stable with fluctuations (top line), increasing height (middle) and decreasing height (bottom) in the dominance tree. On each plot, we indicate the percentage of cities in the urban system that belong to that class. The central node corresponds to the root of the system, i.e. the city with the largest height. The color code corresponds to the one used in Fig. 4(A) and (B).

US

The first type of cities display a constant height ('stable' cities, not shown here). For France this is for example the case of Paris (and New York City for the US), which was always the most important city in the system since 1880. Toulouse in France has a constant height equal to 4 (H -2, which can also be seen in Fig. 4), which indicates the steady importance of this city at a regional level. This class of cities is by far the most important for both countries with a share of 86.7% for France and 72.4% for the US. This is an important indication that the French system reached earlier some kind of more stable state, compared to the US where more active change dynamics were observed in the last century.

The second type comprises cities for which h(t 1 ) = h(t 2 ), but which experienced height fluctuations in between. This is the case for example of Mulhouse, Nantes or Lyon in France, or Omaha in the US. As we can see, Mulhouse saw its height fluctuating between 1870 and 1970 before stabilizing in the last decades. During the whole period its population grew steadily (except during the 2nd world war) and this again shows that the height analysis encodes more than just the population dynamics, by considering also the population dynamics of neighboring cities. The third type correspond to cities whose height decreased, while the fourth type are cities whose height increased and which became more important in the organization of the system. The proportion of cities that belong to one of these two categories is larger in the US than in France, which is consistent with the urban dynamics observed in both countries during the 20th century [START_REF] Bretagnolle | The Organization of Urban Systems[END_REF]. Obviously these results need to be replaced in a geographical context, but these observations prove that non-trivial information can be extracted through this single quantity h, that can shed light on the evolution of cities and their importance in the urban system.

To summarize, the evolution in time of the height of a point in the dominance tree sheds a light on the evolution of its importance in the spatial hierarchy. This information is complementary to the time evolution of its rank according to the value of its mark, as illustrated by Supplementary Figure S8 that shows the height clock of three large US cities alongside their rank clock [START_REF] Batty | Rank clocks[END_REF]. The height and the rank trajectories of a given city can even go in opposite directions, as this is the case for example for Portland and Dallas.

Discussion

An important problem in spatial statistics is to characterize a system of points that have marks. For cities, the simplest mark is the population size and the importance of a city at different geographical scales results from both its location and its rank in the urban hierarchy, determined by its population size. We presented a simple tool based on the dominance tree that encodes the spatial structure of a marked point process.

Starting from an original method proposed by [START_REF] Okabe | An illusion of spatial hierarchy: spatial hierarchy in a random configuration[END_REF], we have introduced new metrics, toy models, and new vizualisations for this problem. In particular, we have conducted a proper analysis of the influence of the spatial distribution of marks in space on the maximal height a point with a given rank can reach in the tree. We have introduced φ, the population size of the attraction basins, and we have also introduced the height clock that allows to monitor the time evolution of the height of a particular city. Finally we have also introduced the height vs. rank, and height vs. population size diagrams (see Supplementary Figures S5 andS6).

The height of a point in this tree appears to be a crucial information, and its monitoring allows to understand the dynamics of spatial hierarchies. We have applied this method to toy models and to empirical analysis of the evolution of the French and the US systems of cities. This method however goes well beyond these examples and could be applied to a wealth of problems, depending on the nature of points and on the quantity chosen to characterize them. It allows for the study of spatial correlations in a simple way, and could help in comparing the dynamics of many different urban systems, discussing their evolution in time, classifying them, etc. and finally contributing quantitatively to urban theory thanks to a shared, non-parametric tool. It could also be useful in connecting distinct bodies of literature, namely the studies of the statistical distribution of settlement and city sizes [START_REF] Verbavatz | The growth equation of cities[END_REF][START_REF] Cottineau | MetaZipf. A dynamic metaanalysis of city size distributions[END_REF], and those that investigate their spatial distribution on the Earth [START_REF] González-Val | The Spatial Distribution of US Cities[END_REF][START_REF] Strano | The agglomeration and dispersion dichotomy of human settlements on Earth[END_REF]. Further studies could also investigate the relation between the height of cities and socio-economical factors, such as their GDP.

Concerning the application of this method to urban systems, several lines of future research can be identified. The first one would be to adapt the method so that it could be applied not only to points but also to cellular tissues, by considering the real geographical envelopes of the geographical units rather than the Voronoi cells computed from their centroids. Since the set of municipalities, departments, regions, etc. in a country naturally form tessellations, the rest of the procedure would remain the same: identify the cells that are the local maxima, aggregate the dominated cells to the one that is a local maxima, and so on. At a smaller, intra-urban scale, it could also be applied to capture the dynamical organisation of hotspots within the city.

Future studies may also further investigate the sensitivity of the method against different definitions of cities, that would correspond to different marked point processes or to different elementary building blocks used for building the dominance tree. In the same vein, it would be interesting to measure differences in the tree structure when considering more elaborated distances than the Euclidean distance. One could for example consider transport networks distances, or average travel times. Such data are however much harder to obtain for many different countries, and for ancient periods of time. Another direction could be to analyze the sensitivity of scaling laws exponent for different hierarchical levels in the dominance tree [START_REF] Xiao | Removing spatial autocorrelation in urban scaling analysis[END_REF].

Beyond urban geography, other problems in spatial statistics could also be addressed with this simple tool. As it encodes in a simple way various correlations, it could shed a new light on the characterization of random tessellations and random patterns in two dimen-sions, and could help in revisiting important problems in the physics of foams, such as the Aboav-Weaire and Lewis laws [START_REF] Chiu | Aboav-Weaire's and Lewis' laws-A review[END_REF][START_REF] Weaire | Soap, cells and statistics -random patterns in two dimensions[END_REF] which in general focus on the distribution of local statistics such as the number of edges, the neighbor statistics, etc. but not on the global spatial organization.

Data

We use census population data for France provided by the French national statistics bureau (INSEE) [START_REF]Recensements de la population 1876-2018[END_REF]. These data give the population sizes of approximately 32, 000 French municipalities (administrative boundaries) of Metropolitan France for every decade during the period 1876-2018. For the US, the dataset [START_REF] Census | Spatial History Project. United States Historical City Populations, 1790-2010, Center for Spatial and Textual Analysis[END_REF] is a compilation of US cities populations between 1790 and 2010 (every 10 years). The data come primarily from the US Census Bureau dataset and contain 7, 500 incorporated cities which at a certain point of their history had a population larger than 2, 500 individuals. Other cities were added from a variety of sources (see [START_REF] Census | Spatial History Project. United States Historical City Populations, 1790-2010, Center for Spatial and Textual Analysis[END_REF] for details). Both datasets are public and can be freely downloaded (see [START_REF] Census | Spatial History Project. United States Historical City Populations, 1790-2010, Center for Spatial and Textual Analysis[END_REF][START_REF]Recensements de la population 1876-2018[END_REF] for the URLs).

We note that the underlying criteria used to delin-eate cities in space for computing their population size may influence our results. Here we use administrative boundaries (municipalities) for France which are known to be problematic in many ways when comparing historical population data from one city to another, or between countries [START_REF] Pumain | Multilevel comparison of large urban systems[END_REF]. However, the dominance tree is constructed in order to highlight the spatial organization of urban systems at various scales and while most methods used to determine relevant spatial boundaries for cities rely on arbitrary thresholds that are different from one country to another (shares of commuting flows, minimal distance between consecutive built-up areas, etc.), this non-parametric method could provide an alternative way to construct urban agglomerations with Voronoi cells at different heights, and could constitute an interesting direction for future studies.

FIG. 1 :

 1 FIG.1: Construction of the dominance tree. (A) Each point is characterized by a mark/value represented here as the size of the circle. We first construct the Voronoi tessellation (shown with dotted lines) for these points and we observe that there are three local maxima: a, l, m, in red and blue. These are the points which do not have any neighbor whose mark is larger than their own. (B) From step (A) we keep the local maxima only, we construct the Voronoi tessellation over this set, and determine the local maxima (which is here the node a). (C) At the end of the process, we are left with the node with the largest mark. (D) The entire process can be described by the 'dominance tree': leaves correspond to the points which are not local maxima, and the children of a given node are the points that belong to its attraction basin. For example, as can be seen in (B), nodes k, j and e belong to the attraction basin of node l, which is a local maximum at level/height h=1. (E-H) Successive Voronoi tesselations obtained when applying this 'decimation' process for US cities marked with their 2010 population sizes[START_REF] Census | Spatial History Project. United States Historical City Populations, 1790-2010, Center for Spatial and Textual Analysis[END_REF] (we didn't show the case h = 0 which corresponds to the Voronoi tessellation of all US cities and which is unreadable). After H = 6 iterations, we obtain a single city. By construction, the node with the largest mark will always be the root of the dominance tree, but for the other nodes their height in this tree will not only depend on their mark but also on their location.

FIG. 2 :

 2 FIG. 2: (A) Initial Voronoi tessellations computed for N = 10 3 points. From left to right, the tunable toy model interpolates from the deterministic model (the population decreases with the distance to the central largest city) to the random model where distance and population to the center are uncorrelated. (B) Rank versus height for the random toy model, computed for N = 10 4 cities and 200 disorder configurations. The straight line is a linear fit of the form log(r) = 8.555 -1.722 * h (R 2 > 0.99). The points show the outliers and the lower and upper sides of the boxes correspond to the second quartile and third quartile respectively (the boxes therefore contain half of the points). The middle bar corresponds to the median. (C) Effective connectivity of the Voronoi tessellation for different values of L in the tunable model.

FIG. 3 :

 3 FIG.3: (A) Two different possible configurations: (left) smaller cities are located between cities 1 and 2 or (right) no cities are in between. In the former case, the city 2 is a local maximum, while in the latter, city 2 belongs to the dominance basin of city 1. (B) Average value of the height h2 of city 2 as a function of its distance d(1, 2) to the largest city (P1) in the random model, normalized by the largest distance possible. Error bars represent here max(h2) -min(h2) observed over 200 random configurations for the cities.

FIG. 4 :

 4 FIG.4: Height in the dominance tree for different cities in the US (A) and in France (B). The amount of transparency corresponds to the time period (following the chronological order from faded to clear) and the star denotes the city's height for the most recent year in our dataset (2010 for the US and 2018 for France). The color of the dots correspond to the class of cities they belong to, defined in Fig.5. (C) Evolution of the Gini coefficient for the population sizes φ of the attraction basins of cities with a given height h. The results are shown for France and the US, and the opacity follows the chronological order (almost transparent points correspond to the earliest time periods).

FIG. 5 :

 5 FIG.5: Height clocks for French (left column) and US municipalities (right column). We illustrate here different classes of possible dynamics with selected cases: stable with fluctuations (top line), increasing height (middle) and decreasing height (bottom) in the dominance tree. On each plot, we indicate the percentage of cities in the urban system that belong to that class. The central node corresponds to the root of the system, i.e. the city with the largest height. The color code corresponds to the one used in Fig.4(A) and (B).

Acknowledgements

The authors thank the anonymous reviewers for their useful comments that helped improving the manuscript. TL thanks the many developers of the free software packages that were used in this study, especially sf [33], GDAL, GEOS and PROJ to perform the geocomputation and ggplot2 [34] to draw the figures. Last but not least, TL warmly thanks the Ramelet family for hosting him in the winter of 2021 during the COVID-19 lockdown, during which most of this work was done.

Data and source code availability. The data that support the findings of this study are available for download at https://github.com/cestastanford/ historical-us-city-populations for the US data, and at https://www.insee.fr/fr/information/ 2414405 for the French data. The source code is available at https://gitlab.huma-num.fr/tlouail/voronoize.

Author Contributions

TL and MB designed the study. TL wrote the code and performed calculations. TL and MB analyzed and interpreted the results, prepared the figures and wrote the manuscript.

Competing Interests

The authors declare that they have no competing interests.


A Morvan

L Calvez

A Vella

A Saiter-Fourcin

"Physical aging of the 62.5GeS2-12.5Sb2S3-25CsCl chalcogenide glass: assessing the mechanisms of equilibration and crystallization"

Keywords: Physical Aging, Crystallization, Chalcogenide glass

The goal of this work is to provide an answer to the important scientific question: will a glass-forming liquid reach its equilibrium state or not before to crystallize during physical aging process, i.e. when it is maintained below its glass transition temperature? For this goal the physical aging process in a 62.5GeS2-12.5Sb2S3-25CsCl glass has been investigated for a long period (5 months) at 30°C below the glass transition temperature Tg. A crystallization was observed during this process, and was characterized through X-Ray Diffraction (XRD), Transmission Electron Microscopy (TEM) and Temperature Modulated Differential Scanning Calorimetry (MT-DSC) analyses. By evaluating the enthalpy loss related to the physical aging and by comparing this experimental quantity to the theoretical one, we have proved that this glass has reached its thermodynamic equilibrium state before crystallizing. Thus, this work evidences experimentally the ability for a glass former to crystallize under Tg after reaching its equilibrium state.

Introduction

Chalcogenide glasses are a centre of interest since many years now, due to their ability to transmit light from visible to far infrared range and their capacity to be easily shaped in different geometries, these advantages offering great potential in the development of new infrared cameras and optical fibres for the mid-infrared regions [START_REF] Rogalski | Infrared Detectors: An Overview[END_REF][START_REF] Adam | Chalcogenide Glasses: Preparation, Properties and Applications[END_REF] . They also have an economic advantage because these glasses are cheaper than currently used single crystal germanium or ZnSe polycrystals.

Moreover, the non-equilibrium nature of these glasses is related to the structural relaxation occurring under the glass transition temperature [START_REF] Hutchinson | Physical Aging of Polymers[END_REF][START_REF] Ma | Physical Aging of Ester-Cross-Linked Hollow Fiber Membranes for Natural Gas Separations and Mitigation Thereof[END_REF][START_REF] Salmerón Sánchez | Influence of the chemical structure on the kinetics of the structural relaxation process of acrylate and methacrylate polymer networks[END_REF][START_REF] Hodge | Enthalpy Relaxation and Recovery in Amorphous Materials[END_REF] . Many studies have been recently reported to understand the glassy state nature and the associated physical aging, with a focus on polymeric materials [START_REF] Schawe | Competition between Structural Relaxation and Crystallization in the Glass Transition Range of Random Copolymers[END_REF][START_REF] Buquet | Physical Ageing and Molecular Mobilities of Sulfonated Polysulfone for Proton Exchange Membranes[END_REF][START_REF] Cangialosi | Effect of Nanostructure on the Thermal Glass Transition and Physical Aging in Polymer Materials[END_REF] and few on chalcogenide glasses [START_REF] Golovchak | Step-Wise Kinetics of Natural Physical Ageing in Arsenic Selenide Glasses[END_REF][START_REF] Shpotyuk | 8 -Physical Ageing of Chalcogenide Glasses[END_REF] . The recent studies on polymer materials highlighted the relaxation and crystal nucleation/growth of several samples once the equilibrium is reached [START_REF] Schawe | Competition between Structural Relaxation and Crystallization in the Glass Transition Range of Random Copolymers[END_REF][START_REF] Androsch | Relaxation and Crystal Nucleation in Polymer Glasses[END_REF][START_REF] Fokin | Effect of Structural Relaxation on Crystal Nucleation in Glasses[END_REF][START_REF] Androsch | Sequence of Enthalpy Relaxation, Homogeneous Crystal Nucleation and Crystal Growth in Glassy Polyamide 6[END_REF][START_REF] Grassia | Complete Set of Enthalpy Recovery Data Using Flash DSC: Experiment and Modeling[END_REF] .

The equilibrium state is reached when the relaxation enthalpy of the glass is equal to the maximum relaxation enthalpy as described in Figure 1 [START_REF] Cowie | The Ageing of Poly(Vinyl Methyl Ether) as Determined from Enthalpy Relaxation Measurements[END_REF] . This maximum of enthalpy loss is noted ΔH∞ and can be estimated by the equation 1:

ΔH ∞ = ΔC 𝑝 (T 𝑔 -T 𝑎 ) [equation 1]
With ΔCp the difference between the heat capacities of the liquid and the glassy states at the glass transition temperature Tg, and Ta the aging temperature. ΔH∞ represents the maximum enthalpy loss when the glass has reached its equilibrium state corresponding to the extrapolated curve of the liquid like state.

ΔHcryst represents the enthalpy loss if the glass crystallizes after reaching the equilibrium state (we may note that the crystal state is not reached in usual studies).

The two images are proposed to illustrate the transparency loss after the physical aging process for the glass studied in the present work. This loss is directly visible with eyes.

Our work focused on chalcogenide glass of the composition 62.5GeS2-12.5Sb2S3-25CsCl, which has proved its capability to produce glass-ceramics with controlled crystallization [START_REF] Ma; Calvez | Crystallization Study of Infrared Transmitting Glass Ceramics Based On[END_REF][START_REF] Xia | Glass Formation and Crystallization Behavior of a Novel GeS2-Sb2S3-PbS Chalcogenide Glass System[END_REF][START_REF] Zhang | A New Class of Infrared Transmitting Glass-Ceramics Based on Controlled Nucleation and Growth of Alkali Halide in a Sulphide Based Glass Matrix[END_REF][START_REF] Zhu | Optical and Mechanical Properties of Far Infrared Transmitting Glass-Ceramics[END_REF] . In a previous work [20], the authors had evidenced a glass transition temperature at the onset point of 280°C, and a crystallization peak with a maximum around 478°C during heating scans. Since physical aging can influence the material performances [START_REF] Fokin | Effect of Structural Relaxation on Crystal Nucleation in Glasses[END_REF][START_REF] Naeem Iqbal | Enthalpy Relaxation of the Glass of Poly (l-Lactic Acid) of Different d-Isomer Content and Its Effect on Mechanical Properties[END_REF][START_REF] Hill | The Effects of Physical Aging in Polycarbonate[END_REF] it is important to study and give a better understanding of this phenomenon. Moreover, the way by which a glass reaches (or not) its equilibrium state after a certain aging time is still matter of debate in the scientific community [START_REF] Cowie | The Ageing of Poly(Vinyl Methyl Ether) as Determined from Enthalpy Relaxation Measurements[END_REF] .

Recently, it has been evidenced that this way is clearly dependent of the gap between Tg and Ta in different glass formers: polymeric systems [START_REF] Cangialosi | Direct Evidence of Two Equilibration Mechanisms in Glassy Polymers[END_REF] , metallic glasses [START_REF] Gallino | Hierarchical Aging Pathways and Reversible Fragile-to-Strong Transition upon Annealing of a Metallic Glass Former[END_REF] , and chalcogenide glasses [START_REF] Golovchak | Dynamics of Structural Relaxation in Bioactive 45S5 Glass[END_REF] . The use of the fast-scanning calorimetry (FSC) has generated many recent studies [START_REF] Schawe | Competition between Structural Relaxation and Crystallization in the Glass Transition Range of Random Copolymers[END_REF][START_REF] Monnier | Distinct Dynamics of Structural Relaxation in the Amorphous Phase of Poly( L -Lactic Acid) Revealed by Quiescent Crystallization[END_REF][START_REF] Androsch | Relaxation and Crystal Nucleation in Polymer Glasses[END_REF] proving that the physical aging can accelerate the ability to crystallize by forming nuclei after reaching the equilibrium state. This phenomenon known from a theoretical point of view is difficult to evidence experimentally because it appears after the complete structural relaxation, but we highlighted it through the study of a selenium glass considered as an ideal glass former according to its simple structure [START_REF] Morvan | Physical Aging of Selenium Glass: Assessing the Double Mechanism of Equilibration and the Crystallization Process[END_REF] . This previous study highlighted that for aging temperature closed to Tg, the molecules constituting the glass have enough mobility to enhance the crystallization process after reaching the thermodynamic equilibrium. But in this last study, the highlight of the crystalline phase presence was not possible through XRD and TEM experiments because of the small quantity of crystals obtained by this way, and the own nature of the selenium, which sublimes under electron beam.

The challenge of this work is to investigate for the first time the physical aging process of a chalcogenide glass with a complex composition (62.5GeS2-12.5Sb2S3-25CsCl)

presenting the ability of changing the optical transmission when a very small quantity of crystals is formed [START_REF] Morvan | Nanoscale Crystallization Mechanisms in a GeSSbCsCl Glass Ceramic and Relationships with Mechanical and Optical Properties[END_REF] , and to provide an answer to the important scientific question:

will a glass-forming liquid reach its equilibrium state (extrapolated curve of the liquid like state) or not before to crystallize during a physical aging process?

Experimental

The chalcogenide glass 62.5GeS2-12.5Sb2S3-25CsCl was synthetized using the meltquenching technique, using high purity materials (5N Ge, Sb, S, Ga and 3N CsCl). The elements were weighted in a glove box and then introduced in a silica tube.

The preparation was evacuated down to 10 -5 Pa, sealed and heated in a rocking furnace at 850°C for 12h. The melt was homogenised and quenched in water at room temperature and placed in a furnace 10°C below Tg in order to avoid residual mechanical stress [START_REF] Adam | Chalcogenide Glasses: Preparation, Properties and Applications[END_REF][START_REF] Zhang | Production of Complex Chalcogenide Glass Optics by Molding for Thermal Imaging[END_REF][START_REF] Wang | Thermal and Optical Properties of GeS2-Based Chalcogenide Glasses[END_REF] .

XRD measurements were performed on a Bruker diffractometer with a Co target (K = 1.79 Å) equipped with a furnace. A step width of 0.04° was used for all the analyses with a scan speed of 1s in a θ-2θ mode, and the measurements were performed on powdered samples.

Unaged and aged materials were characterized by Transmission Electron Microscopy

(TEM) on a JEOL ARM 200F microscope operated at 200kV. Bright field images and selected area electron diffraction (SAED) patterns were recorded with a parallel beam.

Samples were prepared with a three-step procedure: (i) mechanical crushing of the glass in a mortar, (ii) dispersion of crushed pieces in ethanol and sonication, and (iii) deposition of a droplet on a TEM grid with a carbon film. TEM observations were carried out in the thinnest areas of fragments that were transparent to electrons.

Additional data were collected in Scanning TEM (STEM) mode with a probe size of 0.2 nm and a convergence angle of 34 mrad. Images were recorded with both a bright field (BF) detector (collection angles 0-45 mrad) and a High-Angle Annular Dark-Field (HAADF) detector (collection angles 68-280 mrad).

Modulated temperature differential scanning calorimetry (MT-DSC) experiments were performed with the Q100 apparatus from TA instrument. Temperature and energy were calibrated using indium as standard, and heat capacity using sapphire. All the experiments were performed with samples of about 5 mg inside aluminium pans under a 50mL/min nitrogen gas flow. The reference and sample pans had identical masses within +/-0.05mg to optimize the baseline. The temperature modulation was performed with a heating rate of 5K/min, a modulation amplitude of ±0.796K and an oscillation period of 60s. This protocol called "heat-only mode" is the most accurate to analyse samples with an ability to crystalize during a heating ramp [START_REF] Saiter | Probing the Chain Segment Mobility at the Interface of Semi-Crystalline Polylactide/Clay Nanocomposites[END_REF] .

The physical aging was performed during five months in a furnace at 250°C (under air), i.e. 30°C below the glass transition temperature [START_REF] Morvan | Nanoscale Crystallization Mechanisms in a GeSSbCsCl Glass Ceramic and Relationships with Mechanical and Optical Properties[END_REF] , on disk shape samples (diameter of 10 mm and thickness of 2 mm). These disk shape samples were directly cut from the ingot obtained inside the ampoule. After the aging time of 5 months, the samples were cooled down directly at room temperature. During a very long time in a ventilated furnace a very thin film (below 1µm) of oxidation can be observed (irisation) that is then removed by a very slight polishing with alumina powder.

Results and discussion a-Unaged sample

Before the aging process, the microstructure of the sample (unaged sample) was checked through XRD analysis and TEM with SAED patterns. As shown in Figure 2 the X-ray diffraction pattern corresponds to a fully amorphous sample. Even if the X-ray diffraction highlights an amorphous sample, some nuclei can be already present without being detected by XRD analyses. TEM analyses were performed in order to confirm the amorphous state at the nano scale for the unaged sample. The results are presented in Figure 3.

Figure. 3: TEM bright field image on the unaged sample and corresponding SAED pattern (inset).

The TEM image confirms the amorphous state of the unaged sample, with the absence of spot on the SAED pattern.

b-Aged sample

The sample was aged during five months in a furnace at 250°C, i.e. 30°C under the glass transition temperature. Once the aging is over, the transparency modification is directly evidenced and due to a partial crystallization of the glass [START_REF] Zhu | Optical and Mechanical Properties of Far Infrared Transmitting Glass-Ceramics[END_REF][START_REF] Morvan | Nanoscale Crystallization Mechanisms in a GeSSbCsCl Glass Ceramic and Relationships with Mechanical and Optical Properties[END_REF] as shown in Figure 1 (see images).

The aged sample was also imaged by TEM experiments in order to observe the crystals formed during the physical aging. Figure 4 shows an HAADF-STEM image of the aged sample. The TEM image in dark field contrast should exhibit bright spots corresponding to CsCl crystals, but, as explained in details in our previous work [29], the solvent used for the sample preparation dissolves partially or completely these crystals, and only holes (left by the dissolved crystals) or dark spots are observed.

A distribution of crystal size between 10 and 50 nm is observed, with an average value around 30 nm, suggesting that the crystal growth started just after the nucleation. The hypothesis is that the nucleation rate is low (small driving force 30°C below the glass transition temperature) or the nucleation barrier is high, with a fast growth (high atomic mobility).

In order to estimate the crystallinity degree (Xc) of the aged glass, two methods were used. The first one uses the XRD data, by calculating the ratio between the peak area associated to the crystalline part and the full area (which corresponds to amorphous and crystalline parts), as illustrated in Figure 5.

Figure 5: X-ray diffraction pattern of a 5-month aged glass at 250°C with the different fits (highlighted with different colours) used to estimate the area of crystal (peaks) and amorphous parts (halo).

The estimation for the crystallinity degree by XRD analysis is Xc = 10 ± 2 %. But this method is more qualitative than quantitative.

The second method is based on calorimetric analyses, and consists on comparing the heat capacity step value (noted ΔCp) at the glass transition temperature taken at the middle point (noted Tg mid) : ΔCp(Tg mid) = Cp,liq(Tg mid) -Cp,glass(Tg mid) , between the unaged sample (the amorphous one) and the aged one. As already highlighted in our previous work [START_REF] Morvan | Nanoscale Crystallization Mechanisms in a GeSSbCsCl Glass Ceramic and Relationships with Mechanical and Optical Properties[END_REF] a cold crystallization phenomenon occurs in the glass transition range, and MT-DSC analyses were needed in order to separate the contribution related to the heat capacity (shown in the "reversing signal" as the glass transition), to the contribution related to kinetic events (shown in the "non-reversing signal" as cold crystallization, structural relaxation, etc). The reversing signals corresponding to the unaged and aged samples were superimposed and shown in Figure 6. As shown in Figure 6, the impact of the crystallization is clearly observed for the aged sample (green curve) with the ΔCp decrease from 0.36 ± 0.01 J/(g.K) down to 0.32 ± 0.01 J/(g.K). It is well known that the heat capacity step at the glass transition temperature is proportional to the amorphous phase quantity. Thus, whatever the crystallization process, if Xc increases, Cp will decrease. Furthermore, the ratio between Cp(crystallized system) and Cp(amorphous system) is directly proportional to Xc :

𝑋 𝐶 (%) = ∆C p(crystallized) ∆C p(amorphous)
. 100

[equation 2]
According to the values obtained in this work, the crystallinity degree is around 11 ± 2 %, which is in agreement with the crystallinity degree estimated from XRD analysis (Xc,XRD = 10 ± 2 %). Furthermore, the Tg increase from 280°C up to 290°C is also a marker of the presence of the crystalline phase. Indeed, the aging process usually implies a Tg decrease, but if a crystallization process occurs, the relaxation processes are constrained by the presence of crystals, and needs more energy to devitrify, implying an increase of the glass transition temperature.

As explained in the introduction, a glass maintained at a certain temperature below its glass transition temperature, will normally reach its equilibrium state corresponding to the extrapolated curve of the liquid-like state before starting a nucleation process. To confirm experimentally this theoretical scenario, the enthalpy loss of the aged glass should have reached (as a minimum) the ΔH∞ value (see equation 1). In our case this theoretical value is equal to 10 J/g by taking the values of Tg = 280°C, ΔCp = 0.36 J/(g.K), and Ta = 250°C.

As the physical aging process concerns only the amorphous phase, the calorimetric data were normalized to the amorphous fraction, i.e. 90% for the aged sample, in order to estimate the enthalpy loss in J/g. The Figure 7 presents the non-reversing heat flow signal associated to the aging process (and the cold crystallization starting around 275°C as shown in our previous work [START_REF] Morvan | Nanoscale Crystallization Mechanisms in a GeSSbCsCl Glass Ceramic and Relationships with Mechanical and Optical Properties[END_REF] ) of the aged sample. The same curve but corresponding to the unaged sample was superimposed to highlight the difference in terms of calorimetric signature. The enthalpy loss related to the aging process is measured from the green curve presented in Figure 7 (red area), and is equal to H = 11 ± 0.2 J/g. H value being very close to ΔH∞ = 10 J/g, this MT-DSC experiment confirms that the amorphous phase of the aged glass has reached its equilibrium state. Thus, this experimental work is in good agreement with the theory exposed in the introduction part, i.e. a glass can crystallize under its glass transition temperature once its thermodynamic equilibrium state is reached.

Cold crystallization

Conclusions

The goal of this work was to provide an answer to the important scientific question: will a glass-forming liquid reach its equilibrium state or not before to crystallize during a physical aging process? For this goal, this physical phenomenon has been investigated for the first time in the 62.5GeS2-12.5Sb2S3-25CsCl glass for a very long period (5 months) at 30°C below the glass transition temperature Tg. XRD and TEM analyses on the unaged sample confirm the fully amorphous state of this sample. The sample aged 30°C below the glass transition temperature during 5 months has lost its transparency as expected for this specific glass when a certain crystallinity degree is reached. This lost is visible directly with eyes but in order to quantify the crystallinity degree (Xc), data obtained by XRD and MT-DSC analyses were analysed and Xc is estimated around 10%. In order to confirm or not if the chalcogenide glass composition 62.5GeS2-12.5Sb2S3-25CsCl exhibits the expected theoretical behaviour, i.e. a nucleation/growth can be observed only once the aged glass has reached its equilibrium states, the enthalpy loss related to the aging of 5 months is measured by MT-DSC experiments. The value obtained experimentally was found very close to the theoretical one, which corresponds to the enthalpy loss if the glass has reached its thermodynamic equilibrium state. Thus, this new and original experimental work is in good agreement with the theory proposed by the scientific community defending the possibility for a fully amorphous glass to crystallize below its glass transition temperature if it has reached its equilibrium state.

Figure 1 :

 1 Figure 1: Enthalpy of an initially liquid state material as a function of temperature during physical aging at a temperature (Ta) below the glass transition temperature Tg.

Figure 2 :

 2 Figure 2: X-ray diffraction pattern of the unaged sample.

Figure 4 :

 4 Figure 4: HAADF-STEM image of the sample aged during 5 months at 250 °C. Dark discs correspond to partially or totally dissolved CsCl crystals.

Figure 6 :

 6 Figure 6: MT-DSC reversing heat flow signals of the unaged (amorphous) glass (blue curve) and the glass aged for 5 months at 250°C (green curve). The curves have been normalized to the sample mass.

Figure 7 :

 7 Figure 7: MT-DSC non-reversing heat flow signals of the aged and unaged glasses. The area in red corresponds to the enthalpy loss related to the aging process. The curve corresponding to the aged sample was normalized to the amorphous fraction, i.e. 90%.

Acknowledgements

Xavier Sauvage from the GPM laboratory (UMR CNRS 6634, University of Rouen, France) is warmly acknowledged for helping during TEM experiments.