Pierre Dupont
email: pierre.dup@live.fr

Dominique Allainé

Mariona Ferrandiz-Rovira

Roger Pradel

Efficient spatial multi-state capture-recapture model to study natal dispersal. An application to the Alpine marmot

Keywords: Bayesian modelling, Dispersal kernel, Multi-state capture-recapture

1. Studying natal dispersal in natural populations using capture-recapture data is challenging as an unknown proportion of individuals leaves the study area when dispersing and are never recaptured. Most dispersal (and survival) estimates from capture-recapture studies are thus biased and only reflect what happens within the study area, not the population.

2. Here, we elaborate on recent methodological advances to build a spatially-explicit multistate capture-recapture model to study natal dispersal in a territorial mammal while accounting for imperfect detection and movement in and out of the study area.

3. We validate our model using a simulation study where we compare it to a non-spatial multistate capture-recapture model. We then apply it to a long-term individual-based dataset on Alpine marmot (Marmota marmota). 4. Our model was able to accurately estimate natal dispersal and survival probabilities, as well as mean dispersal distance for a large range of dispersal patterns. By contrast, the nonspatial multi-state estimates underestimated both survival and natal dispersal even for short dispersal distances relative to the study area size.

5. We estimated higher inheritance probabilities of female Alpine marmots, which suggests higher levels of philopatry, although the probability to become dominant after dispersal did not differ between sexes. Nonetheless, the lower survival of young adult males suggests higher costs of dispersal for males. We further discuss the implications of our findings in light of the life-history of the species.

6. The flexibility of our spatially-explicit multi-state capture-recapture model makes it adjustable enough to be widely applicable in other species.

INTRODUCTION

Dispersal, and especially natal dispersal, is a fundamental process in biology [START_REF] Dobson | The enduring question of sex-biased dispersal: Paul J. Greenwood's (1980) seminal contribution[END_REF].

In practice, all organisms are faced with the decision to move and spread, or to stay and try to access reproduction on their natal site. This decision entails strong fitness consequences and understanding the causes and consequences of natal dispersal is thus of prime importance in evolutionary and population dynamics studies [START_REF] Bowler | Causes and consequences of animal dispersal strategies: Relating individual behaviour to spatial dynamics[END_REF].

Capture Recapture (CR) models are now widely recognized as essential tools in many fields of biological sciences, including population dynamics studies (Thomson et al. 2009).

The main strength of the CR methodology resides in the ability to disentangle biological processes, e.g. survival [START_REF] Tavecchia | Sex-and Age-Related Variation in Survival and Cost of First Reproduction in Greater Flamingos[END_REF], dispersal [START_REF] Bennetts | Methods for estimating dispersal probabilities and related parameters using marked animals[END_REF] or reproduction [START_REF] Rivalan | Trade-off between current reproductive effort and delay to next reproduction in the leatherback sea turtle[END_REF], from observation processes. However, one limitation of most CR models is the "apparent survival" problem. When moving around, some individuals may leave the study area and settle permanently outside its boundaries. These emigrated individuals are then never recaptured and are undistinguishable from dead individuals. In this situation, the most common in CR studies, survival estimates returned by CR models correspond in fact to the probability for an individual to survive and not leave the study area, i.e. the apparent survival [START_REF] Lebreton | Modeling Survival and Testing Biological Hypotheses Using Marked Animals: A Unified Approach with Case Studies[END_REF]). This limitation is even more salient if the target of the study is to estimate animal movement itself. Since only a fraction of dispersing individuals can be recaptured inside the study area, any measure of dispersal obtained under these conditions will correspond to an "apparent dispersal", thereby underestimating true dispersal. Some solutions exist when individuals can move into and out of the study area (e.g. temporary emigration, [START_REF] Kendall | Estimating state-transition probabilities for unobservable states using capture-recapture/resighting data[END_REF][START_REF] Bird | Estimating population size in the presence of temporary migration using a joint analysis of telemetry and capturerecapture data[END_REF]breeding dispersal, Oro & Doak, 2020). However, when the goal is to study natal dispersal, a single event by definition, these methods do not apply. Even in the case of capture-recapture studies where multiple sites are monitored [START_REF] Lebreton | Multistate recapture models: Modelling incomplete individual histories[END_REF], dispersal estimates only refer to transitions between sites and all individuals dispersing elsewhere will be confounded with dead individuals if no additional information is considered (e.g. dead recoveries or count surveys from other sites; [START_REF] Péron | Studying dispersal at the landscape scale: Efficient combination of population surveys and capture-recapture data[END_REF]. It is then crucial to deal with the issues of "apparent survival" and "apparent dispersal" if one wants to study natal dispersal using CR data.

In most CR studies, ancillary information about individual movements is readily available, as the locations of captures/detections are generally recorded in addition to the identity of the individual. This source of information was usually neglected but the emergence of spatial capture-recapture (SCR: [START_REF] Borchers | Spatially Explicit Maximum Likelihood Methods for Capture-Recapture Studies[END_REF][START_REF] Royle | A hierarchical model for spatial capture-recapture data[END_REF] together with the democratization of highly flexible softwares for Bayesian analyses using Markov Chain Monte Carlo (MCMC) methods (BUGS: [START_REF] Lunn | WinBUGS-a Bayesian modelling framework: Concepts, structure, and extensibility[END_REF], JAGS: Plummer, 2003, NIMBLE: DeValpine et al., 2017) recently motivated several attempts to use this additional information to solve the "apparent survival" problem [START_REF] Schaub | Estimating true instead of apparent survival using spatial Cormack-Jolly-Seber models[END_REF][START_REF] Ergon | Separating mortality and emigration: Modelling space use, dispersal and survival with robust-design spatial capture-recapture data[END_REF]. Whilst originally thought as a way to account for the lower detection probabilities of animals living at the fringe of a given study area, SCR models also present a great potential for studying animal movements (Bischof, Milleret, et al., 2020;[START_REF] Chandler | Spatially explicit integrated population models[END_REF].

Here, we introduce the Natal Dispersal model (ND), a multi-state capture-recapture model designed to simultaneously estimate natal dispersal probabilities, natal dispersal distance distributions and state-specific survival probabilities. Our model was motivated by the long-term study of a territorial and social mammal, the Alpine marmot (Marmota marmota).

We start by presenting the Alpine marmots life cycle -on which the model structure was basedand the Alpine marmot dataset. We then present the Natal Dispersal model and how it differs from a classical multi-state model. We continue with a simulation study to evaluate how both models (ND vs classical multi-state) perform under different natal dispersal scenarios. Finally, we apply the model to an individual-based long-term dataset on Alpine marmots to study sexand age-specific natal dispersal patterns in this species before discussing our results.

MATERIAL AND METHODS

Alpine marmot life cycle

Family groups The Alpine marmot is a hibernating, highly territorial and cooperatively breeding mammal living in family groups of 2 to 20 individuals. A family is typically composed of a dominant couple, adult (age >= 3 years old) and sub-adult (age=2 years old) subordinates, yearlings and juveniles born that year [START_REF] Allainé | Sociality, mating system and reproductive skew in marmots: evidence and hypotheses[END_REF]. Alpine marmots of the same family group share and defend a common territory (ranging between 0.9 ha and 2.8 ha) comprising main and secondary burrows [START_REF] Perrin | Socio-spatial Organization and Activity Distribution of the Alpine Marmot Marmota marmota: Preliminary Results[END_REF]. The main burrows are usually centrally located and used for communal hibernation, sleeping during the active season, and to give birth [START_REF] Perrin | Socio-spatial Organization and Activity Distribution of the Alpine Marmot Marmota marmota: Preliminary Results[END_REF]. Secondary burrows are spread over the entire territory and mainly used to escape potential predators. Family territories are consistently used over the years by the same family (although the family composition may change; see below) and are found adjacent to one another.

Reproduction

In each family, reproduction is monopolized by the dominant couple (apart from a few subordinate males reproducing through extra-pair copulation: [START_REF] Cohas | Extra-pair paternity in the monogamous alpine marmot (Marmota marmota): The roles of social setting and female mate choice[END_REF][START_REF] Ferrandiz-Rovira | Mate choice for neutral and MHC genetic characteristics in Alpine marmots: Different targets in different contexts?[END_REF][START_REF] Hackländer | Reproductive suppression in female Alpine marmots, Marmota marmota[END_REF]. Reproductive status is therefore generally confounded with dominance status in this species. Mating occurs shortly after the end of hibernation (mid-April to early May) and gestation lasts for 30 days. Dominant females give birth to a litter of 1 to 7 pups (median = 4) once a year and pups stay in the main burrow for 40 days, where mothers lactate them until weaning. Pups emerge between mid-June and mid-July and start eating on their own. Marmot juveniles remain in their natal territory at least until sexual maturity, at two years of age.

Dispersal Once they reached sexual maturity at two years of age, subordinate marmots can either stay on their natal territory and help raise subsequent pups or disperse. Dispersal usually takes place shortly after coming out of hibernation (Stephens et al. 2002). Subordinates that stay in the family group may become dominant by inheriting the position in their natal territory following the death of the same-sex dominant while dispersers have to evict the same-sex dominant from another territory [START_REF] Dupont | Litter sex composition influences dominance status of Alpine marmots (Marmota marmota)[END_REF]. A third possibility to access a dominant status is the creation of a new territory by a couple of dispersing individuals, but this is highly unusual in the saturated population under study (only 3 occurrences in 25 years of the study).

Once a subordinate disperses, it can either succeed or fail in securing a dominant position, but it never joins a new family group as subordinate. In case of failure in reaching a dominant status during the active season, the disperser becomes a floater, i.e. a wandering individual forced to hibernate alone and consequently exposed to very high risks of mortality [START_REF] Magnolon | Dispersion natale chez la Marmotte Alpine (Marmota marmota). Modalités et effets de quelques facteurs proximaux[END_REF]. Once the dominant position is secured, a marmot can remain dominant for several years until death or eviction by a new incomer [START_REF] Lardy | Mate change in a socially monogamous mammal: Evidences support the "forced divorce" hypothesis[END_REF]. When evicted, dominant individuals also become floaters subject to a very high risk of mortality [START_REF] Grimm | Modelling the role of social behavior in the persistence of the alpine marmot Marmota marmota[END_REF] and they usually never become dominant elsewhere. Dispersal in the Alpine marmot is therefore exclusively natal dispersal. Hence, we did not consider the possibility for evicted dominants to become dominant elsewhere in this study (fig. 1).

Because Alpine marmots are highly territorial, any intrusion in the territory is met with strong aggressiveness from the same-sex dominant and affiliated subordinates. The cost of dispersal is thus expected to strongly increase with each new territory crossed. Because territories are variable in size, two marmots travelling the same distance may be subjected to very different cost levels. In addition, Alpine marmots have the physical ability to travel long distances and through unsuitable terrain but are rarely observed further than a few hundred meters from their natal territory. For these different reasons, it is thought that the main factor driving dispersal distance of Alpine marmots is the number of agonistic encounters, i.e. the number of territories crossed, rather than landscape characteristics or physiological capacities. Hence, we measured natal dispersal distances as the minimum number of territories crossed between the natal territory and the territory where the marmot became dominant.

Alpine marmot dataset

The Alpine marmot population from the Grande Sassière Nature Reserve (French Alps, 45º29'N, 65º90'E) has been studied since 1990. Between 10 and 35 family groups have been monitored each year using capture-mark-recapture and daily observations. Marmot territories are situated at the bottom of a small valley on either side of a touristic trail crossing the lower part of the reserve. The study area covers an average of 75 ha (1.5 km long and 0.5 km wide), and is situated at approx. 2,340 m.a.s.l. The valley is surrounded by high altitude mountain peaks (> 3,500 m.a.s.l.) and steep slopes on both sides. Whilst the flat part of the valley is suitable for the establishment of Alpine marmot territories, steeper parts are unsuitable given the thickness and instability of snow-packs, rocks and the quasi-absence of vegetation [START_REF] Allainé | Habitat preferences of Alpine marmots, Marmota marmota[END_REF][START_REF] López | Explaining the successful introduction of the alpine marmot in the Pyrenees[END_REF]fig.2).

Marmots were captured each year continuously from mid-May to mid-July, i.e., after most of the dispersal events took place, using two-door live traps baited with dandelions (Taraxacum densleonis) placed near the entrances of the main burrows in the middle of each territory. Since Alpine marmots are highly territorial and only leave their natal territory to disperse, placing traps near the entrances of the main burrows allowed us to assign captured individuals to their family group with great confidence. In addition, field observations during the whole monitoring season allowed us to confirm that marmots were correctly assigned to a family territory. During the juvenile emergence period, territories are carefully observed to identify their emergence date. Juveniles are then counted and captured either using small two-door live traps or by hand.

With this procedure, virtually all juveniles in the study area were captured within three days of their emergence from the burrows.

Once captured, individuals were placed in an opaque bag before being anesthetized with Zolétil 100 (0.1 mL kg -1 ). They were then sexed, aged from their size (up to 3 years) if the exact age was unknown, and their social status was determined based on scrotal development for males and teats development for females. All individuals were individually marked using both a transponder chip injected under the skin of the neck and a numbered metal ear-tag (right ear for females and left for males to facilitate later observation). An additional colored plastic ear tag was placed on the opposite ear for dominant individuals. Captured marmots were also marked on the back with temporary hair dye to help identify them during subsequent observations in the following 7 to 15 days and verify their dominance status as well as their assignment to a family territory. When tranquillized marmots recovered, they were returned to their territory (see [START_REF] Cohas | Extra-pair paternity in the monogamous alpine marmot (Marmota marmota): The roles of social setting and female mate choice[END_REF] for more details on protocol).

We restricted the analysis to the data collected between 1990 and 2015 leading to a dataset composed of 1,270 individuals (577 females / 693 males). During this period, 170 dispersal events (69 females / 101 males) were recorded within the study area with a maximum dispersal distance of five territories crossed.

Multi-State Capture-Recapture models

Based on the life cycle described in section 1, we built two different multi-state models. First, a classical multi-state capture-recapture model which estimates only apparent survival and dispersal (the Apparent Dispersal model; AD). And second, a custom multi-state model that integrates the information available about individual dispersal movements to estimate true survival and dispersal (the Natal Dispersal model; ND). Both models are constructed around three different observable states; a subordinate state S (encompassing juveniles, yearlings and subordinate adults, i.e. pre-dispersal states), a local breeder state LB (individuals that became dominant on their natal territory by inheritance, i.e. a philopatric state), and a disperser breeder state DB (individuals that became dominant after dispersal), and one non-observable state, dead or outside the study area † (fig. 1).

Multi state capture-recapture models can be described as hierarchical models composed of two sub-models; an ecological process model describing how individuals change states between consecutive capture occasions, and an observation process model describing how individual observations collected each year are related to their underlying state. As our focus was on estimating demographic parameters and not population size, we used models that were conditional on first capture.

a. Apparent dispersal model

At first capture, individuals can be in any of the three live states (S, LB and DB) so that zifi the state of individual i at the time of its first capture fi follows a categorical distribution:

𝑧 𝑖𝑓 𝑖 ~ dcat(α) eqn.1
where α is a probability vector of length 3 with a Dirichlet distribution. The Dirichlet distribution enforces the constraint ∑ 𝛼 = 1. From this initial state, each individual has the possibility to transition to the other states or remain in its current one each year. Each year a subordinate individual may survive (probability 𝛷 ′𝑆 ), stay on its natal territory (probability 1 -𝛥 ′ ), and remain a subordinate (probability 1 -ℎ) with overall probability 𝛷 ′𝑆 . (1 -𝛥 ′ ). (1ℎ). It can also survive and inherit a dominant position on its natal territory (transition to state LB) with probability 𝛷 ′𝑆 . (1 -𝛥 ′ ). ℎ, disperse and reach a dominant position elsewhere within the study area (transition to state DB) with probability 𝛷 ′𝑆 . 𝛥 ′ , or die or disperse outside the study area with probability 1 -𝛷 ′𝑆 (transition to state †). Note that here, 𝛷 ′𝑆 is the probability that a subordinate individual survives and does not leave the study area between two consecutive occasions, 𝛥' is the probability that a subordinate individual disperses and becomes dominant within the study area, and h is the probability that a not-dispersing subordinate inherits a dominant position on its natal territory between two consecutive capture occasions.

Note also that we assume the different events (survival, dispersal, and inheritance) to be independent. As we did not consider the possibility for an evicted dominant to survive and become dominant elsewhere, a locally recruited breeder can only survive and remain in state LB with probability 𝛷 𝐿𝐵 or die with probability 1 -𝛷 𝐿𝐵 . Similarly, a disperser breeder can either survive and remain in state DB with probability 𝛷 𝐷𝐵 or die with probability 1 -𝛷 𝐷𝐵 .

Finally, a dead individual remains in its current state † with probability 1. These different state transition probabilities can be represented in a so-called transition matrix where each row corresponds to the departure state and each column to the arrival state of an individual:

𝛺 = [ 𝛷 ′𝑆 (1 -𝛥 ′ )(1 -ℎ) 𝛷 ′𝑆 (1 -𝛥 ′ )ℎ 𝛷 ′𝑆 𝛥′ 1 -𝛷 ′𝑆 0 𝛷 𝐿𝐵 0 1 -𝛷 𝐿𝐵 0 0 𝛷 𝐷𝐵 1 -𝛷 𝐷𝐵 0 0 0 1 ]
The state of an individual a given year can then be modelled as a random draw from a categorical distribution with probability vector corresponding to the individual's state the previous year:

𝑧 𝑖𝑡 ~ dcat(𝛺[𝑧 𝑖𝑡-1 , 1: 4]) eqn.2
Like the transition matrix, we can represent the different detection probabilities linking the state of an individual (in rows) to the different possible observations (in columns) at each capture occasion in an observation matrix (Θ). In our situation, we assumed that we were able to determine the state of an individual with certainty when captured so that the different possible observations were: an individual was captured and identified as a subordinate (1), an individual was captured and identified as a local breeder (2), an individual was captured and identified as a disperser breeder (3) or the individual was not captured (4):

𝛩 = [ 𝑝 𝑆 0 0 1 -𝑝 𝑆 0 𝑝 𝐿𝐵 0 1 -𝑝 𝐿𝐵 0 0 𝑝 𝐷𝐵 1 -𝑝 𝐷𝐵 0 0 0 1 ]
where p x is the recapture probability of an individual in state x. The observation recorded for a given individual in a given year, conditional on the underlying individual state, can then be modelled as a realization from a categorical process with probability vector corresponding to the individual's underlying state:

𝑦 𝑖𝑡 ~ dcat(𝛩[𝑧 𝑖𝑡 , 1: 4]) eqn.3 b. Natal dispersal model
In the AD model presented above, both the survival and dispersal parameters correspond to apparent probabilities since all individuals dispersing outside the study area cannot be observed and are therefore un-distinguishable from dead individuals. We can however reformulate this model by expressing the apparent dispersal probability (𝛥') as 𝛥*r, the product of the dispersal probability (𝛥) and the probability to remain in the study when dispersing (r). We can then construct a new model with a modified state-transition matrix (𝛺):

𝛺 = [ 𝛷 𝑆 . (1 -𝛥). (1 -ℎ) 𝛷 𝑆 . (1 -𝛥). ℎ 𝛷 𝑆 . 𝛥. 𝑟 (1 -𝛷 𝑆 ) + 𝛷 𝑆 . 𝛥. (1 -𝑟) 0 𝛷 𝐿𝐵 0 1 -𝛷 𝐿𝐵 0 0 𝛷 𝐷𝐵 1 -𝛷 𝐷𝐵 0 0 0 1 ]
Under this formulation, 𝛥 is the true dispersal probability i.e. the probability to become dominant after dispersal (in or out of the study area) and 𝛷 𝑆 is the true survival probability.

However, based on classical CR data, only the product as 𝛥*r is identifiable and we need additional information to estimate the probabilities to remain in the study area when dispersing r, and consequently 𝛥 and 𝛷 𝑆 .

The probability to settle within the study area when dispersing is necessarily individualspecific as it depends on i) the location of the natal territory an individual is dispersing from, ii) the dispersal direction, and iii) the distance travelled during dispersal [START_REF] Gilroy | A new approach to the "apparent survival" problem: Estimating true survival rates from mark-recapture studies[END_REF]).

If we consider dispersal as homogeneous in direction, i.e. individuals have the same probability to choose any direction when dispersing, the probability to remain in the study area when dispersing can be calculated for each territory s and each potential dispersal distance d as:

𝑟 𝑠𝑑 = 𝑛 𝑠𝑑 𝑁 𝑠𝑑 eqn.4
where nsd and Nsd are the number of territories inside the study area and the total number of territories situated at a distance d from the natal territory s respectively. In other words, the probability to remain in the study area for an individual dispersing a distance d is equal to the proportion of territories situated at a distance d from its natal territory s where recapture is possible. Note that d refers to the minimum number of territories crossed between the natal and settlement territories, but the same method can be applied to more classical distance measurements such as the Euclidean distance between the centroids of territories.

In our situation, the location of the natal territory s was known for all individuals born within the study area and the dispersal distance d was known for individuals that dispersed and were recaptured within the study area but had to be inferred for i) individuals that dispersed and settled outside the study area and ii) individuals that dispersed and settled inside the study area but were never recaptured. Since we measured dispersal distances as the number of territories crossed between the natal and settlement territories, we modelled them as following a Poisson distribution to which we added 1:

d'i ~ Poisson(τ) eqn.5 di = d'i +1 eqn.6
The addition of 1 was used to enforce that all dispersers moved at least one territory away from their natal one. Note that if dispersal distances are measured in a different spatial unit, and/or depending on the species' characteristics, a different distribution can be used to model dispersal distances, with the only constraint that it is strictly positive. Note also that the parameters of the chosen distribution may be functions of individual characteristics to investigate potential drivers of natal dispersal distances (see Alpine marmot example below).

Simulation study

To compare the two models presented above, we simulated different CR datasets with increasing mean dispersal distances τ from 0 to 12 territories with steps 0.5. We picked values of the other parameters in the model based on our knowledge of the Alpine marmot [START_REF] Berger | Agespecific survival in the socially monogamous alpine marmot (Marmota marmota): Evidence of senescence[END_REF][START_REF] Dupont | Litter sex composition influences dominance status of Alpine marmots (Marmota marmota)[END_REF][START_REF] Ferrandiz-Rovira | The role of the major histocompatibility complex in the wild: The case of the Alpine marmot (Marmota marmota)[END_REF]. We set the true survival probabilities to be equal and high for both breeder states (Φ DB = Φ LB = 0.95), and slightly lower for subordinates (Φ S =0.85). We also considered that recapture probabilities differed between dominants and subordinates (p DB = p LB = 0.65 and p S = 0.95). Dispersal probability Δ was set to 0.4 and inheritance probability h to 0.12. To avoid potential confounding effects due to the shape of the study area, it was randomly generated by sampling 35 cells from a 10 × 6 grid for each simulated data set, i.e. 35 territories where recaptures are possible were randomly selected.

We simulated the dataset for 20 years of captures. Each year, with the exception of the last one, 40 individuals were newly captured leading to a total of 760 individuals in each simulated dataset. We randomly assigned the 760 individuals to one of the 35 territories and individual dispersal distances were sampled following eqn. 5 and 6 with the chosen value for τ. We assumed that the landscape was homogeneous, i.e. the habitat was suitable for settlement anywhere inside or outside the study area. Based on this map, a matrix of the probabilities to remain in the study area when dispersing for all territories and dispersal distances was constructed following eqn.4 (see also Supplementary Material 2). State-transition and observation matrices similar to those in the ND model were then filled with these different probabilities. Individual capture histories were generated using eqn.3 after sampling individual states using eqn.2. Dispersal distances were kept for individuals that dispersed and settled inside the study area only to mimic true capture-recapture data. We repeated the simulation process 100 times for each value of τ, leading to a total of 2,500 simulated data sets. We then fitted both the AD and ND models to all simulated datasets and compared their parameter estimates.

Alpine marmot analysis

To account for the specifics of the Grande Sassière reserve monitoring, the analysis of the Alpine marmot dataset incorporated additional features.

First, the number of family territories monitored each year varied over the course of the study between 10 and 35. Consequently, the probability for an individual to emigrate outside the study area and the probability to observe a dispersal event inside the study area also varied between years. To account for this change in the study area size and shape, we used annual maps representative of the territories' relative positions for each year.

Second, Alpine marmot territories are variable in size and disposition, meaning that the number of territories situated at a distance d can vary depending on the territory taken as reference. The matrix of territory-specific probabilities to remain in the study area when dispersing a distance d (rsd) was thus calculated by hand based on the same annual maps of the territories.

Third, the Grande Sassière study area is surrounded by high altitude mountain tops, where habitat is unsuitable for Alpine marmots. The homogeneous landscape and dispersal assumptions were therefore unrealistic. We have shown in a supplementary analysis that the proportion of suitable habitat outside the study area strongly influences both survival and dispersal estimates (see Supplementary Material 3), and we therefore discarded unsuitable areas surrounding the study area from the potential settlement territories (fig. 2).

Fourth, following results from previous studies on the same population (Dupont et al. by running three independent MCMC chains of 10,000 iterations each with a burning period of 6,000 iterations, leading to 12,000 posterior MCMC samples per simulation. For the Alpine marmot analysis, we ran three chains of 15,000 iterations each with a burning period of 10,000 iterations, leading to a total of 15,000 posterior MCMC samples. We assessed convergence for all models by looking at parameter traceplots and calculating the potential scale reduction factor 𝑅 ̂ [START_REF] Gelman | Inference from Iterative Simulation Using Multiple Sequences[END_REF].

RESULTS

Simulation study

Based on the visual inspection and 𝑅 ̂values (all parameters with 𝑅 ̂< 1.1), we considered all simulation runs to have reached convergence. Both the AD and ND models returned unbiased estimates of the recapture probabilities (AD: posterior mean ± se; p S = 0.95 ± 0.01 and p DB = p LB = 0.65 ± 0.02; ND: p S = 0.95 ± 0.01 and p DB = p LB = 0.65 ± 0.02), inheritance probability (AD: h = 0.12 ± 0.10; ND: h = 0.12 ± 0.12) and dominant survival probabilities (AD: Φ IB = Φ LB = 0.95 ± 0.01; ND: Φ IB = Φ LB =0.95 ± 0.01) under all scenarios.

Subordinate survival and dispersal estimates, on the other hand, differed considerably between models. The AD model consistently under-estimated both parameters and this negative bias increased with mean dispersal distance to reach a lower plateau when the mean dispersal distance was larger than the maximum length of the study area dmax (fig. 3). The estimated apparent dispersal probability Δ' was between 26% and 99% lower than the simulated dispersal probability when the mean dispersal distance increased from 0 to 10 territories. In the same time, the apparent subordinate survival Φ' S was estimated from15% to 40% lower than the simulated value. Additionally, when the mean simulated dispersal distance was over dmax, the mean apparent dispersal estimates tended towards zero.

The pattern was different with the ND model. Although the difference between the simulated and estimated values of both parameters also increased with increasing dispersal distances, it was consistently lower than for the AD model. The mean parameter estimates were at most 84% and 30% lower than the simulated values for the dispersal probability and subordinate survival probability respectively. Most importantly, the dispersal and subordinate survival estimates returned by the ND model were virtually unbiased (< 1%) for simulated dispersal distances up to 0.5 dmax = 5 territories and relatively small (< 10%) up to 0.75 dmax.

By contrast, the negative bias in both dispersal and subordinate survival was at least 10% with the AD model.

The ND model was also able to accurately estimate the mean dispersal distance τ for a large range of simulated values (fig. 4). However, we found that τ estimates decreased drastically when the simulated value reached dmax.

Sex-specific dispersal in the Alpine marmot

The 95% confidence intervals of the ND posterior estimates indicated that recapture probabilities varied with time and decreased with age. Differences in recapture probabilities (on the logit scale) between yearlings and two-years and between yearlings and adults were Δptwo = -1.40 [-1.91; -0.94] and Δpad = -1.90 [-2.37; -1.50] respectively. Recapture probabilities were also higher for females than for males (difference in recapture probabilities between females and males Δpsex = -0.28 [-0.55; -0.02]). Mean dispersal distance was very similar between males and females (τfemales = 1.69 [1.55; 1.85] and τmales = 1.71 [1.56; 1.86] 5c andtable1). Because inheritance is conditional on survival and dispersal, this indicates that females in their third year of life had more chances to inherit a dominant position than males. To help understand these results, figure 6 presents the fate of a cohort of subordinate individuals as the proportions of individuals in each state at each age as predicted by the model (fig. 6).

Survival probabilities of young dominants, i.e. two-year age class, were higher than those of adults but no sex difference was evidenced for both locally recruited and immigrant breeders (table 1). Additionally, no difference between locally recruited and immigrant breeders was evidenced both for two-year-old (ΔΦsub LB = 10.43 [-34.32; 64.69] on a logit scale) and older dominants (ΔΦad LB = 0.14 [-0.36; 0.63] on a logit scale).

DISCUSSION

Simulation study

With this study, we showed how using the spatial information associated with individual detections in a multi-state capture-recapture model can allow the estimation of unbiased natal dispersal and survival probabilities even when a significant proportion of the population settles outside the study area after dispersal. We also showed how the magnitude of the underestimation of these demographic parameters by classical multi-state models depends on the mean dispersal distance (relative to the study area size). The Natal Dispersal model presented here uses the information available on dispersal movements observed within the study area and the location of individual birth territories relative to the boundaries of the study area to estimate individual-specific probabilities to leave the study area when dispersing and discriminate between mortality and emigration. Using a single model to simultaneously estimate survival and natal dispersal patterns presents the double benefit of i) being able to estimate true instead of apparent survival and dispersal probabilities, and ii) account for the imperfect detectability inherent in most field studies while rigorously propagating the resulting uncertainty to all parameter estimates.

The modularity of our approach also makes it flexible enough to be widely applicable.

First, the matrix of probabilities to remain in the study area can be easily derived for any shape and discretization of the study area, whether it is based on territories, as in this study, or regular grid cells of the desired resolution. Second, the model for individual dispersal distances is easily modified to accommodate species-specific dispersal patterns through the use of other distributions (e.g. exponential or log-normal; [START_REF] Ergon | Separating mortality and emigration: Modelling space use, dispersal and survival with robust-design spatial capture-recapture data[END_REF]. It can also easily account for variation in dispersal patterns amongst individuals (e.g. sex-specific dispersal distributions in the Alpine marmot example). These reasons make our approach especially suitable to study natal dispersal in various biological systems including other rodent (e.g.

Eurasian beaver (Castor fiber), [START_REF] Mayer | When to leave: the timing of natal dispersal in a large, monogamous rodent, the Eurasian beaver[END_REF], carnivore (e.g. African wild dogs (Lycaon pictus), [START_REF] Behr | When to stay and when to leave? Proximate causes of dispersal in an endangered social carnivore[END_REF] or bird species (e.g. Egyptian vultures (Neophron percnopterus), [START_REF] Serrano | Phenotypic and environmental correlates of natal dispersal in a long-lived territorial vulture[END_REF].

However, this approach still has some limitations. Notably, and as mentioned above, the performance of the model depends on the ratio of the mean dispersal distance relative to the size of the study area, i.e. on the proportion of dispersal events observable within the study area. When the mean dispersal distance becomes too large, most dispersers leave the study area and the model fails to estimate the correct dispersal distance distribution. In the extreme case, the dispersal distance is much larger than the study area and all dispersers settle outside. In such situation, the model has no information to rely on to estimate dispersal distances and considers that all dispersers are dead instead. In such situation, the choice of prior becomes crucial as in our simulation study, where the choice of a gamma distribution with mass close to 0 led to underestimating the mean dispersal distance. Checking the sensitivity of the estimates to the choice of prior for the dispersal distance parameter is thus highly recommended. As noted by [START_REF] Schaub | Estimating true instead of apparent survival using spatial Cormack-Jolly-Seber models[END_REF]: "if censoring becomes too strong (dispersal distances very large relative to the size of the study area), [the] model is not successful anymore in correctly estimating dispersal and consequently survival. But it is hard to imagine that any model would succeed in obtaining meaningful estimates in this case.". It is thus crucial when studying dispersal to ensure that the dataset at hands contains enough observed dispersal events to accurately estimate dispersal distances. According to our simulation results, the study area should be twice as long as the mean dispersal distance to ensure good estimation of the dispersal pattern.

Another aspect worth-considering is that the efficiency of the model to produce unbiased natal dispersal and survival estimates likely depends on the assumed dispersal distance distribution being correct, which raises the problem of selecting a good dispersal model. In our example, we used a Poisson distribution, which seems like an adequate choice for the Alpine marmots, based on the vast majority of short-distanced and low variance in dispersal events we observed (at least much shorter than the size of the study area). However, it is not possible to rule out the possibility that some individuals display a different dispersal strategy and move much longer distances. In that case, a bimodal distribution or a distribution with a longer tail would probably be a better fit although it could not be detected from the data [START_REF] Turek | Efficient estimation of large-scale spatial capture-recapture models[END_REF]. Alternative dispersal distance distributions, associated with specialized goodness of fit tests (such as Bayesian p-values, e.g. [START_REF] Ergon | Separating mortality and emigration: Modelling space use, dispersal and survival with robust-design spatial capture-recapture data[END_REF] have the potential to help identify a good dispersal distribution, although the power of such analysis to identify the correct dispersal distance distribution based on such data is still an open question.

Finally, our second simulation study stresses out the importance to correctly specify the available habitat both within and outside the study area. If no information is available about the suitability of the surrounding habitat and one assumes dispersers may settle anywhere outside the study area, inferences about the dispersal pattern, and therefore about survival, may be inaccurate (see Supplementary Material). It is therefore crucial to correctly characterize the habitat suitability outside the study area, either by using prior knowledge (as in our study case)

or through a more formal analysis (e.g. habitat selection studies).

Other studies tackled the question of apparent survival in a capture-recapture context [START_REF] Ergon | Separating mortality and emigration: Modelling space use, dispersal and survival with robust-design spatial capture-recapture data[END_REF][START_REF] Gilroy | A new approach to the "apparent survival" problem: Estimating true survival rates from mark-recapture studies[END_REF][START_REF] Schaub | Estimating true instead of apparent survival using spatial Cormack-Jolly-Seber models[END_REF]. [START_REF] Gilroy | A new approach to the "apparent survival" problem: Estimating true survival rates from mark-recapture studies[END_REF] proposed a modelling approach in which they first estimate a dispersal kernel for the study species from which they derive emigration probability for each marked individual to adjust our approach, spatial capture-recapture models present the benefit of being able to infer the movement of individuals between any two consecutive occasions thereby making it theoretically possible to also study breeding dispersal. The downside of this approach is that modelling individual locations at each time step comes at a big computational cost for already computationally intensive models (but see Bischof et al., 2020;[START_REF] Milleret | A local evaluation of the individual state-space to scale up Bayesian spatial capture-recapture[END_REF][START_REF] Turek | Efficient estimation of large-scale spatial capture-recapture models[END_REF].

Sex-specific dispersal in the Alpine marmot

The application of the Natal Dispersal model to the Alpine marmot dataset shed new light on the dispersal process in this species. Notably, we obtained estimates of adult subordinate survival much higher than in previous CR studies of the same population (0.99 against 0.53 in [START_REF] Dupont | Litter sex composition influences dominance status of Alpine marmots (Marmota marmota)[END_REF]. We also estimated the probability to become dominant through dispersal to be higher than previous estimates even though no difference was made between inheritance and dispersal (0.58 against 0.45 in [START_REF] Rézouki | Socially mediated effects of climate change decrease survival of hibernating Alpine marmots[END_REF]. In addition, we estimated the probability to become dominant through dispersal to not differ between sexes contrary to the inheritance probability where we estimated that females had three times as much chances of inheriting the dominant position than males. Based on the estimated mean dispersal distance for both sexes (τfemales = 1.69 [1.55; 1.85] and τmales = 1.71 [1.56; 1.86]) relative to the size of the study area, and the results from the simulation study, we are confident that these estimates are not underestimating the true survival and dispersal probabilities of subordinate Alpine marmots. While we observed a decrease in estimated dispersal distance for simulated dispersal distances larger than the size of the study area in our simulation study, this was accompanied by a decrease in both survival and dispersal probabilities; a pattern we did not observe in the Alpine marmot. Finally, and as expected from the results of the simulation study, the other parameter estimates (i.e. recapture probabilities, juvenile survival and dominant survival) were consistent with previous capture-recapture studies on the same population [START_REF] Cohas | Age-specific effect of heterozygosity on survival in alpine marmots, Marmota marmota[END_REF][START_REF] Dupont | Litter sex composition influences dominance status of Alpine marmots (Marmota marmota)[END_REF][START_REF] Berger | Agespecific survival in the socially monogamous alpine marmot (Marmota marmota): Evidence of senescence[END_REF][START_REF] Rézouki | Socially mediated effects of climate change decrease survival of hibernating Alpine marmots[END_REF]. These results suggest that the vast majority of subordinate individuals are in fact able to become dominant if they survive their two first winters and that the high mortality reported for subordinate sub-adults and adults in previous studies was almost entirely due to dispersal outside the study area (fig. 6).

Contrary to what is generally reported in the literature for mammals, neither the proportion of dispersers (i.e. the age-specific dispersal probability) nor the dispersal pattern (i.e. the mean dispersal distance) differed between sexes [START_REF] Dobson | The enduring question of sex-biased dispersal: Paul J. Greenwood's (1980) seminal contribution[END_REF]. This absence of sexspecific dispersal pattern is not so surprising for a monogamous and monomorphic mammal such as the Alpine marmot [START_REF] Lukas | Cooperative breeding and monogamy in mammalian societies[END_REF]. However, this absence of sex difference in dispersal probability should be considered with caution, given the formulation of our model. Because dispersal probability only reflects successful dispersal, i.e. individuals that dispersed and became dominant in another territory, it does not reflect the dispersal "decision", i.e. the probability that an individual leaves its natal territory, independently of the success of this dispersal. This probability cannot be estimated using our model because no information was available about dispersing individuals that failed to become dominant (i.e. floaters). These unsuccessful dispersers disappeared from the dataset and most likely died in the process [START_REF] Lardy | Mate change in a socially monogamous mammal: Evidences support the "forced divorce" hypothesis[END_REF]. The slight difference in survival of two-year-old subordinates males and females (and the difference in the proportion of dead males and females; fig. 6) is in accordance and tends to indicate that a greater number of males disperse after their second winter compared to females. This would also indicate a lower probability of success and therefore a higher cost of dispersal for males. This alleged male-biased natal dispersal, although unexpected from theory for a monogamous mammal, might be explained in Alpine marmots by the possibility of extrapair paternities. Although reproduction is largely suppressed in subordinate Alpine marmots [START_REF] Arnold | Reproductive suppression in male alpine marmots[END_REF][START_REF] Hackländer | Reproductive suppression in female Alpine marmots, Marmota marmota[END_REF], dominant females have the possibility to mate with multiple males, including subordinates from other families and [START_REF] Cohas | The genetic similarity between pair members influences the frequency of extrapair paternity in alpine marmots[END_REF] found that most extra-pair paternities in the Grande Sassière population were the result of male dispersers (80% of litters containing extra-pair young). This possibility for subordinate males to access reproduction through extra-pair copulation, even when failing to secure a dominant position, suggests that the expected fitness of a dispersing male is likely to be higher than that of a female. This could explain why more male subordinates disperse despite higher levels of competition among males to disperse and access reproduction, as suggested by the lower survival of two-year-old subordinate males. The higher proportion of females inheriting the dominant position from their mother (fig. 6) also suggests female-biased philopatry (or equivalently male-biased dispersal). However, further research on factors triggering dispersal and dispersal costs in this species is needed to confirm these interpretations.

To determine whether males truly dispersed more than females would require taking into account the success of the dispersal event for each individual, i.e. disentangle "natural mortality" from "dispersal related mortality". Such model requires additional information about individuals during dispersal or about their fate. Potential data types useable for this purpose are GPS-collar data and dead recoveries. Incorporating information about when and where individuals go after leaving their territories or where and when they were found dead would potentially allow to estimate what proportion of the dispersing individuals dies in the process. Unfortunately, these solutions are unlikely to work for Alpine marmots as marmots are almost always never recovered (they are either eaten by predators or die inside their burrows) and GPS methods are generally unfit for ground-dwelling species.

To conclude, the ND model we presented here is a new step towards a more complete and realistic assessment of natal dispersal and its drivers for which biologists of many fields have a great interest. As all models, it presents some limitations, but, with a clear understanding of its underlying assumptions, it can contribute to increasing the knowledge about this central phenomenon in population biology that is natal dispersal.

  apparent survival estimates obtained from a separate capture-recapture model. Estimating the dispersal distance distribution within the capture-recapture model as we did here presents the advantage to allow for a better propagation of uncertainty in the different parameters of the model. Ergon & Gardner (2014) and Schaub & Royle (2014) also used a framework that does not require to separately estimate individual emigration probabilities outside the capturerecapture model. Instead, they used open-population spatial capture-recapture models, which explicitly model the location of each individual at each time step, even if outside the study area, thus only implicitly modelling individual emigration probabilities. Whilst relatively similar to

Figure 1 .

 1 Figure 1. Schematic representation of the Alpine marmot's life cycle as used in the natal

Figure 2 .

 2 Figure 2. 3D (a) and schematic map (b) of the Grande Sassière nature reserve and marmot

Figure 3 .

 3 Figure 3. Mean subordinate survival (Φ S ) and dispersal probabilities (Δ) estimates returned by

Figure 4 .

 4 Figure 4. Simulated and estimated dispersal distance for the Natal Dispersal model. Solid line

Figure 5 .

 5 Figure 5. Age-and sex-specific survival Φ(a), dispersal Δ (b), and inheritance h (c)

Figure 6 .

 6 Figure 6. Fate of a cohort of Alpine marmots with age for females (red lines) and males (blue

  models • Spatially explicit capture-recapture • Sex-specific dispersal • 49

Table . 1

 . Age-and sex-specific posterior mean estimates (standard deviation) of the survival, dispersal and inheritance probabilities of Alpine marmots. Sex difference indicates if the 95% credible interval of the sex difference parameter does not overlap the value 0.

	parameter	females	males	sex difference
	Φ S juv	0.56 (0.03)	0.55 (0.02)	NO
	Φ S year	0.71 (0.04)	0.74 (0.04)	NO
	Φ S	two	0.97 (0.05)	0.84 (0.11)	NO
	Φ S ad	0.99 (0.01)	0.99 (0.01)	NO
	Φ LB two	0.99 (0.04)	0.99 (0.07)	NO
	Φ LB ad	0.78 (0.03)	0.70 (0.04)	NO
	Φ IB two	0.99 (0.06)	0.99 (0.10)	NO
	Φ IB ad	0.80 (0.03)	0.77 (0.03)	NO
	Δjuv	0.00 (0.00)	0.00 (0.00)	NO
	Δyear	0.10 (0.03)	0.11 (0.03)	NO
	Δtwo	0.58 (0.06)	0.57 (0.06)	NO
	Δad	0.58 (0.06)	0.58 (0.05)	NO
	hjuv	0.00 (0.00)	0.00 (0.00)	NO
	hyear	0.01 (0.01)	0.01 (0.01)

NO htwo

0.

32 (0.06) 0.09 (0.04)

				YES
	had	0.59 (0.09)	0.58 (0.07)	NO

Acknowledgements:

The authors would like to thank all students involved in the trapping of Alpine marmots and authorities of the Vanoise National Park for allowing us to work in the Grande Sassière Nature Reserve. We would also like to thank T. Ergon and one anonymous reviewer for their helpful comments on a previous version of the manuscript.

Data availability statement: data, R, JAGS and NIMBLE scripts to reproduce the simulation study and marmot analysis are available from the GitHub repository: https://github.com/PierreDupont/marmotDispersal

2015; [START_REF] Berger | Agespecific survival in the socially monogamous alpine marmot (Marmota marmota): Evidence of senescence[END_REF][START_REF] Rézouki | Socially mediated effects of climate change decrease survival of hibernating Alpine marmots[END_REF], we considered recapture probabilities to vary with time, age and sex. In addition to the sex-effect in which we were interested, we accounted for age-specific subordinate survival and transitions probabilities by considering four age classes: juvenile (from zero to one-year-old), yearling (from one to two), subadult (from two to three) and adult (three years old and older). We also accounted for an effect of the logarithm of the number of male helpers on the juvenile survival, as it was shown to be an important driver of the juvenile survival during hibernation in this species [START_REF] Allainé | Is there an optimal number of helpers in Alpine marmot family groups?[END_REF][START_REF] Dupont | Litter sex composition influences dominance status of Alpine marmots (Marmota marmota)[END_REF]. Because no individual attained a dominant status before two years old, only the two last age classes were used for dominant individuals.

Finally, we considered sex-specific dispersal patterns, i.e. two independent Poisson distributions with sex-specific mean. The full set of model parameters was then pt,age,sex the time, age and sex-specific recapture probabilities, Φ IB age,sex and Φ LB age,sex the age and sexspecific survival probabilities of dominant individuals, Φ S age,sex the age and sex-specific survival probabilities of subordinate individuals, Δage,sex the age and sex-specific dispersal probabilities, hage,sex the age and sex-specific inheritance probabilities, ri the individual-specific probability to remain in the study area when dispersing, and τsex the sex-specific mean dispersal distance. The R and JAGS script for the Alpine marmot analysis is provided in Supplementary Material 3.

Model fitting

All models were fitted using Markov Chain Monte Carlo (MCMC) simulations with the computer program JAGS [START_REF] Plummer | JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling[END_REF]) called through R3.2.5 (R core team 2016) via the R package jagsUI (Kellner 2014). All parameter priors were chosen to be weak or un-informative.

All survival, dispersal, inheritance and recapture probabilities were given uniform priors between 0 and 1. The prior for the mean of the Poisson distribution in the ND model was given a vague gamma distribution Γ(0.0001,0.0001). All models in the simulation study were fitted Author Contributions: PD, RP and DA conceived the study. RP and PD developed the models and performed the analysis. DA, MFR and PD contributed to the data collection. PD wrote the manuscript with contribution from all authors. All authors provided editorial advice and final approval on the manuscript.

Conflict of interest:

The authors declare that they have no conflict of interest.

Supplementary material 1:


Akylas Stratigakos
email: akylas.stratigakos@mines-paristech.fr

A robust fix-and-optimize matheuristic for timetabling problems with uncertain renewable energy production

Keywords: optimize, local neighborhood search, renewable energy forecasting, robust optimization, university timetabling

This work presents a fix-and-optimize matheuristic to solve timetabling problems under uncertainty. Specifically, a combined university timetabling and electricity storage scheduling problem is considered, subject to uncertainty stemming from renewable energy production and electricity demand. The problem is formulated as a large Mixed Integer Program (MIP) and the proposed solution combines Large Neighborhood Search coupled with scenario-based robust optimization for handling uncertainty in the objective function. First, an adequate feasible schedule is derived considering only hard problem constraints, in this case scheduling of recurring lecture activities. Next, the solution is improved with a fix-and-optimize heuristic search. In each iteration, the MIP solver explores a large neighborhood by fixing a subset of variables and optimizing over the remaining free variables. The process is repeated several times until a stopping criterion is met. To address uncertainty in the objective, probabilistic scenarios are derived from interval forecasts and the worst-case energy cost is minimized. The results derived from the participation in a technical challenge show that the proposed approach provides competitive solutions relatively fast, even for large problem instances, while also hedging against large forecast errors.

I. INTRODUCTION

The integration of renewables in the generation mix is key to the ongoing decarbonization process. Their stochastic nature, however, poses an important challenge to their integration within larger operational processes. Dealing with this stochasticity typically involves a process where several prediction and optimization components are deployed in sequence. This work presents a methodology to tackle large size scheduling and timetabling problems with uncertain renewable production and electricity demand, based on the participation in the technical challenge organized by the IEEE Computational Intelligence Society (IEEE-CIS) [START_REF] Bergmeir | Ieee-cis technical challenge on predict+optimize for renewable energy scheduling[END_REF].

In this challenge, several instances of a university timetabling problem are considered, which can be modelled as a Mixed Integer Program (MIP). It is known that such problems are NP-hard, therefore, an optimal solution within reasonable time is not attainable. The proposed solution is guided by several challenges that revealed themselves during the early stages of the competition. First, the limited computational resources did not allow to solve the (multiple) problem instances to optimality. Second, the computational cost also hindered the ability to explore different strategies during the validation phase, e.g. how to tackle the parameter uncertainty. Lastly, as the time to be allocated in this challenge was also limited, the decision was made to focus on the optimization component at the expense of the prediction component. Considering the above, the proposed solution adheres to the following: i) can be implemented in a standard machine, ii) provides competitive results relatively fast, iii) and provides hedging against large forecast errors.

To this end, the described solution is based on a fix-andoptimize heuristic search to iteratively improve the solution of the MIP solver (matheuristic). Uncertainty in problem parameters (i.e., renewable production and electricity demand) is handled with scenarios based on marginal predictive intervals. A robust formulation is provided to minimize the worstcase expected cost within the provided set of scenarios, thus offering protection against miscalibrated forecasting models.

II. BACKGROUND

The individual components are briefly described in this section. In short, the problem concerns university timetabling and storage scheduling with the goal to minimize the total energy cost of a university complex, subject to uncertain electricity demand and solar production. For a detailed description of the problem, see [START_REF] Bergmeir | Ieee-cis technical challenge on predict+optimize for renewable energy scheduling[END_REF].

1) Problem formulation: The MIP formulation of the university timetabling problem is omitted due to space limitations (see [START_REF] Daskalaki | An integer programming formulation for a case study in university timetabling[END_REF], [START_REF] Lach | Curriculum based course timetabling: new solutions to udine benchmark instances[END_REF] for relevant formulations). Recurring activities (lectures) are considered "hard" constraints, as they are required for a feasible schedule solution. Once-off activities (corresponding to lab experiments) are considered "soft" constraints, as they are not required for a feasible solution. Recurring activities are repeated in a weekly basis, thus a single set of variables per week suffices. For once-off activities, the defined variables span the whole scheduling period.

In short, a binary variable is assigned for each activity and each required time period. Further, an integer variable is defined for each activity, time period, building, and room type, to determine the number of rooms assigned. Note that this formulation is somewhat excessive. However, modern MIP solvers remove redundant constraints and variables during a presolve phase, thus leading to significantly smaller problems during the branch-and-bound search. For reference, for a large instance the initial problem comprises approximately 9 • 10 3 continuous and 11 • 10 6 integer (4.5 • 10 5 binary) variables. After the presolve phase, these are decreased to 640 continuous and 2 • 10 5 integer (7.8 • 10 3 binary), respectively. To the author's opinion, the requirement of activities to be scheduled in consecutive time slots poses the greater modeling challenge, which is handled with the introduction of auxiliary variables and applying a big-M reformulation of constraints.

2) Forecasting: The uncertain electricity demand and solar production must be predicted for the whole scheduling period; alternatively one could consider the aggregated net load. A separate forecasting model is trained for each series provided. Regarding solar panels, several base models, namely Random Forest, Gradient Boosting Machines, Ridge Regression, and Local Learning Regression, are developed and aggregated via averaging for the final prediction.

Regarding buildings' demand, the prediction task is more challenging due to missing data and restrictions related to the COVID-19 pandemic. Overall, buildings 1 and 4-6 show relatively low demand, which in turn reduces their impact on final energy cost. The decision was made to discard all missing data, as most of them relate to buildings 4-6. All demand forecasts were derived with a Random Forest model. Feature data included the weather data provided by the organisers, alongside categorical variables to model the calendar effect, impact of COVID-19 restrictions, exams period, etc. For buildings 0 and 3 a pre-processing step was included to remove possible outliers. Lastly, a Quantile Regression Forests model is employed to generate probabilistic forecasts for building 3, which are subsequently used to derive a high and low impact scenario.

3) Objective: The solutions are evaluated based on aggregate energy cost, which includes a peak load tariff. A scenariobased robust objective is derived to hedge against large forecast errors. Specifically, two scenarios for aggregated net load are derived based on the marginal predictive density of demand in building 3. The objective is the reformulated in a standard epigraph form, and the worst-case cost over all scenarios is minimized. Ideally the scenarios should incorporate correlation between the different sources of uncertainty. Due to limited time and after a trial-and-error period, it was decided to derive scenarios solely based on the predictive density of building 3, which represented the largest load. Fig. 1

III. METHODOLOGY

This section presents the complete methodology that builds on the individual components discussed above. A fix-andoptimize matheuristic is proposed to solve the MIP problem [START_REF] Lindahl | A fix-and-optimize matheuristic for university timetabling[END_REF], which belongs to the class of Large Neighborhood Search (LNS) heuristics [START_REF] Pisinger | Large neighborhood search[END_REF]. The approach is method based, meaning that the underlying model is the same, but the local neighborhood to search an improved solution changes at each iteration. Note that modern optimization solvers (e.g. GUROBI) require the model to be built once and each iteration is warm-started from the previous solution by adding or removing constraints, which significantly speeds-up computations.

The solution process is described in the flowchart shown in Fig. 2. The derived predictions described in the previous Section are considered as input in this process. First, the problem is solved once considering only "hard" constraints, providing a minimum feasible solution, which translates into optimizing only over recurring activities. The sample solutions provided by the organisers are used to warm-start the solver. Next, the iterative search begins by sampling a subset of activities, considering both recurring (r # ) and once-off (a # ). Variables related to non-selected activities are fixed to their previous solutions and the subsequent problem is solved to optimize over the remaining free variables. The process is repeated several times until the maximum number of iterations is reached (max-iter) or until the solution stops improving after a number of iterations (patience) subject to a predefined tolerance (tol). Variables related to storage remain free during all iterations, as the computational burden increased only marginally. The search algorithm is detailed in Algorithm 1. Note, this search is greedy and random. Due to precedence constraints, subsequent problems with fixed Variables pertaining to the selected activities are free, the rest are fixed to previous solution 6:

Solve MIP, determine objective v new 7:

if v * -v new v * ≥ tol then Update solution, v * ← v new , count = 0 8: else count+ = 1 9:
end if iter+ = 1 10: end while 11: return: Improved solution with objective v *

IV. EXPERIMENTS

This section provides some experimental results and details on tuning the hyperparameters. During Phase 1 (validation), different strategies were explored, mainly concerning how to deal with uncertainty in the objective function. Fig. 3 provides a comparison between a deterministic formulation (i.e., using point forecasts) and the scenario-based robust approach for the Phase 1 data. Note that the validation data includes a large demand outlier, which actually obscures the impact of predictive accuracy in the final cost. Regardless the objective formulation, Fig. 3 highlights that the fix-andoptimize search significantly reduces costs, with results being more pronounced for larger problem instances.

S K D V H B L Q V W D Q F H B O D U J H B S K D V H B L Q V W D Q F H B O D U J H B S K D V H B L Q V W D Q F H B O D U J H B S K D V H B L Q V W D Q F H B V P D O O B S K D V H B L Q V W D Q F H B V P D O O B (QHUJ\&RVW 'HWHUPLQLVWLF%DVH 5REXVW%DVH 'HWHUPLQLVWLF)XOO 5REXVW)XOO
Regarding the Phase 2 (testing phase), all experiments are conducted using a standard machine with an Intel i7@2.30GHz and 32GB of RAM. The problems are modeled and solved with GUROBI using the Python-API. For reference, the required time to build a model is approximately 200-250 sec for small instances and 900-950 sec for large ones. To create the aggregated net load scenarios, the 10% and 90% quantile of the predictive marginal density of building 3 are selected. For the initial feasible solution a hard time limit of 2 hours is set, although a smaller time limit would also suffice. Subsequent iterative solutions are solved with a hard time limit of 20 min. A soft time limit of 5 min, if the optimality gap is below 0.005%, is also set. The total number of iterations is set at max-iter = 150, the patience is set at patience = 15, and the improvement tolerance is set at tol = 0.001%. Lastly, at each iteration the number of sampled activities are r # = 10 for recurring and a # = 5 for once-off, same for both small and large problem instances. The code is provided in [6]. For reference, Fig. 4 shows the evolution of the objective value for a large instance in the testing phase.

V. DISCUSSION

This section provides a brief discussion on how to improve the performance of the proposed solution, insights gained from ,WHUDWLRQV (QHUJ\&RVW the participation in the technical challenge, and promising research directions. Regarding the performance of the final submission, several improvements could be made. First, considering better informed scenarios could improve the solution. The decision to base scenarios solely on the predictive density of building 3 is somewhat arbitrary; considering scenarios that account for the correlation between different sources of uncertainty could improve final performance. Second, improving the formulation of the timetabling problem could reduce computational costs. Third, note that the performance of the iterative search procedure is bounded by the solution to the full MIP problem. Smaller problem instances could probably be solved to optimality within reasonable time, given adequate time and computational resources. Moreover, the search is greedy and random; considering precedence constraints and sampling blocks of activities instead of individual activities could also improve the performance. Overall, it is the author's opinion that the iterative search proposed offers higher utility for larger problem instances. Lastly, the optimization window changed during Phase 2 to match the forecasting window. As a result, additional effort was required to modify the optimization code to the new setting, which somewhat hindered the ability to properly evaluate results in the validation data. Avoiding modeling errors also proved important, as it hindered the ability to exploit and evaluate different solutions during within the timeframe of this challenge.

The impact of forecast uncertainty on optimization efficacy in this challenge, to the author's opinion, remains somewhat obscure. The validation data (Phase 1) included an extremely large demand outlier, which affected the peak demand and the respective tariff. In turn, this mitigated the impact of the objective formulation (deterministic versus robust). Further, examining the results on validation data (Phase 1) showed that, at least for the large instances, the peak demand tariff comprised the biggest part of total energy cost. However, the magnitude of the load to be scheduled during Phase 2 (relating to the respective activities), was significantly smaller. If the problem instances are viewed as data points from a problem distribution to be learned [START_REF] Bengio | Machine learning for combinatorial optimization: a methodological tour d'horizon[END_REF], this could be considered as a change in the underlying distribution. Overall, the peak tariff became less important during Phase 2, which to the author's opinion, somewhat obscures the impact of forecast accuracy in total costs.

Finally, increased predictive accuracy does not always translate into improved optimization performance. An interesting avenue to explore would be to train the model to directly minimize downstream optimization costs. For example, [START_REF] Wilder | Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization[END_REF] proposes training deep learning models in conjunction with combinatorial optimization problems to improve the utility of the predictions (i.e., forecast value).

VI. CONCLUSIONS

This report detailed the participation in a technical challenge organized by IEEE-CIS. Specifically, a robust fix-and-optimize matheuristic to tackle large scale timetabling and scheduling problems under uncertainty was described. Uncertainty in renewable production and electricity demand is modelled with scenarios and the worst-case energy cost is minimized. The results show that the proposed solution provides fast solutions that achieve competitive performance, while also hedging against large forecast errors.

Fig. 2 .Algorithm 1 4 :

 214 Fig. 2. Flowchart of the solution. Yellow boxes denote data input, green boxes denote forecasting models, and blue boxes denote optimization components.

Fig. 3 .

 3 Fig.3. Results on Phase 1. "Base" denotes the initial solution, "Full" denotes the final solution after the fix-and-optimize search is implemented.

Fig. 4 .

 4 Fig. 4. Progression of the local search.

  presents an indicative set of forecasts for a single day. Effectively, this approach creates a small box uncertainty set that aims to protect against larger forecast errors.

	1HW/RDG([SHFWHG
	1HW/RDG/RZ6FHQDULR
	1HW/RDG+LJK6FHQDULR
	%DVH/RDG
	N:
	+RXU

Fig. 1. One day of aggregated net load forecast and high/low scenarios for Phase 2. Base Load refers to the expected buildings demand without considering solar production.

ACKNOWLEDGMENT

The author would like to thank IEEE-CIS for organizing this very interesting technical challenge and, in particular, Dr. Bergmeir for his prompted responses to all the queries.


Anatolii V Mokshin

Jean-Louis Barrat

Shear Induced Structural Ordering of a Model Metallic Glass

Keywords: numbers: 46.35.+z, 05.70.Ln, 64.60.qe, 64.70.pe

We report results of non-equilibrium molecular dynamics simulations of a one-component glassy system under the influence of a shear flow, with the aim of investigating shear induced ordering of this system. In spite of the very low temperature, the system transforms into a strained crystalline state through well defined nucleation events. Various characteristics of the observed ordering at different strain rates and temperatures are discussed. We also define and discuss the transition rates.

Nucleation, clustering and ordering processes in supercooled melts are a subject of considerable interest, both in soft matter and in hard materials such as metallic alloys or molecular glass formers. Crystal ordering is generally thought to proceed through nucleation events, which are described within the framework of classical nucleation theory (CNT) 1,2,3,4 . In this theory, the nucleation is associated with the local crossing of a free energy barrier under the influence of thermal noise. The situation becomes more intriguing if, instead of a supercooled liquid, we consider a glassy system at very deep supercooling [START_REF]Nucleation Theory and Applications[END_REF] . The kinetic slowing down and high viscosity of a glass make the structural rearrangements, including those required for crystal nucleation, so slow that the corresponding transition is practically unobservable on the experimental time scales. Obviously, the devitrification can be accelerated through reheating of the glass. An interesting alternative, which is in fact commonly used in soft matter systems, is to utilize shear as an ordering tool rather than temperature. Crystal-crystal or amorphous crystal transformations are also known to take place in processes such as grinding and milling of molecular materials, which in general involve a complex combination of mechanical stresses and temperature variations [START_REF] Willart | [END_REF] .

Recent numerical studies provide the evidence of crystallization induced by external shear in an amorphous, single-component Lennard-Jones system 7,8 However, such single-component Lennard-Jones systems have a very weak glass forming ability and the transition to a crystalline state at low temperature proceeds through a spinodal, rather than nucleation, route 9,10 .

In this study, we consider the case of a system with a much better glass forming ability than the Lennard-Jones system, but still keeping the simplicity associated with one-component systems. This system is made of particles interacting through the short ranged, oscillating Dzugutov potential 11,12 U (r * )/ǫ = A (r * -m -B) exp c r * -a Θ(a -r * )

+ B exp d r * -b Θ(b -r * ), r * = r/σ, (1)
where σ and ǫ are the length and energy scales, respectively, the parameters A = 5.82, B = 1.28, m = 16, a = 1.87, b = 1.94, c = 1.1, d = 0.27 are taken as described in Ref. 11 and Θ(r) represents the Heaviside step function. So, the potential is characterized by a minimum U (r min ) = -0.581ǫ located at r min = 1.13σ and a maximum U (r max ) = 0.46ǫ at r max = 1.628σ. The presence of these oscillations can be considered as a simplified attempt to reproduce the screened Coulomb repulsion of ions in metals. The potential is tuned in such a way that it prevents the formation of the most common crystal phases at low pressure. At high pressures (P ≥ 5ǫ/σ 3 ) and low enough temperatures 13 , the FCC or HCP crystalline structures can be formed, although the ground state of the system at zero temperature is a BCC structure. A slow quench from a high-temperature state this model favors the formation of quasi-crystals including Frank-Kasper structures, whereas a faster cooling leads to a metastable amorphous state with a modecoupling temperature T MCT estimated to be 0.4ǫ/k B (see Ref. 14 ). A remarkable property is that this is onecomponent good glass former shares many common features with much more complex metallic glasses. With only one component, however, the study of topological reorganization, which excludes such effects as chemical ordering or demixing, is facilitated 15 .

The system consists of N = 19 652 particles contained in a cubic simulation box of the volume L 3 with L = 28.55σ. The structural changes in such a rather large system is expected to be relatively insensitive to the boundary conditions imposed by the simulation cell 16 . A set of well equilibrated independent samples is created at a temperature T = 2.0ǫ/k B . These samples are cooled at a rate ∆T /∆t = 0.000 97ǫ/(k B τ ) and the configurations corresponding to the temperatures T = 0.01, 0.03 and 0.06ǫ/k B are stored 17 . To impose the shear we create two solid walls at the ends of the simulation box in the y direction by freezing all the particles in the x-z plane over the range of three interparticle distances. The top wall moves in the x direction with the instantaneous velocity v (t) = γL y (t)e x at a constant strain rate γ and normal load P yy = 7.62ǫ/σ 3 ; L y is the distance between the walls. This value of the normal pressure is thermodynamically favorable to FCC crystalline forms and allows one to avoid the phase of quasi-crystals for the Dzugutov system 19 . Constant temperature conditions are ensured by rescaling the velocity component in the neutral z direction perpendicular to the flow x and the velocity gradient y direction. The termostatting was implemented only if the temperature has changed more than 1.5% of the considered value. Other details of the simulation protocol are identical to those described in Ref. 8 .

The local order in the system can be characterized by parameters that probe the appearance of orientational order around each atom. In this work, we use both the local orientational 20 bond-order parameter q l (i) and its global 21 counterpart Q l :

q l (i) =   4π 2l + 1 l m=-l N b (i) j=1 Y lm (θ ij , ϕ ij ) N b (i) 2   1/2 , (2a)
Q l =   4π 2l + 1 l m=-l N i=1 N b (i) j=1 Y lm (θ ij , ϕ ij ) N i=1 N b (i) 2   1/2
.

(2b) Here Y lm (θ ij , ϕ ij ) are the spherical harmonics, N b (i) denotes the number of the nearest "neighbors" of particle i, and θ ij and ϕ ij are the polar and azimuthal angles formed by radius-vector r ij and some reference system. The parameter Q l is invariant under rotations. Nonzero values of Q l for the different crystalline cluster structures including those with icosahedral symmetry appear for l ≥ 6 (see Ref. 21 ). In this work we define "neighbors" as all atoms located within a sphere of radius |r ij | < 1.5σ around an atom i. For a fully disordered system one has Q 6 → 0, whereas the ordered geometries are characterized by the next values of Q 6 :

Q f cc 6 = 0.5745, Q bcc 6 = 0.5106, Q ics 6 = 0.6633, Q hcp 6
= 0.4848 and Q sc 6 = 0.3536 (see Refs. 20,22 ). By performing a cluster analysis for the detection of the "solidlike"particles (for details see Refs. 8,20 ), we obtain the cluster size distribution n s (t) and its moments s (m) (t) = s s m n s (t) as a function of time. The zeroth moment s (0) defines the total number of clusters, the first moment s (1) represents the total number of particles involved in an ordered phase (solidlike particles), the ratio of the second and first moments, s (2) /s (1) , is associated with the average cluster size. In a finite sample we also obtain s * = max s, ns =0 [n s ] as the size of the largest cluster. The average spatial cluster size ζ can be also estimated through the distance between maximallyseparated particles within the same cluster.

We begin the discussion of our results with Fig. 1, which represents the time evolution of the orientational order parameter Q 6 (t), of the number of solidlike particles s (1) (t) and of the number of clusters s (0) (t). These curves are computed from a single run in a sample sheared with the rate γ = 0.005τ -1 at the temperature T = 0.03ǫ/k B . All these quantities are correlated and illustrate clearly the transition the initially amorphous system into an ordered state, which takes place rapidly after a system-dependent incubation time. Further, the transition is also seen from the growth in the average cluster size presented in Fig. 2. The clear transition from an amorphous phase to a strong crystalline order after an incubation time, during which the system stays amorphous, distinguishes this case from Lennard-Jonesium 8 . Such a behaviour, in which the incubation time is a statistical, sample-dependent quantity, is typical for a nucleation process. After the transition, the system is transformed essentially into a FCC structure, as is clearly evident from the snapshots and from the radial distribution of particles constitutive of the ordered structures, g cl (r) (inset of Fig. 1(a)). Both the cluster size distribution and the structure of the clusters are essentially independent of shear rate after the transition is completed, indicating that the local structure is not shear Figure 2: Time evolution of two average cluster characteristics in the sheared system with the strain rate γ = 0.005τ -1 at the temperature T = 0.03ǫ/kB . Main: Average cluster size as the ratio between the second and first moments of the cluster size distribution. Inset: spatial correlation length (see text). Relatively insignificant values of both terms for the ordered phase is consequence of the retaining of small crystalline grains (for the range t > 600τ the ratio between the largest cluster s * presented in Fig. 5 and the average cluster size s (2) /s (1) is ∼ 4 times). rate dependent. It should be noticed that ordering takes place primarily away from the walls, as is seen from the density profiles of solidlike particles presented in Fig. 3. Similar profiles were observed for all considered samples under shear and for all considered strain rates and temperatures. Such a behavior points to the homogeneous scenario of the nucleation. The evolution of the largest nucleus as a function of time is shown in Figs. 4 and5. After the transition into an ordered phase, which is accompanied by the growth of amount of solidlike particles s (1) [see Fig. 1(b)], the total number of clusters s (0) [see Fig. 1(c)] as well as the size of the largest cluster presented in Fig. 5 become essentially stationary. The system achieves a high level of structural ordering, which is, however, far from perfect. When ∼ 83 % of the bulk becomes included into an ordered phase, no further ordering is observed. Remaining particles formative the disordered ranges are located mainly near the walls. Although the statistical properties remain essentially constant, the size of the largest cluster displays considerable fluctuations (Fig. 5). The system is made of nanocrystalline grains containing at most a few hundred particles (see inset of Fig. 5), which are prevented by the shear flow from growing further but still exchange particles under the influence of the mechanical noise generated by the shear flow. Until now we have discussed the strain history of a single, typical sample. We now turn to the consideration of the statistical properties of this strain history by collecting statistics on 50 and more independent runs. The evolution of orientational bond-order parameter Q 6 (γ) (with γ = γt) at the temperature T = 0.03ǫ/k B and at various values of the strain rate γ is presented in Fig. 6(a). These curves display a marked transition at a well defined value γ n of the strain, which increases with shear rate. Obviously, the transition can be characterized by a transition rate τ -1 n , which we evaluate using the mean first-passage time (MFPT) method. This method 16,23,24 focuses on the average time τ needed for the system to reach for the first time some state with a defined value of the reaction coordinate φ → φ 1 , φ 2 , φ 3 , ..., φ ∈ R, where R is the domain of states 25 . As a rule, this allows one to evaluate the distribution of time τ (φ) over the reaction coordinate φ. For an activated process, the transition rate τ -1 n is simply determined by the inverse of the MFPT, at which a plateau is reached, i.e. τ -1 n ≡ τ (φ p ) -1 , whereas the total τ (φ) can be fitted by

τ (φ) = τ n 2 (1 + erf((φ -φ * )c)), (3)
where erf(. . .) is the ordinary error function, c is the local curvature around the top of the barrier and φ * corresponds to the critical value of reaction coordinate (for details, see Ref. 16 ).

The transition rates calculated within the framework of the MFPT approach with the order parameter as the reaction coordinate, Q 6 → φ, were used to rescale the strain dependencies presented in Fig. 6(a). So, Fig. 6(b) shows the order parameter [same with Fig. 6(a)] with the rescaled argument ξ = γ/γ n = t/τ n , where γ n = γτ n is the transition strain. As can be seen, all curves collapse onto a single master curve, which is indicative of a common mechanism for the transition scenario. The transition takes place in three stages, which, in turn, reflect the influence of the external drive on the inherent dynamics of the system. The first initial stage, where the system remains mainly in an amorphous phase, is characterized by shear-induced "unjamming" of the system and by the fluidization of the glass [START_REF] Berthier | [END_REF] . Due to the shear the diffusivity grows at this stage and the potential energy increases. The characteristic time scale associated with inherent dynamics is defined here by the relaxation time at a finite shear rate, τ α ( γ) (see Ref. [START_REF] Berthier | [END_REF] ). The second stage corresponds directly to the transition of the amorphous system into an ordered state and is caused by cluster nucleation 8 . Its appearance can be characterized by the transition time scale τ n (or the transition strain γ n ). During this stage the order parameter grows rapidly, following a power law ∝ ξ 3±0.1 . In the third stage, shear initiates rearrangements in the ordered structures and prevents the system from undergoing a complete crystallization, the order parameter increases only very slowly and probably saturates for large strains 27,28 .

The extracted values of the transition rates τ -1 n at different strain rates and three temperatures are presented in Fig. 7. At small values of γ, the transition rates increase with strain rate in a manner close to a power law. Note that the rates are practically independent of temperature at these low values of γ. However, with the subsequent increase of γ the growth of transition rates becomes slower and temperature dependent. Eventually, starting with γ ≈ 0.005τ -1 transition rates decrease with the strain rate γ.

To describe the transition of a system into an ordered phase, the ordinary practice is to use the Arrhenius-like equation

τ -1 n = ν 0 exp(-β∆E), (4) 10 -4 10 -3 10 -2 10 -3
dγ/dt (in units of τ -1 ) The solid curves are the fits using Eqs. ( 5) and ( 6). The short straight line corresponds to a unit slope.

τ n -1 (in units of τ -1 ) T=0.01ε/k B T=0.03ε/k B T=0.06ε/k B
where ν 0 is the attempt frequency (or kinetic prefactor) to overcome a barrier ∆E and β = 1/(k B T ) is the inverse thermal energy; G ≡ β∆E. Two time scales arise naturally in this equation, τ n and ν -1 0 . Identifying a rate τ -1 n with the nucleation rate 20 and the barrier ∆E with the free energy cost to form a cluster of the critical size, Eq. ( 4) is the main equation of CNT. The use of Eq. (4) for systems near thermal equilibrium is straightforward: from the measured transition rates τ -1

n at different temperatures one can easily define the thermodynamic barrier ∆E as well as the kinetic term ν 0 . However, the use of this equation for non-equilibrium cases that violate detailed balance is not generally justified, although there are some indications that a similar approach using an "effective temperature" could be relevant 29,30 .

Nucleation events in systems under shear have been identified in Ref. 31 with the conclusion that the thermodynamic barrier ∆E is an increasing function of shear rate. We propose an interpretation of our results along the same lines by making a number of reasonable assumptions, namely: (i) A formula similar to (4) can be used providing an unknown effective temperature T eff , which is introduced to replace the usual bath temperature T = 1/(k B β).

(ii) The kinetic prefactor is set by the time scale for relaxation in a system under shear, which is inversely proportional to the shear rate in this low temperature regime 32 . Hence ν 0 ∝ | γ|.

While the assumption (i) results in the generalization of the energy barrier height G = β∆E in Eq. ( 4) to the γ-dependent term 7,31 , G( γ, T ) = ∆E/(k B T eff ), assumption (ii) is motivated by the fact that the kinetic processes that can be observed in a glass are dominated by the in-fluence of the shear. As a result, Eq. ( 4) can be rewritten as

τ -1 n ∝ | γ| exp[-G( γ, T )]. (5)
This is justified by the fact that the crystal nucleation in a glassy system free from any drive, e.g. γ → 0, becomes at very deep supercooling extremely rare event, i.e. τ -1 n → 0, practically unobservable within the time-scale of the experiment as structural rearrangements are inhibited by the very high viscosity [START_REF]Nucleation Theory and Applications[END_REF] .

These assumptions are not sufficient to determine independently the effective temperature and the barrier ∆E, both of which, as can be expected, depend on shear rate and on the actual, thermodynamic temperature. However, the variation of a thermodynamic factor G( γ, T ) with shear can be defined. A first clear observation from Fig. 7 is that G( γ, T ) must be typically of order unity. If we make the reasonable assumption that the barrier ∆E is of the same order of magnitude as in usual supercooled liquids (typically of order ǫ, where ǫ is the energy scale of the interaction), it appears that the usual activation at the temperatures under study (T < 0.06ǫ/k B ) is unable to account for our observations. Moreover, any effective temperature T eff , which would be used to explain the observed values of G( γ, T ), should depend on shear rate only weakly. The shear rate dependence of G( γ, T ) is presented in Fig. 8. As the effective temperature can only increase with shear rate, the data clearly shows that the energy barrier must increase with shear rate. This is consistent with the findings of Ref. 31 . However, we found that the parabolic fit (in the variable γ) of G( γ, T ), which was used in Ref. 31 to treat the data, is not sufficient to describe our results. Instead we use a cubic dependence of shear rate:

G( γ, T ) = G (0) + 3 n=1 B n (T )| γ| n , (6)
where G (0) characterizes the limit of the energy barrier hight at zero shear rate and B n (T ) are the correction factors. Equation ( 6) is an expansion in powers of the dimensionless strain rate γτ as a small parameter, γτ ∈ [0.00001, 0.01], for G( γ, T ) about the equilibrium value G (0) . Then, Eq. ( 5) can be extended to the following form:

τ -1 n = A| γ| exp - 3 n=1 B n (T )| γ| n , (7)
where A ∝ exp [-G (0) ] characterizes the zero shear rate limit of the thermodynamic factor. The fits of Eq. ( 7) are shown as continuous lines in Fig. 7 and the obtained numerical values of all parameters involved in Eq. ( 7) are given in Table I. Note that Eq. ( 7) should be considered as an attempt to extend the main equation of CNT for non-equilibrium case 31 of a glassy system under shear at very deep supercooling, where the inherent structural rearrangements are extremely slow and the external shear affects the transport properties directly as was discussed in Refs. 31,33 . The fact that the coefficients B n (T ) are, in fact, almost independent of the temperature shows that the process is activated by shear, rather than by thermal fluctuations.

In order to obtain an independent determination of the effective temperature, it would be necessary to determine an energy barrier E 0 at zero temperature and zero shear rate for the formation of a crystalline nucleus, and to obtain T eff as the ratio E 0 /G (0) . Unfortunately no method is known to obtain E 0 directly. A possibility, which we have not pursued in this work, would be to extrapolate results from high temperatures. As a result, our expression for the nucleation rate involving an effective temperature cannot be seen as a direct validation of this concept, but rather as a plausible extrapolation of other results 29,30 .

Our results demonstrate that the formation of a nanocrystalline phase in a driven, low temperature system proceeds via a nucleation mechanism, not dissimilar to the one observed in the standard liquids at moderate undercooling. Nucleation theory in its classical form is not applicable, in view of the very low temperatures. However, our results are consistent with a nucleation approach that would involve an effective temperature and a nucleation barrier, which increases with shear rate, as already reported in supercooled liquids.

This work could be extended in several directions, in particular:

(i) Metallic glasses often have quasi-crystalline phases. For the Dzugutov system studied here, the quasicrystalline phase can coexist with the liquid (or amorphous) and crystalline phases (see Fig. 4 in Ref. 34 ), and the spontaneous growth of quasi-crystals in this system was recently investigated by Glotzer 35 . It could be interesting to study the competition between quasi-crystal formation and crystal nucleation under an external drive.

(ii) Glassy materials are characterized by a waiting time dependence of their properties, or aging, and a state that depends on the cooling rate. The cooling process and the waiting time could have an influence on the ordering of the system under an external driving.

(iii) Although the one-component system considered in this study is a good glass-former, bulk metallic glasses are, as a rule, multicomponent (five different kinds of atoms, at least). Therefore, it would be of interest to study the influence of structural inhomogeneity appeared from polydispersity on the induced structural ordering.

Figure 1 :

 1 Figure 1: Time-dependence of different characteristics computed for a single sheared sample with the strain rate γ = 0.005τ -1 at the temperature T = 0.03ǫ/kB : (a) Global orientational bond-order parameter Q6(t). Inset: radial distribution of particles formative the ordered cluster structures. It is defined for the moment indicated by arrow. Pronounced peaks of radial distribution appear at the distances ∼ σn 1/2 , n = 1, 2, 3, . . . is the number of a peak, that is typical for FCC structures. (b) Number of solidlike particles s (1) (t). Inset shows the same in logarithmic scale. (c) Number of clusters s (0) (t).

  of τ) 〈ζ〉 (in units of σ)

Figure 3 :

 3 Figure 3: (Color online) Density profiles of solidlike particles for the instantaneous configurations at the different times after shear initiation. The results are typical for all studied samples. Profiles presented here are for the same system as in the previous figure.

Figure 4 :

 4 Figure 4: (Color online) Largest cluster at the different time steps prior to transition. N indicates the size. The sample is the same as for the previous figures.

Figure 5 :

 5 Figure 5: (Color online) Main: Size of the largest cluster versus time. Broken line corresponds to the value 160. The study sample is still the same as for the previous figures. Inset: Three largest clusters at the moment indicated by arrow. The sizes of the clusters are N1 = 177, N2 = 130 and N3 = 126 particles.

1 Figure 6 :

 16 Figure6: (Color online) (a) Orientational bond-order parameter of a sheared system (amorphous initially) versus strain. Each curve corresponds to a fixed strain rate γ ∈ [0.0001, 0.01]τ -1 and is averaged over results obtained for at least 50 independent samples, which were sheared under the same conditions at a temperature T = 0.03ǫ/kB . (b) Main: Orientational bond-order parameter versus rescaled transition strains ξ = γ/γn. Inset: Mean first passage time distributions multiplied by the corresponding strain rates. The solid lines are results of simulations, the dotted lines present the fits to Eq. (3).

Figure 7 :

 7 Figure 7: (Color online) Transition rate versus strain rate at three different temperatures. The points with error bars are the data extracted from our molecular dynamics simulations.The solid curves are the fits using Eqs. (5) and (6). The short straight line corresponds to a unit slope.

Figure 8 :

 8 Figure 8: (Color online) Change of the exponent of thermodynamic factor G( γ, T with the strain rate in regard to the case with γ = 0.0001τ -1 .

Table I :

 I Numerical data for different temperatures T on the parameter A and the correction factors B1, B2 and B3 to the equilibrium barrier height G (0) as obtained from the fitted simulation data.

	T (ǫ/kB)	A	B1 (τ )	B2 (τ 2 )	B3 (τ 3 )
	0.01	2.44	430	-230 × 10 2	500 × 10 3
	0.03	2.78	417	-229 × 10 2	505 × 10 3
	0.06	2.86	390	-210 × 10 2	510 × 10 3

We acknowledge useful discussions with A. Tanguy, T. Biben, R. Khusnutdinoff and R. Yulmetyev and an interesting correspondence with J. Wedekind and F. Varnik. This work was supported by ANR project "SLLOCDYN" and partially by RFBR (Grant No. 08-02-00123-a).