Badreddine

Saidi

Laurence Giraud-Moreau

Abel Cherouat

Optimization of the single point incremental forming process for titanium sheets by using response surface

The single point incremental forming process is well-known to be perfectly suited for prototyping and small series. One of its fields of applicability is the medicine area for the forming of titanium prostheses or titanium medical implants. However this process is not yet very industrialized, mainly due its geometrical inaccuracy, its not homogeneous thickness distribution... Moreover considerable forces can occur. They must be controlled in order to preserve the tooling. In this paper, a numerical approach is proposed in order to minimize the maximum force achieved during the incremental forming of titanium sheets and to maximize the minimal thickness. A surface response methodology is used to find the optimal values of two input parameters of the process, the punch diameter and the vertical step size of the tool path.

Introduction

The single point incremental forming process (ISF) is an emerging process with a high industrial interest. This process is based on the use of a spherical tool, which is moved along CNC controlled tool path (figure 1). During the incremental forming process, the sheet blank is fixed in sheet holder. The tool follows a certain tool path and progressively deforms the sheet. Because incremental forming is a dieless process, it is perfectly suited for prototyping and small volume production [1,2].

Medical implants are one of the potential fields of application of the incremental forming process, due to the need of product customization to each patient. One of the major advantages of the ISF process is the ability to program several punch trajectories on the same machine in order to obtain different shapes perfectly suited to each patient. Titanium alloys are generally used by the biomedical industries because of their corrosive resistance and great biocompatibility with the human body. Titanium is recognized by its high strength-toweight ratio. But its high cost and difficult formability compared with aluminum make that the incremental forming of titanium sheet is little studied. Before the incremental forming process of titanium sheets is really industrialized, studies are necessary. The ISF process indeed suffers from a big slowness, geometrical inaccuracy, a not homogeneous thickness distribution which reduces its industrial suitability [1,3,4,5].

During this process, considerable forces can occur [5]. The control of the forming force is particularly important when titanium sheets are used to ensure the safe use of a CNC milling machine and preserve the tooling and machinery [5,6].

The aim of this paper is to determine input parameters of the process in order to minimize the maximum force achieved during the incremental forming of titanium sheets and to optimize the thickness reduction. Several studies were developed to optimize the incremental forming of aluminum sheets [2,8]. A little less was done for titanium sheets.

In this paper, a numerical procedure is developed in order to study the influence of two input parameters, the punch diameter and the vertical step on the forming force and the minimal thickness. To conduct this study, an experimental design combined to response surface methodology is used [2]. The experimental design is based on finite element analysis realized with ABAQUS/Explicit. Mathematical relations are defined from the response surfaces to predict the forming force and the thinning rate as function of the two input parameters and to define the input parameters allowing to maximize the forming force and maximize the minimal thickness.

Incremental forming of titanium sheets

The single point incremental sheet forming (SPIF) process is illustrated on the figure 1. Different trajectories along which the tool moves have been studied in the literature. The tool path most frequently used and that the authors of the paper have already studied during experimental works [3] is composed of rotational movements in the horizontal plane combined with small steps in the vertical direction (see figure 2). This tool path is used in this study. This paper concerns the incremental forming of titanium ASTM grade 2 sheets with a thickness of 0.5 mm. A truncated cone with a wall angle of 60°, a maximum diameter of 100 mm and a depth of 20 mm was chosen for this preliminary study about the optimization of the incremental forming process with titanium sheets. The sheet has a size of 300 mm x 300 mm.

Numerical simulation

During the incremental forming process, the sheet is held firmly by the clamping system so that there is no slipping. A full model of blank holder modeled by rigid surfaces has been chosen to represent the experimental clamping system. A vertical pressure of 4.3MPA [9] has been applied. Figure 3 presents the numerical model at the beginning of the process. The initial and final mesh of the sheet for a punch diameter of 10mm and a step size of 1 mm is shown on figure 4. In what concerns material modeling, an elasto-plastic model was chosen for the simulation. The titanium behavior was modeled by means of a Swift type hardening law :

V = k (H0+Hp) n (1)
A tensile test allowed generating the strain-stress curve for the titanium sheet. Figure 5 presents this stressstrain curve. An inverse analysis using the data from this tensile test has been used to identify the mechanical parameters. Table 1 summarizes the mechanical properties of the material obtained thanks to the inverse method.

Experimental design

The goal of this present study is to determine the values of the punch diameter (dp), and the vertical step size (Δz) which allow to minimize the forming force and maximize the minimal thickness. As the number of input variables is low, a full factorial design at three levels has been chosen to realize the experimental design. Table 2 presents the forming factors and their respected levels.

The numerical simulation of the incremental forming process has been performed with the ABAQUS explicit solver. The finite elements chosen to mesh the sheet are shell elements S4R with 4 nodes and reduced integration. The spherical tool is modeled by a rigid surface.

Optimization results

A quadratic polynomial equation has been considered in order to fit the numerical results and get a relationship between the response values and the input parameters.

The two polynomial models (2) and (3) of second degree considered for the forming force and minimal thickness are the following [10]:

Fmax = a0 + a1dp + a2 'z + a3 dp 2 + a4 dp 'z + a5 'z 2 (2) STH = b0 + b1dp + b2 'z + b3 dp 2 + b4 dp 'z + b5 'z 2 (3)
Table 4 presents the values of the coefficients, which allow to predict the forming force and the thinning rate according to the two input parameters. These numerical results have been obtained with ANOVA. The 3D response surfaces showing the variations of the forming force and the minimal thickness with respect to the punch diameter and the step size are given in figure 6 and7. These surfaces are useful to understand the interaction properties between both the punch diameter and step size and the output parameters. We can see on these figures that the increase of the tool diameter and the step size induces an increase of force forming. As regards the minimal thickness, a nonlinear response is observed. The maximal thickness for the considered input parameters is obtained for the minimal tool diameter and minimal step size.

These response surfaces have been used to determine the optimal values of the tool diameter and the vertical step size which allow to minimize the forming force and maximize the minimal thickness. The results are presented in the table 5.

STH 0,14-0,156 0,156-0,172 0,172-0,188 0,188-0,204 0,204-0,22 0,22-0,236 0,236-0,252 0,252-0,268 0,268-0,284 0,284-0,3 0,3-0,316 A quasi linear effect of the punch diameter and the step size on the forming force is observed. Figure 9 shows a quadratic effect of the input parameters on the thickness response. Figures 10 and11 presenting the standardized PARETO chart show that the punch diameter and the vertical step size have a similar effect on the forming force but the vertical step size has a more important effect than the punch diameter on the thickness. Figures 12 and 13 present the normal probability plot. The process parameters are considered significance to the experiment as there is a large angled slope of the straight line. We can see that all the points are relatively close to the straight line and thus that the data follow a normal distribution [8].

Conclusion

In this paper, a numerical approach based on finite element analysis and response surface methodology has been applied to optimize the incremental forming process of titanium sheets. A second-order quadratic function has been defined in order to predict the maximal forming force and minimal thickness as function of the punch diameter and the vertical step size.

The impact of the punch diameter and the vertical step size on the maximal forming force and minimal thickness has been studied.

The optimal values of the input parameters allowing minimizing the forming force and maximizing the minimal thickness have been defined from the response surfaces.

This study about the optimization of the incremental forming process with titanium sheets is a preliminary study. Only two input parameters were considered and the material low used in the finite element model doesn't take into account the general anisotropy. This study will be pursued in future works.

Figure 1 .

 1 Figure 1. Incremental forming process.

Figure 2 .

 2 Figure 2. Tool path.

Figure 3 .

 3 Figure 3. Numerical model.

Figure 4 .

 4 Figure 4. Initial and final mesh of the sheet.

Figure 5 .

 5 Figure 5. Stress-strain curve.

Figure 6 .

 6 Figure 6. Response surface associated to the forming force.

Figure 7 .Table 5 .

 75 Figure 7. Response surface associated to the minimal thickness.Table 5. Coefficients of the fit functions. Factor Low High Optimum dp (mm) 10,0 15,0 10,0 Δz (mm) 0,5

DOI: 10 . 68010011 Figure 8 .

 10680100118 Figure 8. Effects of the input parameters on the value of the forming force.

Figure 9 .

 9 Figure 9. Effects of the input parameters on the value of minimal thickness.

Figure 10 .

 10 Figure 10. Standardized PARETO chart for force forming.

Figure 11 .

 11 Figure 11. Standardized PARETO chart for the thickness.

Figure 12 .

 12 Figure 12. Normal probability distribution for force forming.

Figure 13 .

 13 Figure 13. Normal probability distribution for thickness.

Table 1 .

 1 Mechanical properties of the Titanium grade 2.

	Modulus of elasticity : E	100 Gpa
	Yield strength:	210 Mpa
	Hardening coefficient: n	0.084
	Coefficient of proportionality: K	697,7 MPa
	Strain : ε0	0.00157

Table 2 .

 2 Forming factors and levels.Nine numerical simulations of the forming process were carried out. The responses values are the minimal thinning and maximal forming force. Table3presents the numerical results performed according to the full factorial design.

	Factors	Factor	Level	Level	Level
		notation	1	2	3
	dp : Tool	A	10	12.5	15
	diameter (mm)
	Δz : Step size	B	0.5	0.75	1
	(mm)

Table 3 .

 3 Numerical results.

	n 0	Dp	Δz	Minimal	Maximal
		(mm)	(mm)	thickness	forming force
				(mm)	(N)
	1	10	0.5	0,255909	870
	2	10	0.75	0,222006	950
	3	10	1	0,172874	1030
	4	12.5	0.5	0,258584	920
	5	12.5	0.75	0,106811	980
	6	12.5	1	0,197069	1110
	7	15	0.5	0,205474	1060
	8	15	0.75	0,205505	1120
	9	15	1	0,274236	1180

Table 4 .

 4 Coefficients of the fit functions.

	1082,2 -80,7	233,3	5,1	-16,0 186,7
	1,48	-0,12 -1,40	0,0037	0,03	0,6

(2016) © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/).

(2016)


Samson Acoca-Pidolle
email: samson.acoca-pidolle@cefe.cnrs.fr

Perrine Gauthier

Louis Devresse

Virginie Pons

Pierre-Olivier Cheptou

Ongoing convergent evolution of a selfing syndrome threatens plantpollinator interactions

Keywords: POC, SAP, methodology: SAP, POC, lab

Plant-pollinator interactions has evolved early in angiosperm species radiation and is considered as a major trophic link in terrestrial ecosystems. Recent environmental changes however have led to pollinator decline that are susceptible to cause a disruption of plant-pollinator interactions. We used the resurrection ecology methodology to assay ancestral and contemporary descendant phenotypes in natural populations. We assessed patterns of mating system evolution in Viola arvensis, over the last three decades in the Parisian regional basin (France), a depauperate pollinator environment. Population genetic analyses revealed about 15% increase in realized selfing rates in the field during the last three decades. Using a two-generation common garden experiment, we documented trait evolution towards small and inconspicuous corollas, reduced nectar production and reduced attractiveness. Trait shifts were found to be convergent over the four studied populations. Our study demonstrates that plant mating system can evolve rapidly in natural populations in the face of recent environmental changes. Such rapid evolution towards a selfing syndrome illustrates the ongoing disruption of a major trophic interaction of natural ecosystems.

Main Text

Introduction

In the last 30 years, insect populations have declined by more than 75% [START_REF] Hallmann | More than 75 percent decline over 27 years in total flying insect biomass in protected areas[END_REF]. More specifically, we are facing a global decline of pollinators, like bees (Anthophila) [START_REF] Goulson | Bee declines driven by combined stress from parasites, pesticides, and lack of flowers[END_REF][START_REF] Grab | Agriculturally dominated landscapes reduce bee phylogenetic diversity and pollination services[END_REF]. The astonishing diversification of angiosperms, which started 100 billion years ago, is commonly thought to have been stimulated and hastened by the diversity of plant-pollinator interactions [START_REF] Whittall | Pollinator shifts drive increasingly long nectar spurs in columbine flowers[END_REF]. More than 80% of extant angiosperms rely on animals for pollination [START_REF] Ollerton | How many flowering plants are pollinated by animals?[END_REF], thus plants and pollinators have a close relationship. The decline of pollinators is susceptible to affect plant reproduction through the direct lack of pollination, and parallel declines in animal-pollinated plants and insects have been documented, over the last 50 years [START_REF] Biesmeijer | Parallel declines in pollinators and insect-pollinated plants in Britain and the Netherlands[END_REF]. For individual plant species, pollen limitation has been shown to represent a major selective pressure in the evolution of self-fertilization (selfing) in theoretical evolutionary models [START_REF] Lloyd | Self-and cross-fertilization in plants. II. The selection of self-fertilization[END_REF][START_REF] Thomann | Flowering plants under global pollinator decline[END_REF]. Indeed, hermaphroditic angiosperms predominantly mate via outcrossing i.e. with another individual, although recurrent evolution towards selfing is a common macroevolutionary shift [START_REF] Stebbins | Self fertilization and population variability in the higher plants[END_REF][START_REF] Igic | Ancient polymorphism reveals unidirectional breeding system shifts[END_REF]. In recent years, such a shift towards selfing has been experimentally reproduced in the absence of pollinators [START_REF] Bodbyl Roels | Rapid evolution caused by pollinator loss in Mimulus guttatus: rapid evolution after pollinator loss[END_REF][START_REF] Gervasi | Real-time divergent evolution in plants driven by pollinators[END_REF][START_REF] Ramos | Rapid plant evolution driven by the interaction of pollination and herbivory[END_REF] and illustrated in natural populations impoverished in pollinators [START_REF] Brys | Effects of human-mediated pollinator impoverishment on floral traits and mating patterns in a short-lived herb: an experimental approach[END_REF]. This evolution of selfing is commonly associated with changes in reproductive traits such as a decrease in the distance between male and female organs (herkogamy) [START_REF] Bodbyl Roels | Rapid evolution caused by pollinator loss in Mimulus guttatus: rapid evolution after pollinator loss[END_REF][START_REF] Ramos | Rapid plant evolution driven by the interaction of pollination and herbivory[END_REF][START_REF] Brys | Effects of human-mediated pollinator impoverishment on floral traits and mating patterns in a short-lived herb: an experimental approach[END_REF] but also a decrease in floral traits mediating plant-pollinator interactions, such as a decrease in corolla and petal size, and decreased scent and nectar production [START_REF] Gervasi | Real-time divergent evolution in plants driven by pollinators[END_REF][START_REF] Ramos | Rapid plant evolution driven by the interaction of pollination and herbivory[END_REF][START_REF] Brys | Effects of human-mediated pollinator impoverishment on floral traits and mating patterns in a short-lived herb: an experimental approach[END_REF]. A shift in traits leading to less conspicuous and attractive flowers has been well described at the phylogenetic level [START_REF] Sicard | The selfing syndrome: a model for studying the genetic and evolutionary basis of morphological adaptation in plants[END_REF]. This group of trait changes characterizes the evolution of a selfing syndrome [START_REF] Sicard | The selfing syndrome: a model for studying the genetic and evolutionary basis of morphological adaptation in plants[END_REF][START_REF] Ornduff | Reproductive biology in relation to systematics[END_REF].

Faced with pollinator decline, natural populations are expected to evolve towards reduced plant-pollinator interactions and increased autonomous selfing [START_REF] Thomann | Flowering plants under global pollinator decline[END_REF]. Demonstrating trait evolution requires the distinction of the genetic and environmental components of phenotypes. The "resurrection ecology" approach, which involves growing dormant seeds collected several or a few decades previously alongside their descendants [START_REF] Franks | Rapid evolution of flowering time by an annual plant in response to a climate fluctuation[END_REF][START_REF] Weider | Evolutionary aspects of resurrection ecology: Progress, scope, and applications-An overview[END_REF], is a powerful methodology to detect such evolution [START_REF] Weider | Evolutionary aspects of resurrection ecology: Progress, scope, and applications-An overview[END_REF]. Indeed, this approach allows for the measurement of traits of ancestral and recent genotypes in a common environment, thus revealing the genetic components of trait variation.

Using this approach, we analyzed patterns of mating system trait evolution in the field pansy (Viola arvensis). To do so, we compared floral traits in ancestral seed population (from the 1990's and early 2000's) with the contemporary descendant population collected in 2021 from four localities in the Parisian regional basin (fig. S1 and table S1, [START_REF] Scoppola | Viola section Melanium (Violaceae) in Italy. New data on morphology of Viola tricolor-Group[END_REF]). The study species is an annual weed of intensive crops (mostly wheat and oilseed rape) that has a zygomorphic, showy flower and mixed selfing rates [START_REF] Scoppola | Viola section Melanium (Violaceae) in Italy. New data on morphology of Viola tricolor-Group[END_REF]. The length of the nectar spur and observations in natura are indicative of pollination by long-tongued bees. Populations were collected across a range of 100 kilometers in the study region. In general the study region has degraded pollination services compared with the rest of France [START_REF] Martin | New indices for rapid assessment of pollination services based on crop yield data: France as a case study[END_REF]. The status of pollinators or pollination services in such intensive agricultural areas has been reported as declining in several studies over our timescale [START_REF] Goulson | Bee declines driven by combined stress from parasites, pesticides, and lack of flowers[END_REF][START_REF] Grab | Agriculturally dominated landscapes reduce bee phylogenetic diversity and pollination services[END_REF][START_REF] Potts | Global pollinator declines: trends, impacts and drivers[END_REF][START_REF] Burkle | Plant-pollinator interactions over 120 years: loss of species, co-occurrence, and function[END_REF][START_REF] Raven | Agricultural intensification and climate change are rapidly decreasing insect biodiversity[END_REF][START_REF] Janousek | Recent and future declines of a historically widespread pollinator linked to climate, land cover, and pesticides[END_REF]. Parts of Belgium (some 100 to 200 km north of our study area) have experienced similar agricultural dynamics as our study region, and studies there indicate that 32.8% of bee species are threatened or extinct based on comparisons with their presence before and after 1970 [START_REF] Drossart | Belgian Red List of Bees[END_REF]. These studies also highlighted a more severe decline or higher vulnerability of longtongued bees [START_REF] Goulson | Bee declines driven by combined stress from parasites, pesticides, and lack of flowers[END_REF][START_REF] Biesmeijer | Parallel declines in pollinators and insect-pollinated plants in Britain and the Netherlands[END_REF][START_REF] Janousek | Recent and future declines of a historically widespread pollinator linked to climate, land cover, and pesticides[END_REF][START_REF] Drossart | Belgian Red List of Bees[END_REF].

In this study, we combined genetic and phenotypic approaches. Screening of thirty-two individuals per population with microsatellite markers on the first generation (F1) allowed us to assess population genetic metrics in natural conditions. These metrics provide information on allelic richness of populations and any changes in the mating systems as measured by the inbreeding coefficient, FIS [START_REF] Hartl | Principles of Population Genetics[END_REF]. The F1 progenies, obtained in common garden, enabled us to establish a half-sib family design for the second generation (F2), with similar maternal environments for the ancestral and descendant populations. Each F2-population is represented by twenty families composed of five half-siblings, for a total of 800 individuals. Phenotypic traits were measured in the second generation (F2) in 2022 in an air-conditioned greenhouse. We quantified the evolution of traits mediating plant-pollinator interactions and morphological traits not directly involved in plant-pollinator interactions. We also performed a plant pollinator preference experiment by monitoring bumblebee visitation in mixed artificial populations containing genotypes of the two different source generations. Together, these approaches document on the evolution of plant-pollinator interactions and evolution of plant mating system induced in recent environmental changes.

Results

Characterization of changes in the genetic structure of populations

Our localities exhibited no significant differences in allelic richness nor genetic diversity between ancestral and descendant populations based on fifteen microsatellites markers (Table 1). In the case of a sampling bias of our populations or a strong demographic bottleneck or genetic drift between the two samples, we predicted a decrease in allelic richness and genetic diversity. Our results rule out these scenarii. On the other hand, we observed a consistent increase in selfing rates in all localities between ancestral and descendant populations (Table 1), with an average increase of around 15%, which is consistent with results found in other locations [START_REF] Cheptou | Rapid evolution of selfing syndrome traits in Viola arvensis revealed by resurrection ecology[END_REF].

Characterization of changes of morphological and reproductive traits

Descriptive multivariate analyses performed on eight floral traits associated with the mating system (corolla length, labellum width, floral area, number of guides, anthesis duration, flowering date, volume of nectar and average floral display) and the ability of plant to set seeds in the absence of pollinators (i.e. autonomous self-pollination), corrected for differences between localities, showed differentiation between ancestral and descendant populations. We choose these traits because they are usually associated with shifts in mating systems [START_REF] Bodbyl Roels | Rapid evolution caused by pollinator loss in Mimulus guttatus: rapid evolution after pollinator loss[END_REF][START_REF] Gervasi | Real-time divergent evolution in plants driven by pollinators[END_REF][START_REF] Ramos | Rapid plant evolution driven by the interaction of pollination and herbivory[END_REF][START_REF] Sicard | The selfing syndrome: a model for studying the genetic and evolutionary basis of morphological adaptation in plants[END_REF]. The descendent populations clustered together below the first bisector while ancestral populations clustered above the first bisector (Fig. 1). This result highlights consistent directional changes of floral traits in all four localities. Linear models on the nineteen traits measured revealed consistent patterns (Table 2). We detected a significant difference between ancestral and descendant populations for all the floral traits, except sepal length (Table 2, table S2). Plants have evolved a 10% decrease in floral area (Fig. 2A) through a decrease in the labellum width and flower length (Table 2). Flowers have fewer nectar-guides (Fig. 2B), which are assumed to promote recognition of the path to floral rewards. Indeed, a correlation between the number of nectar-guides and the rate of visits inducing pollination has been reported elsewhere [START_REF] Hansen | Floral signposts: testing the significance of visual 'nectar guides' for pollinator behaviour and plant fitness[END_REF]. In contrast, none of the morphological traits supposedly unlinked to pollinator attraction (vegetative traits and sepal length) showed consistent changes (Fig. 2C, fig. S2, Table 2). Based on the resurrection approach, patterns of trait changes between ancestral and descendant populations provide a signature of mating system and associated trait evolution. The lack of changes in vegetative traits rules out the possibility that these changes are due to inbreeding depression, which is highly unlikely to only affect floral traits [START_REF] Winn | Analysis of inbreeding depression in mixed-mating plants provides evidence for selective interference and stable mixed mating[END_REF]. In addition, evolution in petal length and not in sepal (fig. S2) or leaf length (Fig. 2C), despite their homologous origins [START_REF] Breuninger | Control of tissue and organ growth in plants[END_REF], provides strong evidence for specific selection against costly traits that mediate interactions with pollinators.

Regarding reproductive traits, the results revealed a marginally significant increase over time in the ability to set seeds in the absence of pollinators (Fig. 2D) and a significant decrease in seed weight (Table 2). Seed weight represents maternal investment in individual progeny. This shift towards increased reproductive insurance is associated with reduced seed investment, as observed in previous studies [START_REF] De Jong | Hamilton's rule, imprinting and parentoffspring conflict over seed mass in partially selfing plants: Parent-offspring conflict over seed mass[END_REF].

Our study also revealed that locality and locality x age interaction are major causes of trait variations. The locality effects indicate spatial differentiation in measured traits and the locality x age interaction reveals a site dependent evolution in time. The latter interaction effects were different in magnitude, but not in sign, between ancestors and descendants for all the traits, except anthesis duration (Table 2). Evolution towards a selfing syndrome in the four localities is in the same direction but differs in magnitude. These effects indicate low gene flow between localities. The consistent evolutionary pattern (Fig. 1) highlights convergent shifts in fairly independent populations.

Characterization of changes of plant-pollinator interactions

Measurements of nectar volume, which is a fundamental rewarding trait, revealed a consistent decrease of 20% on average (Fig. 3A, Table 2). Because nectar production is costly [START_REF] Pyke | What does it cost a plant to produce floral nectar?[END_REF] and only dedicated to rewarding pollinators, the breakdown of plant-pollinator interactions is expected to be associated with such reduced nectar production. Nectar production is a trait that contributes to the maintenance of stable plant-pollinator interactions [START_REF] Brandenburg | Hawkmoth pollinators decrease seed set of a low-nectar Petunia axillaris line through reduced probing time[END_REF], and its decrease in our study species may loosen the interaction between the plant and its pollinators.

How this evolution impacts plant-pollinator interactions from the point of view of pollinators is addressed in our preference experiment during visits of the long-tongued bumblebees Bombus terrestris to ancestral and descendant plants mixture. Bumblebees visited more frequently ancestral than descendant plants i.e., they have a consistent preference of ancestral plants that is associated with their higher attractiveness and reward (Fig. 3B, Table 2). The evolution towards reduced attractiveness is consistent with reduced nectar production and the evolution of a selfing syndrome.

Discussion

Overall, floral display, floral area, floral attractiveness, and nectar production decreases, enhance selfing and increase of realized selfing rates are thus consistent with the ongoing evolution of a selfing syndrome [START_REF] Sicard | The selfing syndrome: a model for studying the genetic and evolutionary basis of morphological adaptation in plants[END_REF][START_REF] Ornduff | Reproductive biology in relation to systematics[END_REF]. Importantly, trait evolution is convergent across the four populations investigated. This is a unique demonstrate for such an evolutionary change over a short-time scale. While angiosperm-pollinator interactions have evolved over the long-term, our study thus shows that current environmental changes can drive rapid evolution towards a breakdown of such interactions.

Recurrent evolution towards selfing has been repeatedly documented, but to our knowledge, only over the long-term [START_REF] Stebbins | Self fertilization and population variability in the higher plants[END_REF][START_REF] Barrett | The evolution of plant sexual diversity[END_REF]. This evolutionary transition is classically considered to be irreversible (the "selfing as an evolutionary dead-end" hypothesis, 9,35) and higher extinction rates of such phyla have been reported [START_REF] Goldberg | Species selection maintains self-incompatibility[END_REF]. Evolution towards selfing could thus be driven by natural selection over the short-term but could impede long-term plant population survival [START_REF] Cheptou | Does the evolution of self-fertilization rescue populations or increase the risk of extinction?[END_REF]. The increase in selfing rate is indeed expected to have genomic impacts on individuals and populations through inbreeding (homozygosity), reduced effective population sizes and increased genetic drift. Consequently, selfers may have a reduced potential response to natural selection [START_REF] Glémin | Adaptation and maladaptation in selfing and outcrossing species: new mutations versus standing variation[END_REF][START_REF] Noël | Experimental evidence for the negative effects of self-fertilization on the adaptive potential of populations[END_REF]. While we acknowledge that short-term evolution of selfing in Viola arvensis is no-doubt facilitated by its initial mixed mating system (see for instance 40), our results illustrate how mixedmating species may evolve towards a higher reliance on self-pollination and reduced plantpollinator interactions in the face of global change, in particular pollinator decline. If a rapid transition towards selfing syndrome is shared among Angiosperms, many lineages may be in the situation of an extinction debt [START_REF] Goldberg | Species selection maintains self-incompatibility[END_REF].

Direct effect of environmental changes (habitat destruction and fragmentation, agricultural land pollution and species introduction) currently affect pollinator abundance and cause their decline [START_REF] Goulson | Bee declines driven by combined stress from parasites, pesticides, and lack of flowers[END_REF][START_REF] Potts | Global pollinator declines: trends, impacts and drivers[END_REF]. A decline in floral resources is also caused by these same changes and as a result becomes a driver of pollinator decline [START_REF] Goulson | Bee declines driven by combined stress from parasites, pesticides, and lack of flowers[END_REF]. Over time changes in plant community composition are known to be accompanied by changes in nectar resources [START_REF] Baude | Historical nectar assessment reveals the fall and rise of floral resources in Britain[END_REF]. In addition, as we illustrate here, nectar production can also rapidly evolve within species due to a relaxed selective pressure associated with pollinator decline. Decreases in nectar production may in fact reinforce the pollinator decline if levels fall below those necessary to sustain wild bee populations. The impact on pollinator populations is thus doubled as they become victims of both environmental changes and plant trait evolution [START_REF] Weinbach | Eco-evolutionary dynamics further weakens mutualistic interaction and coexistence under population decline[END_REF]. There is thus a potential eco-evolutionary positive feedback loop that furthers pollinator decline and reinforces plant evolution towards a selfing syndrome. Our results could in fact explain plant-pollinator network degradation (documented in a previous study 22) and how such degradation may generate cascading effects in trophic networks in general, and not just plant-pollinator interactions.

Pollinator decline is presumably ongoing and, even in the best scenarios, will continue in future [START_REF] Janousek | Recent and future declines of a historically widespread pollinator linked to climate, land cover, and pesticides[END_REF]. There is thus an urgent need to investigate the reversibility of plant evolution toward a selfing syndrome and assess whether or not the evolution of the plant-pollinator interactions has reached the threshold below which the breakdown of the interactions is inevitable [START_REF] Weinbach | Eco-evolutionary dynamics further weakens mutualistic interaction and coexistence under population decline[END_REF].

Materials and Methods

Study system and sampling

Viola arvensis is an annual species of wild pansy (Violaceae). It is commonly found in fields as a weed or in meadows. This species is described as self-compatible with a mixed mating system [START_REF] Scoppola | Viola section Melanium (Violaceae) in Italy. New data on morphology of Viola tricolor-Group[END_REF][START_REF] Cheptou | Rapid evolution of selfing syndrome traits in Viola arvensis revealed by resurrection ecology[END_REF]. Plants were sampled twice in four locations of the Parisian Basin (see table S1). The first sampling was done and conserved by the Conservatoire Botanique National de Bailleul for Crouy and Lhuys in the 1990's and by the Conservatoire Botanique National du Bassin Parisien for Commeny and Guernes in the 2000's. Seeds were collected by sampling fruits on a minimum of 100 mothers picked randomly in each population. The second sampling was performed in exactly the same locations in February 2021 by randomly sampling 40 or more seedlings. This sampling scheme produced four couples, one for each location, of ancestral and descendant populations.

Experimental design

To avoid maternal and potential storage effects, a refresher generation (F1) was grown in common garden for plants issued from the different populations. As we have seeds and seedlings, common conditions started at seedling stage. For ancestral populations, in mid-November 2020, seed lots were placed in germination chambers in petri dishes with vermiculite under light at 15°C for 12h and in the dark at 6°C for 12h. 32 seedlings of ancestral and descendant populations were transplanted in late February 2021 in open ground in the Baillarguet site of the experimental platform "Terrains d'Expériences" of the Labex Cemeb, 10 km north of Montpellier, France (43°40′53.92′′N, 3°52′28.74′′E). Populations were separated from each other by insect-proof mesh greenhouses (fig. S1). During full blooming, in mid-April, we introduced commercial bumblebee hives, Natupol from Koppert France® (Cavaillon, France), into the greenhouses to produce open-pollinated fruits that we collected. This refresher generation was used to generate seeds produced in the same environment for the experimental generation (F2) [START_REF] Franks | Using the resurrection approach to understand contemporary evolution in changing environments[END_REF]. This F1 generation allows us to generate a family design in the F2 generation. We randomly selected 20 mothers per population in the F1 generation to construct 20 families in the F2 generation. Each family is composed of five siblings (fig. S1). In December 2021, 30 seeds per family were placed in germination chambers in petri dishes with vermiculite under light at 15°C for 12h and in the dark at 6°C for 12h. Plants were placed in individual pots (11 × 11 × 11 cm 3 ) filled with standardized soil in late February and randomized in a greenhouse of the experimental platform "Terrains d'Expériences" of the Labex Cemeb, Montpellier, France. Watering was manual and done when needed. Three families did not germinate enough to obtain 5 plants and thus were composed of one, two or four siblings. In total, 792 plants were followed during this experiment. When all the morphological measurement were performed, plants were moved outside to finish their life cycle.

Population genetics F1 plants were genotyped using fifteen microsatellite markers, 8 design for Viola arvensis [START_REF] Cheptou | Rapid evolution of selfing syndrome traits in Viola arvensis revealed by resurrection ecology[END_REF] and 7 design on Viola tricolor [START_REF] Latron | Polymorphic nuclear markers for coastal plant species with dynamic geographic distributions, the rock samphire (Crithmum maritimum) and the vulnerable dune pansy (Viola tricolor subsp. curtisii)[END_REF], a closely related species [START_REF] Scoppola | Viola section Melanium (Violaceae) in Italy. New data on morphology of Viola tricolor-Group[END_REF]. DNA was extracted from fresh leaves using Dneasy Plant Mini Kit from QIAGEN®. Analyses of metrics such as allelic richness, expected heterozygosity (He), using the R package genepop [START_REF] Rousset | genepop: Population Genetic Data Analysis Using Genepop[END_REF], were performed to look at a change in genetic diversity over time. A drastic decrease in allelic richness (or He) is indicative of genetic drift, strong population bottlenecks or a sampling bias. Heterozygosity deficiency (FIS) analyses were performed, using the R package HierFstats [START_REF] Goudet | hierfstat: Estimation and Tests of Hierarchical F-Statistics[END_REF], on fourteen microsatellite markers. FIS was used to estimate the selfing rate using the following formula: .

Confidence intervals were obtained by bootstraps over loci of population (10,000). To perform an overall test on FIS, we pooled individuals of the same location, independently of their age and performed bootstraps to construct again two populations of 32 individuals again. We calculated the FIS of these two reconstruct populations and reproduce the procedure 10 000 times per locality on each locality. Then, we recorded the number of cases where the absolute value of the difference in reconstruct FIS is higher than the true difference in FIS measured. This number divided by 10 000 gives us a P-value for H0: "True difference = 0". We performed a Fisher's combined probability test on the P-values obtained for each locality [START_REF] Fisher | Statistical Methods for Research Workers[END_REF].

Measurements of plant traits

Plants traits were measured on the test generation (F2). All the morphometric measurements (fig. S3) were performed by a single experimenter during all the experimentation, to avoid any experimenter effect, using a digital caliper to 0.01 mm. Corolla length, labellum width, spur length, sepal length, number of guides and anthesis duration were measured on the five first opened flowers of the main shoot of each plant. We chose the first five opened flowers of the main shoot because more exhausted measurements in the F1 showed that they well capture the floral variance of an individual. Measurements were made day+2 after anthesis, or day+3 when falling on weekends. The opening of flowers was followed every working day to assess the anthesis duration. Floral area approximation is the product of corolla length and labellum width. Sepal/corolla lengths ratio is the ratio of sepal to corolla lengths. Around 4,000 flowers were followed. Flowering date is the date of anthesis of the first flower of a plant. The number of flowers opened per plant was counted each working day, from the opening of the first flower, for one month and half, and two more times with one week and then two weeks delay between counts. This number was integrated using trapezoidal rule and divided by the duration between the flowering date of the plant and the last day of measurement, giving the average floral display. Nectar volume was collected and cumulated from three flowers per plant using 0.5 µL microcapillaries. It was measured randomly on individuals of the two populations of the same location simultaneously to avoid differences in nectar production caused by the time of measurement. Thus, for this measure, Locality effect (see "Statistical analysis") is confounded with the day of measurement. It was measured on half of the plants from each population, representing all the families, thus 400 plants were sampled.

Vegetative traits such as rosette diameter, leaf length and leaf width were measured at the flowering date of the plant. Leaf length and leaf width were measured on the biggest leaf at this time, and the product of these two measures gives leaf area approximation. Biomass was measured on dried naturally dead plants at the end of their life cycle.

s Fis = 2F IS 1 + F IS 7
Reproductive traits were obtained while counting and weighing all the seeds contained by fruits to the nearest 0.1 mg using a precision balance. After all the measurements plants were placed by groups of twenty in insect-proof mesh greenhouses which contained commercial bumblebee hives, Natupol from Koppert France® (see Bumblebee preferences experiment below). After 2 hours minimum, plants were moved outside in open pollination context. Youngest receptive flower and oldest non-receptive ones on the same shoot were marked before the exposition to pollinator to differentiate fruits issued from flower in zero-pollinator context (selfpollination) versus fruits produced in open pollination. Thus, open pollination fruits are issued from flowers produced later in the shoot. For each plant, just before the dehiscence of the fruit produced by self-pollination, we sampled the two fruits at the same time and recorded the difference in the flower position number separating collected fruits (minimum 1, by definition, to maximum 5, with a mean of 1.6). The sampling of the two fruits was performed at the same time because it was too complicated to follow 1,600 fruits independently. Thus, open pollination fruits are not fully developed and comparing weights between open pollination and self-pollination fruits has a slight bias. However, comparing each measure between populations is fully possible while differences in floral positions are distributed randomly between populations. The seed weight reported is the measure of the weight of seeds contained by each fruit divided by the number of seeds counted in the fruit. The ratio of the number of seeds produced by self-pollination to the number of seeds produced by open pollination is a good proxy of the self-pollination ability and is comparable among plants if we consider that in open pollination treatment we have a full seed set and an incomplete seed set after self-pollination treatment.

Bumblebee preferences experiment

Plants were exposed to commercial bumblebee hives, Natupol from Koppert France®, in insectproof mesh greenhouses once all the trait measures had been made. In each insect-proof greenhouse, we placed, for the entire experiment, a bumblebee hive of around 30 workers with a queen, with nutrient reserves. We randomly disposed plants from the same location in a group of 20 plants, composed of 10 plants from each ancestral and descendant population randomly selected in one of the two insect-proof greenhouses. Pots were equally spaced at 35 cm. The observer, stationed in the entrance of the insect-proof greenhouses, selected a foraging worker bumblebee, and recorded each flower visit during 10 to 15 minutes. Greenhouse number (1, 2), position in the greenhouse (1-20) and plant identities were also recorded. Experiments were performed between the 2 nd and the 19 th of May 2022. We performed the experiment on 8 groups for each locality but just 6 and 7 observations were kept for Lhuys and Crouy respectively, because no pollinators foraged enough time in the discarded observations. On the 580 plants in our dataset, 3 are individuals exposed two times for experimental purpose. During a 10-15 minutes test, 4.28 plants were not visited in average and 124 in total for all the tests. In the analyses, we sorted out plants visited and unvisited (see below). We calculated the proportion of visits per plant as the ratio between the number of visits on a given plant divided by the total number of visits during a flight.

Statistical analysis

Changes in traits were evaluated using linear discriminant analysis (LDA) and univariate linear mixed models (LMMs). We performed the LDA over nine traits supposed to be linked with reproductive system; corolla length, labellum width, floral area, number of guides, anthesis duration, flowering date, volume of nectar, average floral display and ratio of seeds produced in self-pollination to open pollination. For traits with repeated measures, we used the means of the trait per individual. We normalized values by locality in order to focus on Age effect. We predefined only the population (8 groups) in the LDA and not the age. LDA gives linear combinations of traits (discriminant functions) that attempt to capture differences among groups.

In our case, it produced seven discriminant functions. The two first functions represent 55% and 22% of the proportion of the trace.

Univariate analyses were conducted using LMMs with each trait as response variable, Age (ancestral or descendant), Locality (4 localities), and the interaction between Age and Locality as fixed factors, and Family as random factor. In order to account for repeated measures on the same individual for floral measurement, we add Individual as random factor nested in Family. For seeds traits, we accounted for the difference in flower position between the two fruits collected adding this difference as random factor nested in family. When significant Age effect, analysis was followed by T-test comparison, on the means of the trait per individual for repeated measures, of ancestral and descendant populations of each locality. p-values were corrected using Bonferroni correction. Rosette diameter, leaf length, width and area, nectar volume and ratio of seeds without/with pollinators were log-transformed to approach normal distribution of residuals. In these analyses, Age significant effect represents a signal of global evolution signal of traits between ancestral and descendant populations, Locality significant effect, represents a difference in traits between populations and so population polymorphism spread in this species. A significant effect of the interaction between Age and Locality represents a difference in the magnitude or in the sign of the evolution in time between localities. In the case of a significant Age effect and a significant interaction effect, if evolution is consistent (same sign of evolution) among localities and only the magnitude change, we can hypothesize a directional selective force that selects a trait in a direction convergently in every locality. In other cases, a significant interaction between Age and Locality indicates a local selective pressure on the traits or, perhaps more conservatively, an evolutionary drift of the trait. These cases are beyond the scope in this study.

For bumblebee preference analyses, we first performed a GLMM to test the effect of Age on the number of unvisited plants (proportions: CoA: 0.24, Co21: 0.35; CrA: 0.26, Cr21: 0.26; GuA: 0.15, Gu21: 0.18; LhA: 0.13, Lh21: 0.12). We modelled it with Age and Locality as fixed effect, Family, Greenhouse, Position nested in Greenhouse and Date as random effects, with binomial distribution and logit link, as response variable was a Boolean variable with TRUE for visited plants and FALSE otherwise. There was no significant effect of Age in our dataset (N = 580 plants, estimates: AgeA-Age21 = 0.223, P = 0.296). Thus, we excluded unvisited plants (N = 124) in the rest of the analysis. We performed a LMM with log-proportion of visits per plant as response variable (N = 456), and the same structure as for other traits with Age, Locality, and the interaction between Age and Locality as fixed effect, Family as random effect and we added Greenhouse, Position nested in Greenhouse and Date as random effects. We then performed Ttest comparisons between Age for given Locality and corrected p-values using Bonferroni correction.

All our statistical analyses were performed using the R software (R Core Team, Version 4.2.1). LDA was performed using the "lda" function implemented in the package "MASS" [START_REF] Venables | Modern applied statistics with S[END_REF]. LMMs and GLMM were performed using the "lmer" and "glmer" functions implemented in the package "lme4" [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF]. 2). (B) The number of attractive guides showed an overall decrease in time (Table 2). (C) Rosette diameter, measured on each plant at the start of flowering, shows no consistent or significant differences between ancestors and descendants, like other vegetative traits (Table 2). (D) The log-ratio of seeds produced in selfpollination compared to open pollination is a proxy of selfing ability. We measured it by collecting one fruit in self-pollination and one in open pollination per plants (N = 693). This showed a trend towards an increase in selfing ability between ancestors and descendants (Table 2). The first two letters are the name of the locality (Co = Commeny; Cr = Crouy; Gu = Guernes; Lh = Lhuys). "A" (▲) ancestral population (collected in 2000 for Co, 1993 for Cr, 2001 for Gu and 1992 for Lh) and "21" (•) descendant population (all collected in 2021). There is an overall trend of a decrease in the volume of nectar over time (Table 1), which is significant for Cr and Gu (table S2). (B) Test of the preference of bumblebees. 10 plants of the ancestral and the descendant populations of a single locality were exposed together to bumblebees. We recorded the number of visits on each plant by a unique flying bumblebee for 10 to 15 minutes in 6 to 8 replicates per location and divided it by the total number of visits during the flight. Only visited plants are represented. There is an overall preference of ancestral plants (Table 1), significant for Co (table S2). The first two letters are the name of the locality (Co = Commeny; Cr = Crouy; Gu = Guernes; Lh = Lhuys). "A" (▲) ancestral population (collected in 2000 for Co, 1993 for Cr, 2001 for Gu and 1992 for Lh) and "21" (•) descendant population (all collected in 2021).

Figures and Tables

Table 1. Allelic richness, genetic diversity and selfing rates in our populations (mean ± SD or [CI]

). Average allelic richness and genetic diversity are based on fifteen microsatellites markers and selfing rate on fourteen. These markers were genotyped on the thirty-two individuals per population of the F1 generation issued from natural populations. For selfing rate, confidence intervals (CI) were obtained by bootstraps on loci (10,000). Fisher's combined probability test gave a significant overall increase in FIS ( = 35.04, P < 0.001).

The first two letters are the name of the locality (Co = Commeny; Cr = Crouy; Gu = Guernes; Lh = Lhuys). "A" ancestral population (collected in 2000 for Co, 1993 for Cr, 2001 for Gu and 1992 for Lh) and "21" descendant population (all collected in 2021).

1 .

 1 Figure 1. Multivariate differences between our population. Multivariate analysis (linear discriminant analysis performed on nine traits; N = 357 plants, almost equally distributed among groups, Nmin = 36 individuals; 7 points are not represented because out of scale) used to explore differences between ancestral and descendant populations. In the analysis, we predefined the eight different groups but not the status of the age factor. The first two functions represent 55% and 22% of the proportion of the trace respectively. The first two letters are the name of the locality (Co = Commeny; Cr = Crouy; Gu = Guernes; Lh = Lhuys). "A" (▲) ancestral population (collected in 2000 for Co, 1993 for Cr, 2001 for Gu and 1992 for Lh) and "21" (•) descendant population (all collected in 2021). Enlarged symbols are centroids of groups. The gray dashed line represents the first bisector.

Figure 2 .

 2 Figure 2. Evolutionary changes in floral, reproductive and vegetative traits. (A and B) Evolution of floral traits. We measured the first five developed flowers per individuals (N ≈ 4,000). (A) Floral area is the multiplication of labellum width and corolla length. It showed an overall decrease between ancestors and descendants (Table2). (B) The number of attractive guides showed an overall decrease in time (Table2). (C) Rosette diameter, measured on each plant at the start of flowering, shows no consistent or significant differences between ancestors and descendants, like other vegetative traits (Table2). (D) The log-ratio of seeds produced in selfpollination compared to open pollination is a proxy of selfing ability. We measured it by collecting one fruit in self-pollination and one in open pollination per plants (N = 693). This showed a trend towards an increase in selfing ability between ancestors and descendants (Table2). The first two letters are the name of the locality (Co = Commeny; Cr = Crouy; Gu = Guernes; Lh = Lhuys). "A" (▲) ancestral population (collected in 2000 for Co, 1993 for Cr, 2001 for Gu and 1992 for Lh) and "21" (•) descendant population (all collected in 2021).

Fig. 3 .

 3 Fig. 3. Decrease of plant-pollinator interactions. (A) Decrease of nectar production. We measured nectar volume as the sum of the volume in three flowers per plants on fifty plants per population (N = 400).There is an overall trend of a decrease in the volume of nectar over time (Table1), which is significant for Cr and Gu (tableS2). (B) Test of the preference of bumblebees. 10 plants of the ancestral and the descendant populations of a single locality were exposed together to bumblebees. We recorded the number of visits on each plant by a unique flying bumblebee for 10 to 15 minutes in 6 to 8 replicates per location and divided it by the total number of visits during the flight. Only visited plants are represented. There is an overall preference of ancestral plants (Table1), significant for Co (tableS2). The first two letters are the name of the locality (Co = Commeny; Cr = Crouy; Gu = Guernes; Lh = Lhuys). "A" (▲) ancestral population (collected in 2000 for Co, 1993 for Cr, 2001 for Gu and 1992 for Lh) and "21" (•) descendant population (all collected in 2021).

Table 2 . Linear models: effects of time, locality, and their interaction on the changes of floral, vegetative and reproductive traits

 2 . Values in parentheses are sample sizes. All traits are shown and the degree of significance of each effect on each trait is reported (LMM, *** P < 0.001, ** P < 0.01, * P < 0.05, • P < 0.1, n.s. P > 0.1, see tableS2for full statistical values). (+) a positive effect (increase), (-) a negative effect (decrease), (~) a difference of effect sign between localities (see tableS2for exact values).

			3.07 ± 1.67	2.73 ± 1.39	2.53 ± 1.60
		1.80 ± 0.86	3.20 ± 1.90	3.33 ± 1.68	2.00 ± 1.13
		0.274 ± 0.278	0.411 ± 0.265	0.311 ± 0.218	0.257 ± 0.269
	Genetic diversity (He)	0.095 ± 0.122	0.399 ± 0.272	0.322 ± 0.251	0.193 ± 0.211
	Selfing rate [95% CI] derived from heterozygosity	0.75 [0.69-0.81]	0.73 [0.64-0.83]	0.54 [0.40-0.71]	0.44 [0.29-0.55]
	deficiency, FIS	0.85 [0.81-0.94]	0.77 [0.70-0.84]	0.77 [0.74-0.84]	0.87 [0.77-0.94]

Acknowledgments

We thank Pauline Durbin and Fabien Lopez for help in taking care of the plants during the experiment, Marie-Pierre Dubois for lab-work guidance, the Conservatoire Botanique National de Bailleul and the Conservatoire Botanique National du Bassin Parisien for the supply of seeds and information on sampling locations, Antoine Deverge Merdrignac and Anaïs Vignerol for help in data collection, and John D. Thompson for improving this manuscript. Genotyping data was obtained using technical facilities of GenSeq of the "Institut des Sciences de l'Evolution de Montpellier" with the support of LabEx CeMEB an ANR "Investissements d'avenir" program (ANR-10-LABX-04-01).


Mohammed Guerboub

Steve Dave

Wansi Wendji

Carlo Massobrio

Assil Bouzid

Mauro Boero

Guido Ori

Evelyne Martin
email: evelyne.martin@unistra.fr

Impact of the local atomic structure on the thermal conductivity of amorphous Ge2Sb2Te5

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Thermal properties are expected to be sensitive to the network topology and yet no clearcut information is available on how the thermal conductivity of amorphous systems is affected by details of the atomic structure. To address this issue we use as a target system a phase change amorphous material (GST, i.e. Ge 2 Sb 2 Te 5 ) simulated by first-principles molecular dynamics (FPMD) combined with the approach-to-equilibrium molecular dynamics technique (AEMD) to access the thermal conductivity. Within density-functional theory (DFT), we employed two models sharing the same exchange-correlation functional but differing in the pseudopopotential (PP) implementation (namely, Trouiller-Martins, TM and Goedecker, Teter and Hutter, GTH). They are both compatible with experimental data and yet the TM PP construction results in a Ge tetrahedral environment largely predominant over the octahedral one, while the proportion of tetrahedra is considerably smaller when the GTH PP is used. We show that the difference in local structure between TM and GTH models impacts the vibrational density of states while the thermal conductivity does not feature any appreciable sensitivity to such details. This behavior is rationalized in terms of extended vibrational modes.

I. INTRODUCTION

Phonons act in crystals as heat carriers experiencing scattering events. However, at short scale, some phonons behave ballistically, their mean free path becoming larger than the dimensions of the system. This results in a reduction of the total thermal conductivity, as observed experimentally and theoretically [START_REF] Henry | Spectral Phonon Transport Properties of Silicon Based on Molecular Dynamics Simulations and Lattice Dynamics[END_REF][START_REF] Regner | Broadband phonon mean free path contributions to thermal conductivity measured using frequency domain thermoreflectance[END_REF][START_REF] Alvarez | Memory and nonlocal effects in heat transport: From diffusive to ballistic regimes[END_REF][START_REF] Palla | Interface thermal behavior in nanomaterials by thermal grating relaxation[END_REF] . The impact of downscaling on the thermal conductivity of amorphous materials and glasses is more elusive. The thermal conductivity is reduced at short scales in amorphous silicon, as shown by experiments [START_REF] Wingert | Thermal transport in amorphous materials: A review[END_REF] and modeling [START_REF] He | Heat transport in amorphous silicon: Interplay between morphology and disorder[END_REF] . In addition, first-principles molecular dynamics (FPMD) simulations of thermal transients (approach-to-equilibrium molecular dynamics, AEMD hereafter) have evidenced a reduction of thermal conductivity in other disordered materials (amorphous GeTe 4 7 , Ge 2 Sb 2 Te 5 [START_REF] Duong | First-principles thermal transport in amorphous Ge 2 Sb 2 Te 5 at the nanoscale[END_REF] and SiO 2 [START_REF] Martin | Thermal conductivity of amorphous SiO 2 by first-principles molecular dynamics[END_REF] ). A view commonly accepted for the heat transport in amorphous materials includes propagative heat carriers, the propagons (equivalent to phonons in the crystals) along with heat carriers that do not have well defined wave-vectors but are extended, the diffusons, and localized vibrations, the locons [START_REF] Allen | Diffusons, locons and propagons: Character of atomie yibrations in amorphous Si[END_REF] . The reduction of thermal conductivity mentioned above highlights the impact of extended modes in the heat transport of disordered materials.

Yet, it is still unclear if the thermal conductivity is sensitive to the local atomic environment, and if so, by which amount it is modified. In principle, this can be achieved by realizing different model systems of the same material, the difference being nature and connectivity of their building motifs.

While this strategy is somewhat cumbersome to implement experimentally, it falls entirely within the capabilities of quantitative atomic-scale modelling. This is exactly the guideline pursued in this paper, in which we focus on amorphous Ge 2 Sb 2 Te 5 (aGST), a phase-change material used in resistive memory devices for in-memory computing [START_REF] Joshi | Accurate deep neural network inference using compu-tational phase-change memory[END_REF] or as reconfigurable light management solution in nanophotonics [START_REF] Berguiga | Ultimate phase sensitivity in surface plasmon resonance sensors by tuning critical coupling with phase change materials[END_REF] . In a previous work [START_REF] Bouzid | Atomic-scale structure of the glassy Ge 2 Sb 2 Te 5 phase change material: A quantitative assessment via first-principles molecular dynamics[END_REF] , we have studied the structure of bulk aGST via FPMD by implementing a thermal cycle on a periodic box of 144 atoms. Several choices of the pseudopotential (PP) constructions within the PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] and BLYP [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF] exchange-correlation functionals were adopted. The combination of the BLYP exchange and correlation functional with the Troullier-Martins (TM) [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF] PP resulted in an excellent agreement between the calculated and measured total neutron structure factor, leading to a predominance of Ge atoms in tetrahedral en-

vironnment. An entirely different structure, yet still compatible with experiments in terms of total neutron structure factor, was obtained by employing BLYP along with GTH [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF] PPs, the proportion of Ge in tetrahedra decreasing from 69% (TM) to 22% (GTH). As a first purpose of this work, we intended to rule out size effects and confirm (or disprove) the difference of local structure be- into the relationship between a set of macroscopic properties and the network structure of a PCM material. We are able to exploit two available atomic structures for aGST, each one corresponding to a network organization in line with contrasting experimental descriptions (defective octahedra vs tetrehedra), to assess their impact on a fundamental macroscopic thermal property.

The paper is organized as follows. As a necessary prerequisite to the understanding of the issues inherent in the structural determination aGST, we briefly recall in Sec. II the context of the long lasting debate on the relative importance of tetrahedral vs octahedral basic units. In Section III, we provide a description of our calculation methodology and we introduce the FPMD models.

Section IV is devoted to the analysis of their structure. In Section V, an analogous comparative analysis is performed for the vibrational properties. Conclusions are contained in Section VI.

II. THE ATOMIC STRUCTURE OF AMORPHOUS GST: TETRAHEDRAL VS

OCTAHEDRAL BONDING GEOMETRIES

Having established that our investigation is fully rooted on the existence of different kinds of networks for aGST, it is worth contextualizing it within the longstanding debate on the relative importance of tetrahedral vs octahedral structural motifs. From the standpoint of modelling, we shall refer almost exclusively to structural descriptions obtained in the framework of FPMD, either employed as a self-contained tool of quantitative prediction or in combination with experimentally guided structural refinements as reverse Monte-Carlo (RMC). Early contributions addressing the presence of both tetrahedral (TE) and defective octahedral (DO) structural motifs were made available in Refs. 19-27. Vibrational properties were also considered in connection with the corresponding crystalline behavior. [START_REF] Akola | Density functional study of amorphous, liquid and crystalline Ge 2 Sb 2 Te 5 : homopolar bonds and/or ab alternation?[END_REF][START_REF] Caravati | First-principles study of crystalline and amorphous Ge 2 Sb 2 Te 5 and the effects of stoichiometric defects[END_REF][START_REF] Sosso | Raman spectra of cubic and amorphous Ge 2 Sb 2 Te 5 from first principles[END_REF] These findings have been complemented and enriched over the years by additional pieces of information focusing on the assessment of the relative percentage of TE and DO Ge sites. The predominance of octahedral sites over tetrahedral ones is a common feature of a number of independent FPMD calculations, [START_REF] Akola | Density functional study of amorphous, liquid and crystalline Ge 2 Sb 2 Te 5 : homopolar bonds and/or ab alternation?[END_REF][START_REF] Akola | Structural phase transitions on the nanoscale: The crucial pattern in the phase-change materials Ge 2 Sb 2 Te 5 and GeTe[END_REF][START_REF] Caravati | Coexistence of tetrahedraland octahedral-like sites in amorphous phase change materials[END_REF][START_REF] Caravati | First-principles study of crystalline and amorphous Ge 2 Sb 2 Te 5 and the effects of stoichiometric defects[END_REF][START_REF] Cai | Structural fingerprints of electronic change in the phasechange-material: Ge 2 Sb 2 Te 5[END_REF][START_REF] Hegedüs | Microscopic origin of the fast crystallization ability of Ge-Sb-Te phasechange memory materials[END_REF][START_REF] Akola | Experimentally constrained density-functional calculations of the amorphous structure of the prototypical phase-change material Ge 2 Sb 2 Te 5[END_REF][START_REF] Akola | Amorphous structures of Ge/Sb/Te alloys: Density functional simulations[END_REF][START_REF] Caravati | Influence of the exchange and correlation functional on the structure of amorphous[END_REF] for which only one third of Ge atoms are in tetrahedral geometries. Interestingly, even smaller TE Ge sites percentages are found when resorting to a Machine Learning model, [START_REF] Mocanu | Modeling the Phase-Change Memory Material, Ge 2 Sb 2 Te 5 , with a Machine-Learned Interatomic Potential[END_REF] built within the same FPMD framework of Refs. 24 and 31. These results are not entirely in agreement with experimental evidence pointing toward more extended tetrahedral order (EXAFS, x-ray absorption fine structure) and, for the nearest neighbor distances, a concomitant sizeable overestimate affecting FPMD values. [START_REF] Kolobov | Understanding the phase-change mechanism of rewritable optical media[END_REF][START_REF] Baker | Application of Bond Constraint Theory to the Switchable Optical Memory Material Ge 2 Sb 2 Te 5[END_REF][START_REF] Hosokawa | Existence of tetrahedral site symmetry about Ge atoms in a single-crystal film of Ge 2 Sb 2 Te 5 found by x-ray fluorescence holography[END_REF] Observed "wrong" Ge-Ge bonds are also more compatible with a tetrahedral arrangement. [START_REF] Jóvári | wrong bonds' in sputtered amorphous Ge 2 Sb 2 Te 5[END_REF] To solve this inconsistency by invoking the notion of polymorphism in phase-change materials, Akola et al. [START_REF] Akola | Polymorphism in phase-change materials: melt-quenched and as-deposited amorphous structures in Ge 2 Sb 2 Te 5 from density functional calculations[END_REF] proposed that the asdeposited (AD) models, closer in spirit to experimental realizations, could lead to a larger number of tetrahedral sites with respect to melt-quenched (MQ) ones, thereby promoting shorter FPMD distances. As touched upon in Sec. I, an alternative manner to recover shorter interatomic distances while extending the impact of tetrahedral configuration for MQ structure was presented in Ref. 13 (this reference contains also a detailed comparison with previous structural data obtained from experiments and FPMD calculations). In that work, a specific choice for the description of the nonlocal part in the Troullier-Martins 17 pseudopotential construction for Ge resulted in unprecedented agreement for the total neutron structure factor and the interatomic distances. Ge atoms are found in a predominant tetrahedral network coexisting with a non-negligible fraction of Ge in defective octahedra. However, it was argued that this atomic-scale picture is not consistent with Raman spectral features, tetrahedral coordination arising from Ge atoms being absent in experiments and first-principles calculations. [START_REF] Sosso | Raman spectra of cubic and amorphous Ge 2 Sb 2 Te 5 from first principles[END_REF] Overall, while an unambiguous structural determination for aGST remains elusive, it is worth mentioning that anomalous x-ray scattering [START_REF] Hosokawa | Key experimental information on intermediate-range atomic structures in amorphous Ge 2 Sb 2 Te 5 phase change material[END_REF] focussing on both short and intermediate range order have somewhat reconciled TE and DO views on the Ge coordination, both manifesting themselves in the experimental patterns with close percentages.

On the theoretical side, it appears that refinements in the DFT-based description could be instrumental to an improved structural assessment. With this purpose in mind, inclusion of dispersion forces has attracted considerable attention for both amorphous GeTe and GST. [START_REF] Micoulaut | Towards accurate models for amorphous GeTe: Crucial effect of dispersive van der Waals corrections on the structural properties involved in the phase-change mechanism[END_REF][START_REF] Cobelli | Metal-semiconductor transition in the supercooled liquid phase of the Ge 2 Sb 2 Te 5 and GeTe compounds[END_REF][START_REF] El Kheir | Density functional simulations of decomposition pathways of Ge-rich GeSbTe alloys for phase change memories[END_REF] Ref. 41 provides a very recent review of state-of-the-art ideas and open issues on the atomic structure of amorphous GeTe and GST.

Finally, given the above context, we reiterate that one of the strength of the present work is to take advantage of two atomic structure for aGST, both compatible with distinct experimental determinations and produced in a concomitant fashion, in order to ascertain their sensitivity to a specific macroscopic thermal property, i.e. the thermal conductivity.

III. METHODOLOGY OF CALCULATIONS

Calculations were performed within the Car-Parrinello 42 method as implemented in the developers version 4.3 of the CPMD code 43 . We selected for the exchange-correlation part of the Kohn-Sham total energy functional the exchange formula proposed by Becke [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF] and the correlation one of Lee, Yang and Parr 16 (BLYP). The valence-core interactions are described by either numerical norm-conserving TM [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF] or by analytical norm-conserving GTH [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF] PPs. For the choices inherent in the Kleinmann-Bylander expression of the local and non-local parts of the pseudopotentials we refer to the rationale given in Ref. 13. Valence electrons are represented by a plane-wave basis set compatible with periodic boundary conditions, with a cutoff of 30 Ry, and the sampling of the Brillouin zone is limited to the Γ point. A fictitious electron mass of 500 a.u. and a time step ∆t = 0.12 fs (5 a.u.) are chosen to ensure optimal numerical conservation of the constants of motion.

In all canonical NVT simulations, the ionic temperature was controlled with a Nosé-Hoover [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] thermostat chain, [START_REF] Martyna | Nosé-Hoover chains: The canonical ensemble via continuous dynamics[END_REF] whereas for the fictitious electronic kinetic energy we used a Blöchl-Parrinello thermostat [START_REF] Blöchl | Adiabaticity in first-principles molecular dynamics[END_REF] with a target kinetic energy depending on the ionic temperature. K (5 ps), T =300 K (8 ps), T =600 K (10 ps) and T =900 K (8 ps). At T = 900 K, the mean square displacements of the Ge, Sb and Te atoms reported in Fig. 1 lead to diffusion coefficients larger than 10 -5 cm 2 s -1 (Table I). These values ensure loss of memory of the initial configuration and effective randomization. Quench to room temperature occurred via a first step of 30 ps at T = 600

D Ge (cm 2 s -1 ) 2.8×10 -5 2.2×10 -5 D Sb (cm 2 s -1 ) 2.2×10 -5 1.8×10 -5
D Te (cm 2 s -1 ) 2.3×10 -5 1.2×10 -5

K followed by a second step of 35 ps at T = 300 K. Statistics were collected on this final part of the thermal cycle.

IV. STRUCTURAL PROPERTIES A. Partial pair correlation functions

Partial pair correlation functions g GeGe (r), g GeSb (r), g GeTe (r), g SbSb (r), g SbTe (r) and g TeTe (r)

for the four models considered (144-TM, 504-TM, 144-GTH, 504-GTH) are reported in Fig. 2, with the corresponding bond lengths and coordination numbers given in Table II.

Focussing on g GeGe (r), g GeSb (r) and g SbSb (r) one observes more regular profiles for 504-TM and 504-GTH compared to the smaller size, indicative of noise reduction. In the case of g GeGe (r)

the amplitude of the first peak, corresponding to homopolar Ge-Ge bonds, is significantly more intense for 504-TM than for 504-GTH, with a ratio between the relative heights increased in com- parison with 144-TM and 144-GTH. The position of this peak (2.46 Å for 504-TM vs 2.47 Å for 144-TM, 2.54 Å for 504-GTH vs 2.75 Å for the 144-GTH) confirms the better performances of TM when compared to GTH for this specific feature, while the increase in size has the effect of approaching the GTH value to the experimental data (2.47 Å). Turning to g GeSb (r), we can observe that the positions of the first peak (see also Table II) for TM and GTH models are closer in the larger boxes (2.72 Å for 504-TM vs 2.74 Å for 504-GTH) as well as their heights. While the TM approach is in better agreement with experiments, this is a further sign of increased reliability of GTH when compared to TM within the DFT-FPMD framework. In the case of g GeTe (r), the first peak is located at a larger distance for both 504-GTH and 144-GTH models, although a shift toward smaller values is present in the case of the larger size (2.74 Å for the 504-GTH vs 2.80 Å for the 144-GTH). TM values compare very well with experiments (see Table II). The partial pair correlation functions g SbSb (r), g SbTe (r) and g TeTe (r) obtained via 504-TM and 504-GTH are very n GeGe 0.6±0.2 0.42@3.1Å 0.29@3.1Å 0.63@2.9Å 0.26@3.1Å n SbSb 0.41@3.2Å 0.36@3.4Å 0.55@3.3Å 0.54@3.2Å n TeTe 0.18@3.1Å 0.20@3.2Å 0.27@3.1Å 0.28@3.2Å n GeSb 0.59@3.3Å 0.43@3.3Å 0.41@3.1Å 0.31@3.1Å n GeTe 3.3±0.5 3.02@3.2Å 4.27@3.5Å 2.75@3.1Å 3.27@3.3Å n SbTe 2.8±0.5 2.99@3.4Å 3.85@3.4Å 3.19@3.5Å 2.60@3.3Å close over the whole range of distances with the exception of the first minimum in the g SbSb (r), II T = 300 K T = 600 K T = 900 K n GeGe 0.26@3.1Å 0.32@3.1Å 0.52@3.1Å n GeSb 0.31@3.1Å 0.36@3.2Å 0.43@3.3Å

lower
Finally, the coordination numbers n SbSb and n TeTe obtained via TM and GTH do not change appreciably with the system size, while n GeSb remains higher for both (N= 144 and N= 504) TM models. These considerations reflect the shapes of the partial pair correlations functions and are globally in line with the better quality of TM models in comparison with experiments.

A global view on the distribution of the coordination numbers for Ge, Sb and Te is shown in % of the fourfold coordinated Ge do not form a bond with a Ge or a Sb atom, to be compared with ∼ 50 % for the GTH models. These findings are in agreement with previous evidence linking the occurrence of tetrahedral motifs and the presence of Ge-Ge and Ge-Sb connections. [START_REF] Caravati | Coexistence of tetrahedraland octahedral-like sites in amorphous phase change materials[END_REF][START_REF] Hegedüs | Microscopic origin of the fast crystallization ability of Ge-Sb-Te phasechange memory materials[END_REF][START_REF] Akola | Experimentally constrained density-functional calculations of the amorphous structure of the prototypical phase-change material Ge 2 Sb 2 Te 5[END_REF][START_REF] Akola | Amorphous structures of Ge/Sb/Te alloys: Density functional simulations[END_REF] We have seen that Ge atoms are mostly fourfold coordinated for both sizes regardless of the Impact of the local structure on the thermal conductivity of aGST choice for the pseudopotential scheme. The question arises on whether or not this feature means a similar network organization, as it could be tentatively concluded. For N = 144, the results of Ref. 13 showed that this was not the case since the proportion of Ge atoms in tetrahedral units (69 vs 22 % for TM and GTH, respectively) was a major difference in the Ge environment between the BLYP-TM and BLYP-GTH schemes. This was obtained by resorting to the integral of the order parameter 51 q = 1 -3 8 ∑ i>k ( 1 3 + cos θ i jk ) 2 between 0.8 and 1. The same analysis performed for N = 504 gives 69% of tetrahedral Ge units for 504-TM model and 37 % for 504-GTH one, confirming the stronger tetrahedral character of the network obtained by adopting the TM BLYP pseudopotential formulated in Ref. 13. These values can be compared to those obtained in Ref. 36 where as-deposited (AD) and melt-quenched (MQ) models were created to demonstrate that the former method is more effective in favoring tetrahedral bonding geometries than the latter. Values obtained in Ref. 36 for tetrahedral coordination are 58% for AD and 36% for MQ, both smaller than the corresponding 69% reported above for 504-TM within MQ modelling.

Ever since the first FPMD modelling efforts on aGST, rings structures have been considered

to play an important role in driving the crystallization process. [START_REF] Hegedüs | Microscopic origin of the fast crystallization ability of Ge-Sb-Te phasechange memory materials[END_REF] Our calculated distribution of primitive rings, [START_REF] Le Roux | Ring statistics analysis of topological networks: New approach and application to amorphous GeS 2 and SiO 2 systems[END_REF] reported in Fig. 4, is normalized on the total number of rings containing between three and ten atoms. Considering the fourfold rings, the ABAB configuration with A being a Ge or Sb atom, and B a Te one are more numerous (85%) in 504-GTH than in 504-TM (74%) ( see insets of Fig. 4). This is consistent with the well known tendency of tetrahedral configurations to favor Ge-Ge homopolar bonds, [START_REF] Caravati | Coexistence of tetrahedraland octahedral-like sites in amorphous phase change materials[END_REF][START_REF] Akola | Amorphous structures of Ge/Sb/Te alloys: Density functional simulations[END_REF] as confirmed by the data presented in Tab. IV.

To close this section, a final word is in order on the comparative analysis of TM and GTH systems for differents sizes (N= 504, N= 144). Based on the structure of the various GST networks considered here, it appears that TM and GTH schemes are closer when the system size is larger (N= 504 against N= 144) and yet an unambiguous difference in terms of overall structural organization does persist. to the observation that S T (k) is reproduced to an acceptable extent by both 504-TM and 504-GTH, with TM remaining preferable over GTH. This rationale is substantiated by the calculation of the Wright parameter 53 R χ :

R χ =    ∑ i S Exp T (k i ) -S T (k i ) 2 ∑ i S Exp T (k i ) 2    1/2 (1)
that is inversely proportional to the degree of agreement with the experimental values S Exp T (k). The values for the four systems are given in Table V. While the BLYP-TM recipe remains the best performing, for larger size the discrepancy affecting the GTH-based result decreases significantly. TM remains the best choice for reproducing experimental data on the structure of bulk amorphous GST, pointing to a predominant Ge-centered tetrahedral organization. The size of the supercell (504 vs 144 atoms) has a moderate effect on the comparative performances of the two schemes, TM vs GTH, in regard of their adequacy to describe quantitatively the structure of bulk aGST.

V. VIBRATIONAL PROPERTIES A. Density of states

A first insight into the impact of the GST network organization on the vibrational properties is given by the Fourier transform of the velocity-velocity autocorrelation function, providing the vibrational density of states (vDOS) (Fig. 6). The total vDOS differ at frequencies larger than 170 cm -1 . The contributions of Sb and Te show little sensitivity to the type of network generated by either TM or GTH PPs. The vDOS of Ge instead, features a shift at a higher frequency of the main peak for 504-TM (∼ 200 cm -1 against ∼ 150 cm -1 for 504-GTH) together with a range of values extending up to ≈ 275 cm -1 that explains the differences observed in the total vDOS. This result is in line with the observation that the local environment of Ge is the most sensitive to the modelling framework (TM or GTH), showing that higher frequencies in the vDOS of Ge appear when Ge is mostly in the TM tetrahedral environment.

We have employed the vibrational eigenvectors (obtained by diagonalization of the Hessian) to quantify the localization of the vibrational modes by calculating the inverse participation ratio (IPR) defined as: 54,55

1 p n (ω n ) = N N ∑ i=1 e i (ω n ) 2 2 / N ∑ i=1 e i (ω n ) 2 2 (2)
for a mode of eigenvalue ω n and eigenvector e i (ω n ). The IPR is equal to 1/N for a fully delocalized vibrational mode, while it reaches 1 at the extreme opposite (localization on a single atom). The IPR are reported for both models in Fig. 7. The IPR is significantly higher in 504-TM at high frequency, in the region where the vibrational states differ. These localized modes are not expected to contribute to the heat transport as it will appear clearly below when considering the results obtained for the thermal conductivities. The link between modes localization and Ge tetrahedral sites was already pointed out in Ref. 22.

B. Thermal conductivity

By calculating the thermal conductivity we aim at understanding the impact of the atomic structure on relevant quantities such as heat carriers (to be intended in their most general expression) not directly arising from averages over the time evolution of the atomic positions, as done up to now.

To this purpose, we resort to the approach-to-equilibrium molecular dynamics (AEMD) methodology, already applied to several amorphous systems. [START_REF] Duong | Thermal conductivity and transport modes in glassy GeTe 4 by first-principles molecular dynamics[END_REF][START_REF] Duong | First-principles thermal transport in amorphous Ge 2 Sb 2 Te 5 at the nanoscale[END_REF][START_REF] Martin | Thermal conductivity of amorphous SiO 2 by first-principles molecular dynamics[END_REF][START_REF] Bouzid | Thermal conductivity of glassy GeTe 4 by first-principles molecular dynamics[END_REF]57 AEMD is based on the implementation within molecular dynamics of a particular regime of the heat equation, where periodic and transient temperature profiles evolve in one dimension. The period is fixed at the largest possible value, i.e. the direction in the z direction of the elongated supercell (hereafter called L and equal to 40.6 Å). The transient decay time τ is related to the thermal conductivity κ by the expression :

κ = L 2 4π 2 ρC τ (3)
where ρ = 1/V is the number density or inverse of the volume, and C is the heat capacity. For this quantity, the deviation γ from the Dulong and Petit C DP = 3Nk B value (C = γC DP ) is calculated from the temperature derivative of the total energy averaged over the thermal cycle. For the TM scheme, a value of γ TM =1.365±0.015 was obtained in our previous work. [START_REF] Duong | First-principles thermal transport in amorphous Ge 2 Sb 2 Te 5 at the nanoscale[END_REF] The value obtained for the GTH scheme is γ GTH =1.130±0.018.

The transient and periodic temperature profiles are established by setting up an out-of-equilibrium initial configuration via the application of two local thermostats while the periodic boundary conditions are maintained (it is the so-called phase 1). The first thermostat, at T =400 K, acts on the atoms lying between 0 and L/2, while the second at T =200 K on those between L/2 and L. At the end of phase 1, lasting typically 5 ps, the temperature profile has a periodic form alternating between hot and cold regions, with a period equal to the length L in the z direction. The temperature profile in the range [0 -L], averaged during phase 1, is shown in Fig. 8.

During the transient that follows (phase 2), the temperature profile becomes sinusoidal (Fig. 8)

as predicted by the heat equation, and its amplitude decreases with time. The difference of temperature between the hot block (0 < z < L/2) and the cold one (L/2 < z < L) decreases following exp(-t/τ) as depicted in Fig. 9 on a semi-log plot. 1 st AEMD 6.1 0.12 5.0 0.12 2 nd AEMD 6.2 0.11 3.4 0.17

VI. CONCLUSION

The establishment of a correlation between atomic structure and thermal properties presented in this work has been done on the basis of two DFT-FPMD model of amorphous Ge 2 Sb 2 Te 5 , differing only in the type of PP used for the core-valence interaction. In Ref. 13 we demonstrated that the best agreement with neutron scattering data was obtained for the Trouiller-Martins (TM) [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF] PP, in conjunction with appropriate choices for the local and non local part of the Ge PP. We had also shown, and confirmed in the present work for a larger simulation box (N= 504), that the GST network topology issued by the TM model favors the predominance of Ge tetrahedrally coordinated, unlike the atomic structure obtained via the alternative GTH PP construction. [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF] The existence of two legitimate models for aGST has opened the way to a study of the sensitivity to the network atomic structure of the thermal conductivity κ. Recent FPMD-AEMD studies [START_REF] Duong | First-principles thermal transport in amorphous Ge 2 Sb 2 Te 5 at the nanoscale[END_REF] have demonstrated a sizeable reduction of κ in the nanoscale regime and a correct asymptotic behavior to the bulk value. As a first information on vibrational properties, we have calculated the corresponding density of states, exhibiting a more important tail at higher frequencies in the TM case. This change is not reflected by the values of the thermal conductivity that are found insensitive to the relative proportion of tetrahedral vs octahedral structural units. We can conclude that extended mode of vibrations, responsible for changes of thermal conductivity in the nanoscale range, are not affected by local variations in the atomic structure as those found when comparing TM and GTH models.

  tween TM and GTH by considering a larger supercell made of 504 atoms (112 Ge, 112 Sb and 280 Te). To this end, we have constructed a new GST-GTH model and extended the trajectory of the GST-TM already available in Ref. 8. Persisting structural differences are exploited to ascertain the sensitivity of the vibrational density of states and of the thermal conductivity to the atomic scale organization of aGST. Therefore, a second, major outcome of our work is a fundamental insight

6 Å 3 )

 63 Our models made of N= 504 atoms (termed concisely hereafter 504-TM and 504-GTH, meaning 504-TM model and 504-GTH model respectively) have an elongated shape (20.3×20.3×40.designed to study thermal transport as a function of the length along the z direction via the AEMD technique as explained in Section III. This non-cubic shape remains fully compatible with the objective of representing the structure of a bulk material. Periodic boundary conditions are applied in the three directions. Thermal cycles are employed to randomize the initial atomic configuration and to quench afterwards into the amorphous state. For the models containing 144 atoms, we implemented in Ref. 13 five trajectories of respectively 40 ps at T = 300 K, 20 ps at T = 600 K, 70 ps at T = 900 K, 50 ps at T = 600 K and 42 ps at T = 300 K. The 504-TM model is the same as model D2 of Ref. 8. However, the trajectory at T = 300 K has been extended by 35 ps to increase the statistical accuracy of the average properties. The initial configuration used for the 504-GTH is taken from the atomic positions of 504-TM relaxed within the BLYP-GTH computational setup. Then, a thermal ramp has been implemented as follows : T = 10 K (5 ps), T =100

FIG. 1 .

 1 FIG. 1. Mean square displacement of Ge (red line), Sb (blue line) and Te (green line) in 504-GTH at T = 900 K.

FIG. 2 .

 2 FIG. 2. Partial pair correlation functions of amorphous Ge 2 Sb 2 Te 5 . The color codes for the four models are as follows: 504-TM: red, 504-GTH: blue, 144-TM: orange and 144-GTH: green.

Fig. 3 .FIG. 3 .

 33 Fig. 3. The local environnement of Sb and Te is very similar within the TM and the GTH schemes at N=504, as an indication that the differences observed at N=144 decrease with increasing size. 504-TM and 504-GTH feature, as 144-TM, a dominant threefold coordination for Sb and twofold for Te. Ge atoms are mostly fourfold coordinated in all realizations. Worth of interest is the change in the nature of the second most frequent population for Ge. In 144-GTH there was a significant amount (34 %) of fivefold coordinated Ge atoms, while threefold atoms are instead found in a similar proportion (32 %) in 504-GTH. Another difference appears when looking at the chemical species connected to fourfold coordinated Ge atoms, shown in Table IV. Whereas Ge-Ge and Ge-Sb bonds are observed for all cases, for example taking the form of GeTe 3 or SbTe 3 structural

C.FIG. 4 .FIG. 5 .

 45 FIG. 4. Distribution of rings of size n containing up to n = 10 atoms. For each model, the percentage of fourfold rings found in ABAB (A=Ge or Sb and B=Te) configuration is given in the inset.

FIG. 6 .

 6 FIG.6. Total and partial vibrational densities of states for the models 504-TM (red lines) and 504-GTH (blue lines).

FIG. 7 .

 7 FIG.7. Inverse participation ratio of the vibrational states for the models 504-TM (red lines) and 504-GTH (blue lines). The total vibrational DOS are reported in light grey to recall the range of frequencies involved in each case.

  The transient time τ and thermal conductivity κ obtained in this way for 504-TM and 504-GTM models are reported in TableVI. The error bar results from the account of several AEMD runs via extensions of the phase 1 by 5 ps. This allow reaching a new initial point of the transient regime fully decorrelated from the previous one. Within the statistical error, the thermal conductivity does not differ when comparing the two models. This result exemplifies the intriguing behavior of the vibrational properties of amorphous materials, with local vibrational modes impacting the vDOS, while extended ones, involved in the heat transport, are insensitive to the local structural features. It is worth noting that these conclusions hold true for two systems very much similar in terms of comparison with the experimental total structure factor and yet quite different when focussing on the local network organization (as shown by the percentages of tetrahedral vs octahedral structural motifs).

FIG. 8 .

 8 FIG.8. Temperature profiles averaged during phase 1 (black curve) and phase 2 (blue curve) and fit by a sinusoidal variation in phase 2 (red curve). AEMD of the 504-GTH model.

FIG. 9 .

 9 FIG. 9. Temperature difference between the hot and cold blocks during phase 2 (blue line) and exponential fit to extract the transient time τ (here equal to 3.4 ps). Data refer to the 504-GTH model.

TABLE I

 I

. Diffusion coefficients at T = 900 K 504-GTH 504-TM

TABLE II .

 II Bond lengths (r αβ ) and partial coordination numbers (n αβ ) extracted from the partial pair correlation functions g αβ (r) of amorphous Ge 2 Sb 2 Te 5 . Experimental results are from Refs. 49 and 50. The results of the four models are reported: 504-TM, 504-GTH and the two models of Ref. 13, 144-TM and 144-GTH. Bond lengths are taken as the position of the first maxima of the pair correlation functions.The coordination numbers n αβ are obtained by integration of the pair correlation functions up to the first minima. n αβ are given in the table together with the position of the first minimum r min αβ at which they were determined (notation: n αβ @r min αβ ). For a comparison with previous FPMD data obtained in Refs. 20, 22, and 29 see Ref.13

		Exp.	144 atoms	504 atoms
	Ref. 49 Ref. 50 144-TM 144-GTH	504-TM 504-GTH
	r GeGe (Å)	2.47	2.47	2.75	2.46	2.54
	r SbSb (Å)		2.92	3.00	2.95	2.95
	r TeTe (Å)				2.95	2.89
	r GeSb (Å)		2.73	2.82	2.72	2.74
	r GeTe (Å) 2.61	2.63	2.66	2.80	2.63	2.74
	r SbTe (Å) 2.85	2.83	2.89	3.00	2.89	2.89

  in 504-TM. A slight overestimate is found when comparing Sb-Te interatomic distances (for The main effect of the larger size (504-GTH) in comparison with the smaller 144-GTH counterpart is a better agreement with available measurements for n GeTe and n SbTe (3.27 vs 3.3±0.5for n GeTe and 2.60 vs 2.8±0.5 for n SbTe ). For these coordination numbers, TM models have comparable accuracy. Concerning n GeGe , a very good agreement with experimental data is found for 504-TM (0.63 vs 0.6±0.2), while n GeGe values of 144-GTH and 504 -GTH are very close and differ substantially from experiments. TableIIIprovides the evolution of n GeGe and n GeSb during the ramp down of the thermal cycle for 504-GTH. The concentration of homopolar Ge-Ge bonds at T = 900 K reduces by a factor of two upon quenching, while a less pronounced decrease of the

both 504-TM and 504-GTH) to experimental values. We conclude that the partial pair correlation functions calculated with N= 504 within two different frameworks (TM and GTH) are less dissimilar than what they were found to be for N= 144 (see Ref.

13

). Yet, this confirms the better performances of the TM PPs as implemented in Ref.

13

for this class of systems. B. Coordination number and local environment number of Ge-Sb bonds is observed. Similar trends are obtained for the other models (504-TM, 144-TM and 144-GTH).

TABLE IV .

 IV Decomposition of the fourfold Ge environment in structural units with their chemical identity.

	For each structural unit, we give its relative weight (percentage) with respect to the total number of Ge
	atoms fourfold connected
	Fourfold-coord. Ge (%) 504-GTH 504-TM 144-GTH 144-TM
	Te 4	43.7	16.2	46.5	26.2
	GeTe 3	19.2	35.4	11.5	22.2
	SbTe 3	19.7	19.0	23.0	26.7
	GeSbTe 2	9.8	14.5	13.0	8.1
	Ge 2 Te 2		5.6
	Sb 2 Te 2	7.7	5.6	6.0	9.0
	Ge 2 SbTe				7.8
	Ge 3 Te		3.6

TABLE V .

 V Wright parameter R χ (%) for the four GST models.

	DFT recipe 144 atoms 504 atoms
	BLYP-TM	6.8	9.2
	BLYP-GTH	12.9	11.5

ACKNOWLEDGEMENTS

This work was funded by the French ANR via the project n. ANR-17-CE09-0039-02 "SIRENA". Calculations were performed by using resources from GENCI (Grand Equipement National de Calcul Intensif) (Grants No. 0910296 and 0905071). The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data. 57 T.-Q. Duong, C. Massobrio, M. Boero, G. Ori, E. Martin, Heat transport in disordered network forming materials: Size effects and existence of propagative modes, Computational Materials Science 177 (2020) 109607. doi:10.1016/j.commatsci.2020.109607.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

AUTHOR DECLARATIONS

The authors have no conflicts to disclose.