Évelyne Martin

Assil Bouzid

Guido Ori

Mauro Boero

Evelyne Lampin

Carlo Massobrio

Impact of dispersion forces on the atomic structure of a prototypical network-forming disordered system: The case of liquid GeSe2

Keywords: numbers: 61.25.Em, 61.20.Ja, 71.15.Pd

I. INTRODUCTION

The proper account of dispersive (van der Waals, vdW) contributions is a challenging task for atomic-scale computational material science. In this context, A x B 1-x (A = Ge, Si; B = O, Se, S) disordered network-forming materials are worth investigating, since made of atomic species qualitatively prone to be denoted as cations and anions and yet bearing close electronegativities. Such combination of physical effects results in cohesive forces of iono-covalent nature. The question arises on whether there might be, as a function of the composition x and for given families of compounds, bonding situations for which dispersion forces cannot be merely neglected by invoking efficient screening of other predominant bonding contributions.

Based on a well-established series of achievements, it appears that first-principles molecular dynamics (FPMD) in the framework of density functional theory (DFT) is a reliable tool to gain insight into the properties of glasses and liquids such as Ge x Se 1-x , Ge x Te 1-x , or their ternary counterparts [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] . For these systems, the number of results obtained without an explicit account of dispersion forces largely exceeds those for which such forces have been explicitly considered.

However, it is nowadays widely accepted that the mere use of the available exchange-correlation (XC) functionals cannot systematically ensure a proper account of long-range electron correlations responsible for van der Waals contributions. For instance, this was shown to be the case for liquid Ge 15 Te 85 , 17 pointing out the need of additional studies extended to other telluride and/or chalcogenide disordered materials.

Along these lines, and focusing on glassy GeTe 4 , the impact of dispersion forces and its interplay with the choice of the XC functionals within DFT was carefully analyzed. In Ref. 18, we compared the levels of agreement on structural properties between experiments and theory (FPMD), as predicted by a set of different modelling schemes. Four different theoretical recipes were considered, each one corresponding to a distinct structural model. These were obtained by selecting either the PBE (Perdew, Burke and Ernzerhof) 19 or the BLYP (Becke, Lee, Yang and Parr) 20,21 XC functionals and, for each one of the two, by disregarding or including van der Waals dispersion forces via the Grimme contribution. 22,23 The inclusion of the vdW forces lead to better performances and realistic structural predictions in the case of both the PBE and the BLYP exchange-correlation functionals.

In view of this result, the main motivation of the present work is to assess further the impact of dispersion forces on a representative case of disordered chalcogenide, liquid GeSe 2 . The presence of a predominant structural unit (the GeSe 4 tetrahedron) coexisting with homopolar bonds makes the atomic structure of liquid GeSe 2 challenging to elucidate for both experiments and atomicscale modelling. [24][25][26][27][28][29][30][31] As shown by FPMD, the description of both short and intermediate range order of liquid GeSe 2 depends on the choice of a specific exchange-correlation functional. 32,33 In particular, generalized gradient approximated (GGA) functionals favoring electronic localization bring calculations in better agreement with structural data. Substantial improvements were found when introducing the Perdew-Wang (PW) recipe at the place of the local density approximation (LDA) and, in a further step, the BLYP scheme at the place of the PW scheme. 20,21,[32][33][34][35] In the case of liquid GeSe 2 , the partial structure factors obtained via DFT-FPMD were able to reproduce experimental data for an extended range of wavevectors in reciprocal space, the same consideration holding for pair correlation functions in real space. 27,33 However, a number of differences between theory and experiments did persist, as the underestimate of the first sharp diffraction peak (FSDP) in the Ge-Ge partial structure factor. Such shortcoming is responsible for the absence of the FSDP in the Bhatia-Thornton concentration-concentration structure factor 36 and it was found to be largely independent on the periodic cell dimensions, as shown by taking a number of atoms N = 480 instead of N = 120. 37 Since the occurrence of such FSDP could not be ascribed to trivial system size effects, the consideration of other plausible origins for the residual lack of agreement with experiments is very much in order. Therefore, one of our main motivation is to focus on the dispersion forces as a missing ingredient in the theoretical description of this chalcogenide. 36 Our approach complements available indications on the impact of dispersion forces on Se/S made disordered chalcogenides, as those rationalized in Ref. 38. In the present work, we produce a new set of FPMD data on liquid GeSe 2 by including dispersion forces in the FPMD Kohn-Sham hamiltonian. Two different methodologies are employed.

One consists in following the formalism due to Grimme and, specifically, the D2 formula. 22 We also employed a scheme more profoundly rooted into the density functional theory in the sense that the electronic structure, obtained within the selected DFT-GGA framework, is used to compute the maximally localized Wannier functions (MLWF) on which the present van der Waals scheme is based. [39][40][41] Our results show that the account of dispersion forces does not alter significantly the structural properties of liquid GeSe 2 . However, the Grimme scheme is unable to reproduce correctly the experimental profile of the Ge-Ge pair correlation function, resulting in a perceptible underestimate of the number of corner-sharing connections. The level of agreement is restored when using the MLWF formalism, confirming its superior performances in terms of correct reproduction of physical properties closely related to chemical bonding. We found no impact of vdW forces on the concentration-concentration structure factor over the full range of wavevectors in reciprocal space.

This work is organized as follows. Section II provides details on our theoretical model, that is based on the first-principles molecular dynamics framework. Section III contains structural information extracted from the partial structure factors in the Faber-Ziman and Bhatia-Thornton formalisms. Section IV is devoted to the partial pair correlation functions and to a description of the network topology via the counting of the edge-sharing and corner-sharing connections.

Conclusive remarks are collected in Section V.

II. THEORETICAL MODEL

The theoretical framework of our calculations is identical to the one employed and detailed extensively in previous papers devoted to liquid GeSe 2 . 27,33 In the present case, FPMD simulations have been performed at constant volume on a system containing 240 (80 Ge and 160 Se) atoms positioned in a periodically repeated cubic cell of size 19.76 Å, this value allowing to recover the experimental density of the liquid 42 at the desired temperature (T = 1050 K). N = 240 turns out to be a suitable system size to obtain reliable structural properties of liquid GeSe 2 at both short and intermediate ranges. This choice is substantiated by our comparative work on system sizes N = 120 and N = 480, clearly demonstrating that all results obtained in the past with N = 120 were not significantly altered when moving to N = 480. 37 We recall that the electronic structure was described within density functional theory and evolved self-consistently in time by Car-Parrinello molecular dynamics (CPMD) 43 and the cpmd code. 44 The BLYP generalized gradient approximation was used, based on the exchange energy obtained by Becke 20 , and the correlation energy according to Lee, Yang and Parr. 21 Valence electrons have been treated explicitly, in conjunction with norm conserving pseudopotentials of the Trouiller-Martins type to account for core-valence interactions. 45 The wave functions have been expanded at the Γ point of the supercell on a plane wave basis set with an energy cutoff E c = 20 Ry. Long-range dispersion forces were accounted for in two distinct manners. The first consists in using the DFT-D2 formula proposed by Grimme 22 and it will be termed vdW1 in what follows. We stress the fact that, although empirical, such a vdW correction is a thorough DFT-based formulation in which the parameters are self-consistently tuned on different functionals, including the one used in this work, and benchmarked on a wealth of different systems from simple molecules to complex reactive surfaces and chalcogenides. [46][47][48] No experimental parameters are considered in the construction of this specific vdW correction and its inclusion does not affect at any stage the Kohn-Sham equations, 49,50 thus preserving the first-principle character of all electronic structure calculations. However, the Grimme formalism does not allow for any update of the vdW interactions as a result of the changes in the electronic structure that occur during the dynamical evolution. To go beyond this limitation, we employed a second vdW scheme in which the vdW coefficients are recalculated for each ion of the system, becoming explicitly time-dependent. [39][40][41] This scheme (referred to as vdW2 in what follows) makes use of the electronic wavefunctions via the unitary transformation of the Kohn-Sham orbitals giving the maximally localized Wannier functions. Such strategy has rapidly acquired an extended record of reliability for several test cases. 41,51,52 In the FPMD approach, a fictitious electron mass of 2000 a.u. (i.e. in units of m e a 2 0 where m e is the electron mass and a 0 is the Bohr radius), and a time step of ∆t = 0.12 fs have been used to integrate the equations of motion. The control of the temperature has been implemented for both the ionic and electronic degrees of freedom by using Nosé-Hoover thermostats. [53][54][55] The initial coordinates have been obtained by selecting from a periodic cell of N = 480 atoms (previously employed to produce the results of Ref. 37) the positions of the inner 240 atoms. As a first step, the system has been brought to T = 0 K by structural optimization via forces minimization, so as to avoid any artifact due to the presence of spurious interactions that might be due to our initial system construction. FPMD runs lasting in between 30 ps and 100 ps have been produced at increasing temperatures without accounting for vdW contributions to bring the system at the production temperature of T = 1050 K in the liquid phase. After collecting statistical averages over a period of 50 ps, the final configuration is stored to ensure the same initial starting point for simulations runs involving the vdW1 and vdW2 dispersion forces. In both cases, the liquid is first cooled at T = 300 K and then heated up over a time interval of 115 ps for vdW1 and 170 ps for vdW2. Statistical averages pertaining to the vdW1 and vdW2 situations were collected over 30 ps and 42 ps respectively. The different length of time for the collection of the statistical averages has no special significance, both trajectories being sufficiently extended to access thermal equilibrium. The statistical errors for the models employed have been found to be very similar.

We make inferences on the impact of the vdW1 and the vdW2 schemes for a set of properties in reciprocal and direct space. Our references for comparison are the experimental data of Ref. 24 obtained from isotopic substitution in neutron diffraction and the present FPMD results with no consideration of dispersion forces (NvdW hereafter). For all sets of calculations, unless otherwise stated, statistical errors are limited to 5% at most.

III. STRUCTURE FACTORS

A. Partial structure factors: Faber-Ziman and Bhatia-Thornton results

The calculated Faber-Ziman (FZ) partial structure factors for liquid GeSe 2 are shown in Fig. 1.

The FPMD sets of results are very close to each other in the case of the partial structure factors S FZ GeSe (k) and S FZ SeSe (k). However, the intensity of the main peak is significantly higher in the vdW1 case for S FZ SeSe (k). By focusing on the FSDP region around k ∼ 1 Å -1 , there is a vanishing impact of the vdW corrections on the intermediate range properties associated with Ge-Se correlations.

In the case of S FZ GeGe (k), there are some discernible variations in the peak intensities at k ∼ 1 Å -1 . Such differences are of little relevance in the context of the present investigation since mostly related to statistical errors comparable to those found when moving from N = 120 to N = 480.

Worth of notice for S FZ GeGe (k) are the intensities of the maximum and of the minimum located at k ∼ 2 Å -1 and k ∼ 3 Å -1 , respectively. These are more pronounced in the case of vdW1 when compared to the vdW2 and NvdW models, the three values lying within the experimental statistical errors. Taken together, it remains true that the analyses of the partial structure factors in the FZ formalism do not bring any striking evidence on the possible occurrence of structural changes due to dispersion forces. We recorded some minor effects in S FZ SeSe (k) and S FZ GeGe (k). It is of interest to seek further evidence on the impact of dispersion forces by considering the Bhatia-Thornton partial structure factors 56 S NN (k) (number-number), S NC (k) (numberconcentration) and S CC (k) (concentration-concentration) (Fig. 2).

These can be obtained by linear combinations of the Faber-Ziman partial structure factors as follows:

S NN (k) = c Ge c Ge S FZ GeGe (k) + c Se c Se S FZ SeSe (k) + 2c Ge c Se S FZ GeSe (k), (1)
S NC (k) = c Ge c Se [c Ge (S FZ GeGe (k) -S FZ GeSe (k)) -c Se (S FZ SeSe (k) -S FZ GeSe (k))], (2)
S CC (k) = c Ge c Se {1 + c Ge c Se [(S FZ GeGe (k) -S FZ GeSe (k)) +(S FZ SeSe (k) -S FZ GeSe (k))]}. (3)
where c Ge and c Se denote respectively the atomic fraction of Ge and Se.

We recall that due to very close values of the coherent scattering length of the chemical species Ge and Se (b Ge = 8.185 fm, b Se = 7.97 fm) 57 and to the limited range of variation of S NC (k) and S CC (k), S NN (k) is a very good approximation of the total structure factor, i.e.

|S T (k)-S NN (k)| < 0.025. Therefore, the considerations developed hereafter on S NN (k) apply equally well to the total neutron structure factor S T (k).

There are no discernible effects of the dispersion forces on the Bhatia-Thornton partial structure factors S NC (k) and S CC (k). In particular, one notice the persistent absence (or the large underestimate) of the FSDP in all calculated S CC (k) strongly contrasting with the prominent feature found in neutron diffraction experiments. 24 Given these results, one can safely rule out any impact of vdW forces on the intensity of the FSDP, thereby adding another piece of information toward precise understanding of this longstanding open issue. The Bhatia-Thornton partial structure factors S NN (k) features higher peak intensities at k ∼ 2 Å -1 and k ∼ 3.5 Å -1 when adopting the vdW1 recipe. Although these variations are minimal and within the statistical error of the experimental measurements, this is in line with the previous observation of the S FZ GeGe (k) and S FZ SeSe (k) cases, indicating that the Grimme approach deviates more from the experimental pattern than the MLWF scheme.

IV. REAL SPACE PROPERTIES

A. Partial pair-correlation functions

Fig. 3 provides a comparison between the three sets of calculations NvdW, vdW1 and vdW2 of the pair correlation functions g GeGe (r), g GeSe (r) and g SeSe (r), and the corresponding experimental results. Integration of the pair-correlation function g GeSe (r) over the first shell of neighbors leads to coordination numbers nGeSe lying within 6 % ( nGeSe (exp) = 3.5, nGeSe (NvdW) = 3.7, nGeSe (vdW1) = 3.7, nGeSe (vdW2) = 3.6). These numbers exemplify the lack of sensitivity of Ge-Se correlations to the presence of dispersion forces, in agreement with the predominant ionic contribution that can be assumed to be responsible of most of the Ge-Se bonding.

In the case of g SeSe (r) the experimental value for nSeSe is nSeSe (exp) = 0.23. With our calcula-tions we obtain nSeSe (NvdW) = 0.20, nSeSe (vdW1) = 0.19, nSeSe (vdW2) = 0.25. Consideration of these numbers and visual inspection of Fig. 3 are indicative of a broad level of agreement among the three sets of data. However, the Grimme recipe (vdW1) stands out as the farthest apart from the experimental prediction when focusing on the second main peak at about 4 Å.

The case of g GeGe (r) is worth dwelling on in deeper detail. This partial pair correlation function has been recognized as the one featuring the least satisfactory agreement with experiments, despite an overall quite acceptable behavior in terms of coordination numbers due to the use of the BLYP exchange correlation functional. 33 The partial pair correlation function g GeGe (r) features a three peak structure for r < 4 Å, representative of homopolar bonds, edge-sharing (ES) and corner-sharing (CS) connections. The pattern corresponding to these marks is reproduced by all FPMD models, with unequal intensities of the first peak corresponding to homopolar bonds.

These differences are partially compensated by a broad agreement on the values of the coordination numbers nGeGe (NvdW) = 0.18, nGeGe (vdW1) = 0.19, nGeGe (vdW2) = 0.23 all of them lying close to the experimental result, 0.25. We recall that sizeable variations on the intensity and width of the first peak are somewhat expected in view of the extended temporal persistence of Ge-Ge homopolar bonds that might form on specific trajectories. For the purpose of the present work, it is of interest instead to focus on the heights of the second and of the third peak of g GeGe (r). The relative intensities of these peaks is clearly different in the vdW1 case, with a third peak standing out well beyond any expected statistical error. This corresponds to an overestimate of corner-sharing connections induced by the Grimme approximation. As a valuable alternative, MLWF succeed in leaving unaltered the quite acceptable profile of the calculated g GeGe (r), thereby providing a first, clear-cut evidence of its better performances for this kind of network structures.

B. Bond angle distributions and network topology

In Fig. 4, we show the Se-Ge-Se (θ SeGeSe ) and Ge-Se-Ge (θ GeSeGe ) bond angle distributions.

The shapes of θ SeGeSe for the three cases are nearly identical and highly symmetrical, pointing out a predominant tetrahedral arrangement, with angles distributed around the expected value 109 • . In this respect, the dispersion forces do not alter the regularity of the tetrahedral structures. Turning to θ GeSeGe , the two distinct peaks visible at about 80 o and 100 o can be assigned to the formation of edge-and corner-sharing tetrahedral connections, respectively. 27,28,33 Changes in the edge-sharing/corner-sharing peak intensities ratio suggest that corner-sharing connections are more numerous in the vdW1 case, while the vdW2 scheme has little effect on the relative proportion of these quantities when compared to NvdW results. To quantify these observations we have calculated the number of Ge atoms belonging to corner-sharing and edge-sharing connections, N Ge (CS) and N Ge (ES) respectively, together with the number of Ge atoms found in homopolar bonds N Ge-Ge . To this end we adopted the proposal of Ref. 58, i.e. N Ge (CS) = 1-N Ge (ES)-N Ge-Ge , which holds in the absence of extended chains. 58 By using the values of N Ge (CS) and N Ge (ES) for the three models NvdW, vdW1, vdW2 (see this meaning that the number of CS connection cannot be higher than the number of ES connections. Even though the values obtained for CS and ES in the vdW1 case are not too dissimilar when accounting for the statistical uncertainties, it remains true that a tendency to overestimate the CS/ES ration is noticeable. Therefore, the result corresponding to vdW1 for R vdW 1 CS/ES (1.1) is indicative of an artificial and yet, moderate contribution introduced in the FPMD model by the Grimme recipe. This contribution favors the CS arrangement and longer Ge-Ge distances for the inter-tetrahedral connections.

V. CONCLUSION

Liquid GeSe 2 is a disordered network-forming material well suited to be exploited as benchmark system for model improvements and assessment of performances. In this case, we are interested in the impact of dispersion forces on the structural properties of this network, due to the existence of such effects in other disordered chalcogenide materials. Also, dispersion forces are a missing ingredient in the theoretical treatment of this prototypical liquid, for which numerous investigations by using different first-principles atomic-scale models have become available. To account for van der Waals dispersion forces, we have employed two schemes. In the first, due to Grimme, the coefficient expressing the intensity of the long range attractive forces does not evolve (i.e. is not updated) following the temporal evolution of the electronic structure. In the second, this requirement is met by resorting to the maximally localized Wannier functions. Overall, the set of structural properties considered here (partial structure factors and pair correlation functions, as well as standard fingerprints of network topology as the coordination numbers and the kind of the tetrahedral connections) are not significantly affected by the account of the dispersion forces.

In particular, the maximally localized Wannier scheme leads to structural data very close to those obtained in the absence of dispersion forces. On the contrary, the Grimme approach introduces some spurious contributions clearly noticeable at the level of the partial pair correlation functions and the Ge-Se-Ge bond angle distribution. In terms of structural features, these correspond to a moderate overestimate of the number of corner-sharing connections at the expenses of the edgesharing connections. Our results provide convincing evidence on the superior reliability of the MLWF scheme in accounting for dispersion forces. It appears of paramount importance to rely on a methodology that does not modify structural properties when the dispersion forces have a negligible effect in comparison to other bonding effects. Further studies of this kind for disordered network systems at concentrations differing markedly from the stoichiometric one would be instrumental to confirm this conclusion.

Figure 1 .

 1 Figure 1. (Color online) The Faber-Ziman partial structure factors S FZ GeGe (k) (top panel), S FZ GeSe (k) (middle panel) and S FZ SeSe (k) (bottom panel) for liquid GeSe 2 at T = 1050 K. Experimental data from Ref. 24 (crosses with error bars) are compared to results for the models NvdW (green line), vdW1 (red line) and vdW2 (blue line). The results of calculations are shown by including (lower curve) and by not including (higher curve) the experimental results in the comparison.

Figure 2 .

 2 Figure 2. (Color online) The Bhatia-Thornton partial structure factors S NN (k) (top panel), S NC (k), (middle panel) and S CC (k) (bottom panel) for liquid GeSe 2 at T = 1050 K. Experimental data from Ref. 24 (crosses with error bars) are compared to results for the models NvdW (green line), vdW1 (red line) and vdW2 (blue line). The results of calculations are shown by including (lower curve) and by not including (higher curve) the experimental results in the comparison.

Figure 3 .

 3 Figure 3. (Color online) The partial pair correlation functions g GeGe (r) (top panel), g GeSe (r) (middle panel) and g SeSe (r) (bottom panel) for liquid GeSe 2 at T = 1050 K. Open small circles are the experimental data

Figure 4 .

 4 Figure 4. (Color online) The calculated bond-angle distributions θ GeSeGe and θ SeGeSe for liquid GeSe 2 at T = 1050 K. The results are given for the models NvdW (green line), vdW1 (red line) and vdW2 (blue line).These distributions have been calculated by including neighbors separated by less than 3 Å.

  Tab. I) one obtains R NvdW CS/ES = N Ge (CS,NvdW)/N Ge (ES,NvdW) = 0.8, R vdW 1 CS/ES = N Ge (CS,vdW1)/N Ge (ES,vdW1) = 1.1 and R vdW 2 CS/ES = N Ge (CS,vdW2)/N Ge (ES,vdW2) = 0.87. The experimental prediction of Ref. 24 points toward close percentages of edge-sharing and corner-sharing connections. Likewise, all FPMD models produced in the past are in line with values for R CS/ES moderately smaller than 1,

Table I .

 I Percentages of Ge atoms in edge-sharing, corner-sharing and homopolar bonding configurations

		N Ge (ES)	N Ge (CS)	N Ge-Ge
	l-GeSe 2 (NvdW)	45	36	18
	l-GeSe 2 (vdW1)	39	43	18
	l-GeSe 2 (vdW2)	41	36	23


Assil Bouzid

Hayat Zaoui

Pier Luca Palla

Guido Ori

Mauro Boero

Carlo Massobrio

Fabrizio Cleri

Evelyne Lampin

Thermal Conductivity of glassy GeTe 4 by First-Principles Molecular Dynamics

A transient thermal regime is achieved in glassy GeTe 4 by first-principles molecular dynamics following the recently proposed "approach-to-equilibrium" methodology. 1 The temporal and spatial evolution of the temperature do comply with the time-dependent solution of the heat equation. We demonstrate that the time scales required to create the hot and the cold parts of the system and observe the resulting approach to equilibrium are accessible to first-principles molecular dynamics. Such a strategy provides the thermal conductivity from the characteristic decay time. We rationalize in detail the impact on the thermal conductivity of the initial temperature difference, the equilibration duration, and the main simulation features.

Molecular dynamics (MD) simulations of thermal properties, e. g. thermal conductivity and interface thermal resistance, give access on equal footing to all the orders of phonon-phonon scattering via the anharmonicity of the interatomic energy landscape, together with relevant atomicscale details such as surface roughness, point or extended defects, pores or nano-inclusions and interfaces. Using interatomic potentials to describe the interactions, classical molecular dynamics (CMD) has been applied to the calculation of thermal conductivity either via the equilibrium molecular dynamics (EMD) (Green-Kubo 2 ) method or the so-called direct method. 3 As a viable and faster alternative, we have developed the approach-to-equilibrium molecular dynamics (AEMD) 1,4 based on the establishment and analysis of thermal transients. This strategy does not require a stationary state to be attained, all information leading to the description of the thermal process being readily available within the molecular dynamics trajectory of the transient. A gain of one order of magnitude in computational time can be estimated for AEMD with respect to stationary approaches.

It appears that the description of the energy landscape and its anharmonic character might not be properly addressed by empirical potentials, resulting in poor predictions of the thermal conductivity. For instance, in crystalline silicon, the thermal conductivity can be largely overestimated. [5][6][7] The situation can be even worse in the case of interfaces between organic and inorganic materials 8 , when the anharmonic regime of the potential energy is sampled by atomic vibrations far from equilibrium, at high pressure/temperature 9 , or in the case of complex materials like glasses. [START_REF] Boero | Springer Series in Materials Science 215: Molecular Dynamics Simulations of Disordered Materials. From Network Glasses to Phase-Change Memory Alloys[END_REF] To go beyond these limitations and achieve quantitative predictions, a first-principle description of thermal transport via density-functional theory (by using FPMD, first-principles molecular dynamics) is highly desirable. However, FPMD is often characterized by a prohibitive computational cost, limiting typical applications to several hundreds of picoseconds at most for simulations cells hardly containing more than 500 atoms. So far, FPMD has been used in combination with lattice dynamics, under the approximations employed to treat the thermal conductivity of a kinetic gas. 9 The question arises on whether the available stationary approaches are well suited for the study of heat transport within a fully dynamical FPMD scheme. In the case of the direct method, the supercell dimensions have to be carefully determined either to converge to a bulk thermal conductivity or to allow approaching the convergence via a reliable extrapolation. As for the Green-Kubo method, the definition of energies and currents quantummechanically is far from straightforward requiring specific methodological efforts [START_REF] Marcolongo | [END_REF] . For these reasons, one would like to turn to AEMD to circumvent these stumbling blocks and obtain a quantitative assessment of thermal transport properties.

On the side of transient methods, early attempts to study a thermal regime by FPMD are due to Luo and Lloyd. 12 These authors have considered crystalline structures of Si, Ge and SiGe containing 64 atoms, showing that a temperature difference between two halves of the system can be obtained by applying a double thermostat. Having in mind AEMD as a viable transient method worth applying, these results are quite instructive since for Si and Ge the distribution of phonon mean free paths (MFPs) extends above the micrometer. 13,14 Therefore, for Si and Ge, no meaningful comparison between experimental thermal conductivities and those stemming from such small supercells can be obtained. Accordingly, to be feasible and rewarding, the applications of AEMD within FPMD have to be not only computationally affordable but also fully compatible with the characteristic spatial lengths inherent in the thermal processes.

In the present letter we achieve this goal for glassy GeTe 4 (g-GeTe 4 ), a system widely considered in the context of phase change memory devices. In this material, the heat carriers have very short mean free paths, namely 4-5 Å in the temperature range [100-200] K. 15 This means that the minimal box lengths allowing for the study of the thermal processes are largely within the reach of FPMD, thereby granting legitimacy to the application of a transient method like AEMD.

The structure of g-GeTe 4 has been recently studied 16 by Car-Parrinello molecular dynamics (CPMD) 17,18 providing atomic models where the characteristic structure factor and correlation function are in excellent agreement with experimental measurements. The pristine structure of g-GeTe 4 adopted in this work contains N at =185 atoms in a cubic box with a side L =18 Å (see Fig. 1) and it was obtained using the CPMD approach. Within the framework of a study devoted to the impact of different exchangecorrelation functionals, 16 with or without account of the dispersion forces, two models of g-GeTe 4 were produced by adopting the Perdew, Burke, and Ernzerhof (PBE) 19 and the Becke, Lee, Yang and Parr (BLYP) 20,21 exchange-correlation functionals, respectively. The two schemes were found to have similar performances, the BLYP one being moderately closer to experimental data. The complete set of details for the simulation of g-GeTe 4 are available in Ref. 16 . We used the BLYP exchange-correlation functional with a norm-conserving pseudo-potentials 22 describing the valence-core interaction. Long-range dispersion forces are considered by following the formula by Grimme. 16,23 Thermostats are applied to the fictitious electronic degrees of freedom by following the Blochl and Parrinello guidelines. 24 We have applied the AEMD methodology to g-GeTe 4 , with a target temperature of T t = 130 K since the MFPs are measured to be the lowest around this temperature. 15 We started by equilibrating the simulation box at T t . Then, we divided the simulation box into two sub-parts along the ydirection (Fig. 1). These blocks are intended to be the cold and hot parts allowing for the establishment of the thermal transient. Each block contains around 90 atoms. During a first phase, two distinct Nosé-Hoover 25,26 thermostats are applied to the blocks, with a temperature difference ∆T 0 equal to 50, 100 or 200 K. This step does not require any code development since the multiple thermostats option is already implemented in the CPMD code. Fig. 2 shows the time evolution of the temperature of the hot and the cold blocks for a temperature difference of 200 K. The temperature gap establishes in a few ps and it is stabilized by the application of the two thermostats along a time slot lasting 35 ps. Then, the thermostats are switched off to allow for the second phase of AEMD to begin. As a consequence, the temperature of the hot (cold) block decreases (increases). Both reach 130 K in 10 ps.

The temperature profile and the temporal evolution of the temperature difference during phase 2 of AEMD are analyzed to check whether they comply with the transient solution of the heat equation. In the conditions of the AEMD simulations, the solution T (t, y) is a Fourier series 1,27 :

T (t, y) = T t + ∞ ∑ m=0 2∆T 0 (2m + 1)π sin 2π(2m + 1) L y e -(2m+1) 2 t/τ .
(1) The dominant contribution of the Fourier series (m=0) is a decaying sine :

T (t, y) -T t ∝ ∆T 0 sin 2πy L e -t/τ . (2)
The temperature difference between the two blocks that corresponds to:

∆T (t) = 2 L L/2 0 T (t, y)dy - 2 L L L/2 T (t, y)dy (3)
has a dominant contribution of the form:

∆T (t) ∝ ∆T 0 e -t/τ . (4)
The temperature difference, shown in the inset of Fig. 3 follows Eq. 4 as exemplified by the linear variation on the semi-log graph. The decay time τ is extracted by a suitable fit. The temperature profile, averaged along 2 ps (in between 37 and 39 ps of the total simulation time) follows Eq. 2. Since both the temporal evolution of the temperature difference and the temperature profile are consistent with the transient solution of the heat equation, we can exploit the relation between the decay time and thermal conductivity derived from the heat equation:

1 κ = L 2 4π 2 C V .ρ τ ( 5
)
where ρ = N at /L 3 is the number density. The specific heat C V = 3k B (0.78 ± 0.02) is calculated from the energy fluctu- ations per atom ∆E in a {NV T } run at 130 K: 28

C V = ∆E 2 k B T 2 (6)
The thermal conductivity is equal to 0.02 W m -1 K -1 for the corresponding best fit value of τ = 4 ps (Fig. 3). It is of interest to ascertain the impact of parameters like the initial temperature difference or the duration of phase 1 on the value of the decay time τ. The results are presented in Fig. 4. For each configuration, ∆T (t) is fitted to an exponential decay and five different lengths of the time interval are adopted to obtain averages of ∆T (t). The resulting values of τ (and the standard deviations) are given in Fig. 4. On the right axis, τ is converted into a thermal conductivity using Eq. 5.

For instance, in configuration 1, used to obtain the results of Figs. 2 and3, we applied a single thermostat at T = 130 K to the whole box, before starting the phase 1 with a double thermostat. The difference of the thermostat temperatures is set to ∆T 0 : 50, 100 or 200 K. The magnitude of the temperature decay is considerably reduced with ∆T 0 =50 K, affecting the precision in the determination of τ. Therefore, we have taken either ∆T 0 =100 or 200 K in the other configurations. In configuration 2, the whole box is first entirely heated at T f + ∆T 0 /2, before a half is cooled at T f -∆T 0 /2. In this way, the temperature fluctuations are reduced in the hot block at the end of phase 1. To ensure that our results do not depend significantly on the simulation framework, we have also considered the impact of the Nosé thermostat associated with the fictitious electronic degrees of freedom by switching it off during the phase 2 (configuration 3). Finally, in configuration 4, the frequency of the ionic Nosé thermostats ω i is increased from 200 to 350 cm -1 in order to optimize the onset of the temperature difference between hot and cold parts. Fig. 4 shows that the decay time takes values included between 4.0 and 10.5 ps without any clear dependence on the different parameters. The corresponding range for the thermal conductivity is [0.008:0.021] W m -1 K -1 . The thermal conductivity averaged over all configurations standing for the different simulation conditions is equal to 0.013 ± 0.003 W m -1 K -1 . This averaged value is far from being unrealistic, as it can be understood from the following considerations involving mostly experimental pieces of evidence and (a few) possible details inherent in our simulation approach. Along these lines, it is worth mentioning that Zhang et al. 15 obtained an experimental value of 0.14 W m -1 K -1 at 130 K using a parallel temperature conductance technique, reportedly leading to an overestimate by 20-25 % of the thermal conductivity compared to a laser flash measurement. This latest experimental framework is closer to our simulation methodology. As an additional issue, the structure and chemical order of such disordered systems are known to be strongly dependent on the conditions of synthesis and processing, resulting in network structures differing by the amount of chemical disorder. A reduced chemical disorder acts to increase the thermal conductivity up to a maximum reached for an ordered crystalline network. For example, in GeTe films, the thermal conductivity goes from 0.23 W m -1 K -1 for the glass to 3.08 W m -1 K -1 for the crystal. [START_REF] Fallica | [END_REF] Zhang et al. 15 claim that the glassy state of GeTe 4 is unstable, and the materials has a tendency to crystallize. Again, in the presence of grains, the overall thermal conductivity is likely to be larger than the one pertaining to the glass. On the FPMD simulation side, enhanced chemical disorder is likely to happen since a tendency to overestimate the departures from heterogeneous bonding have been found to characterize chalcogenide glasses, [30][31][32] in a way that could contribute to obtain lower values for the thermal conductivity. As a final remark on modelling, one cannot neglect the possible impact of size dependent effects due to heat carriers with mean free paths larger than the present box length. Future calculation on large cell sizes will be highly instrumental in this context.

In conclusion, we have exploited the approach-to-equilibrium methodology to calculate the thermal conductivity of a glass with first-principles molecular dynamics in the framework of density functional theory. The setup of the methodology is tractable and straightforward, without requiring any additional implementation in the first-principles CPMD code. We have shown that it is possible to create in a few ps a hot and a cold block by using a double thermostat, and this for an affordable number of atoms (N at =185). The approach-to-equilibrium takes a short time, around 10 ps. Both the decay of the temperature difference during this second phase, and the related temperature profile do comply with the heat equation. This allows an access to the thermal conductivity via the decay time. In the specific case of glassy GeTe 4 , the present approach leads to a calculated thermal conductivity lower than the experimental value. This difference can be rationalized in terms of uncertainties in the experimental structural determination, different amounts of chemical disorder in the corresponding experimental vs. theoretical structures, and specific modelling issues. However, none of them affect the validity and impact of the present approach that proved to be a powerful tool to better understand and quantify thermal phenomena.

a

  Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504, F-67034 Strasbourg, France. † Present address: Chaire de Simulation à l'Echelle Atomique (CSEA), Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland. b Univ. Lille, CNRS, Centrale Lille, ISEN, Univ. Valenciennes, UMR 8520 -IEMN, F-59000 Lille, France. E-mail: evelyne.lampin@univ-lille1.fr

Fig. 1

 1 Fig. 1 Simulation box containing 185 atoms of amorphous GeTe 4 .

Fig. 2

 2 Fig. 2 Temperature in each block during the two AEMD phases.

Fig. 3

 3 Fig.3 Temperature profile and time evolution (from the beginning of phase 2) of the temperature difference (inset) between hot and cold blocks. The dark lines are the CPMD results. The red line is a sine curve adjusted to the CPMD results. The green line in the inset is a decaying exponential adjusted to the CPMD results.

Fig. 4

 4 Fig.4 Sensitivity of the decay time τ in phase 2 and thermal conductivity to the parameters of the simulation. In blue, the initial temperature difference is equal to 200 K, in green to 100 K and in red to 50 K. Configuration 1: (+), 2: (o), 3: ( ), 4: ( ).

J o u r n a l N a me , [ y e a r ] , [ v o l . ] , 1-5 | 1

J o u r n a l N a me , [ y e a r ] , [ v o l . ] , 1-5 | 4

Acknowledgements

Calculations were performed on the facilities of IEMN (clustphy) and on the computational Mesocenter of the University of Strasbourg.


Assil Bouzid

Guido Ori

Mauro Boero

Evelyne Lampin

Carlo Massobrio

Atomic-scale structure of the glassy Ge 2 Sb 2 Te 5 phase change material: A quantitative assessment via first-principles molecular dynamics

I. INTRODUCTION

Ge 2 Sb 2 Te 5 (GST hereafter) is the most established and successful phase change material (PCM) [1,2] due to its high switching speed and its remarkable optical contrast between the amorphous and the crystalline phases. To exploit the properties of GST and optimize applications, a precise understanding of its atomic scale network and bonding properties is highly desirable. This can be highly beneficial to unveil the cause of current technological limitations involving resistance drift and threshold switching phenomena [3,[START_REF] Wong | Proc. IEEE[END_REF]. Despite the known dependence of the local order of amorphous GST on the sample processing [START_REF] Carria | [END_REF]6], the structure of the GST network (mostly octahedral, tetrahedral, or containing both) and its bonding nature remain elusive [7,8]. Experimentally, x-ray absorption fine structure (EXAFS) and x-ray absorption near-edge structure (XANES) measurements [9] point toward predominant tetrahedral Ge motifs, in contrast to the octahedral arrangement of the crystalline phase. Further refinement by EXAFS and bond constraint theory [10,12] revealed a high concentration of Ge-Ge homopolar bonds, along with overcoordinated Te sites. However, 125 Te nuclear magnetic resonance (NMR) spectroscopy [13] proposed a mixed Te coordination environment that consists solely of heteropolar Ge/Sb-Te bonds. Also, x-ray fluorescence holography performed on an epitaxial GST layer highlighted the occurrence of a cubic structure with a tetrahedral symmetry around Ge sites [14]. The displacements of Ge atoms from octahedral to tetrahedral sites was found at the very origin of a nonthermal crystal-to-amorphous phase transition via ultrafast time-resolved electron diffraction experiments [15].

The question arises on the capability of a modeling scheme based on density functional theory (DFT) to elucidate the atomic structure of GST in a way consistent with experiments. We are referring here to the structural optimization or evolution of temporal trajectories at finite temperatures, as performed within first-principles molecular dynamics (FPMD). While the amorphous state [16] was first rationalized in terms of a spinal-like geometry (tetrahedral geometry), subsequent FPMD simulations revealed the coexistence of tetrahedral and octahedral Ge sites [17][18][19], in puzzling contrast with EXAFS data. Overall, it appears that FPMD approaches have been instrumental to unravel the fundamental features of several processes inherent in the behavior of GST (as the crystallization mechanism, for instance). However, available FPMD models overestimate Ge-Ge, Ge-Te, and Sb-Te bond lengths obtained by EXAFS measurements leading to not entirely satisfactory predictions for neutron structure factors and pair correlation functions [17][18][19][20][START_REF] Akola | European Phase-Change and Ovonics Symposium[END_REF].

This observation prompts the search of alternative approaches to describe glassy GST at the atomic scale. Such motivation is substantiated by recent achievements obtained for glassy GeTe 4 [START_REF] Bouzid | [END_REF] using the Becke-Lee-Yang-Parr [23,24] (BLYP) exchange-correlation (XC) functional and appropriate pseudopotentials (PPs) [25].

In this paper, we exploit the modeling framework of Ref. [START_REF] Bouzid | [END_REF] to carry out FPMD calculations of glassy Ge 2 Sb 2 Te 5 . Our results are in very good agreement with x-ray [9,10], neutron scattering [26,27], and EXAFS [10,14,28] measurements, as shown by comparison with published FPMD data.

The present paper is organized as follows. In Sec. II, we provide details of our calculation methodology, built on the use of two different XC functionals and well calibrated choices for the PPs and their construction. Section III is devoted to our results and it is organized in subsections. The first features the comparison between theory and experiments for the total structure factor S(k), followed by the analysis of the total and partial pair correlation function g(r). Then, we move to the coordination numbers. In a final section, we highlight correlations between structural and electronic properties by relying on the Wannier functions formalism. Conclusions are contained in Sec. IV.

II. METHODOLOGY OF CALCULATIONS

We adopted the Car-Parrinello [29] approach (CPMD) [30] by making use of the BLYP XC functional or the Perdew-Burke-Ernzerhof (PBE) XC functional [31]. The valence-core interaction is described by either norm-conserving Troullier-Martins (TM) PPs (PBE-TM and BLYP-TM) or by normconserving Goedecker, Teter and Hutter (GTH) PPs [32,33] (BLYP-GTH). Valence electrons are represented by a planewave basis set expanded at the point of the simulation cell with an energy cutoff of 30 Ry. A fictitious electron mass of 500 a.u. and a time step t = 0.12 fs ensured optimal conservation of the constant of motion. The ionic temperature was controlled with a Nosé-Hoover [34][35][36] thermostat chain [37], whereas for the fictitious electronic kinetic energy we used a Blöchl-Parrinello thermostat [38] with the target kinetic energy set to 0.04 au. We adopted the Kleinman-Bylander (KB) construction [39] based on a local component V loc (r) and a pseudoatomic eigenstate |φ l,m expansion for the nonlocal term. We remind that the Kleinman-Baylander construction of PPs is based on a separable form where we distinguish between a local and a nonlocal part:

V KB = V loc (r) + l max lm |V l (r)φ lm V l (r)φ lm | φ lm |V l (r)|φ lm , ( 1
)
where V loc (r) is the local part of the potential, V l (r) = V PS (r) -V loc (r) is the nonlocal part of the potential, and V PS (r) is the atomic potential of the reference state. φ lm denotes the eigenstate of the atomic pseudo-Hamiltonian. While the choice of the local part V loc (r) is arbitrary, one typically truncates the sum over the angular momentum lm at a given value l max . BLYP-GTH PPs were found insensitive to the KB construction choice. In the case of TM PPs, for both Sb and Te, our choice for the local part and for the maximum angular momentum component (l max ) of the nonlocal one are loc = p and l max = d. This is identical to what is commonly employed in the literature. Instead, for Ge, we made an alternative selection, namely, loc = p, l max = p. This option features an extended record of reliability for a large variety of disordered chalcogenide materials [40][41][42][43][44].

The initial configuration, common to all simulations, consisted of a disordered system of 144 atoms (32 Ge, 32 Sb, and 80 Te) in a cubic simulation cell of side 16.86 Å, corresponding to a density of 0.030 atom Å-3 , on which periodic boundary conditions were applied. NVT simulations were performed at T = 300 (40 ps), 600 (20 ps), 900 (70 ps), 600 (50 ps), and 300 (42 ps) K. The overall features of the network topology (tetrahedral versus octahedral coordination) are not substantially modified when releasing the residual pressure on the initial volume via a slight expansion. Statistics were collected on this last part of the thermal cycle over the last 30 ps. We remark that at 900 K the system shows a high mobility and broad pair correlation functions typical of a liquid state. Given these conditions, atomic displacements covered distances of several atomic bond lengths ensuring effective loss of memory of the initial configuration. Three systems were produced, differing by the XC functional and the type of valence-core PPs. Namely, a first system generated within a BLYP XC functional and TM PPs (BLYP-TM), a second one within a PBE XC functional and TM PPs (PBE-TM), and a third one (BLYP-GTH) within the identical BLYP XC functional of the first case but making use of GTH PPs (BLYP-GTH).

To establish correlations between the structural features of the GST network and its chemical bonding properties, we calculated maximally localized Wannier functions [45,46] averaged over 200 configurations sampled along the last 25 ps. The analysis of the Wannier centers with respect to the nuclear positions proved instrumental in gaining insight into the correlation between structural and bonding properties for several disordered chalcogenides [42,47,48]. The disagreement with experiments is much more severe in the BLYP-GTH case. Knowing that the performance of a given model is at its best for smaller values of the Wright parameter R χ [49] [Eq. ( 2)], BLYP-TM features values of R χ roughly half (R χ equal to 6.84%) of those found previously by FPMD-DFT models,

R χ ≡ i [S exp (k i ) -S FPMD (k i )] 2 i S 2 exp (k i ) 1/2 . ( 2
)
Here, S exp (k i ) and S FPMD (k i ) are the experimental and the simulated structure factors at a given wave vector k i . Worth of note is also the R χ values for the PBE-TM case, equal to 7.31. Figure 2 exemplifies the different performances of the DFT schemes employed, by focusing separately on (i) those for which the agreement on the total neutrons structure factor is less satisfactory [Fig. 2(a)] and (ii) on our BLYP-TM results [Fig. 2(b)]. All theoretical models are compared to neutron diffraction measurements from Ref. [27]. The results considered in Fig, 2(a) refer to our BLYP-GTH model and to those of Ref. [18] [Born-Oppenheimer simulation (BOMD) with PBE XC functional and TM PPs, called hereafter BOMD-PBE-TM] and Ref. [19] [second generation Car-Parrinello simulation (SGCPMD) with PBE functional and GTH PPs, called hereafter SGCPMD-PBE-GTH]. We found R χ equal to 12.93% for BLYP-GTH, to 12.68% for BOMD-PBE-TM, and to 12.09% for SGCPMD-PBE-GTH. The origin of the I, we use l max = p as the maximum channel in the KB construction of the Ge PP, while the TM PP from Ref. [18] as well as our GTH and the GTH scheme from Ref. [19] have d as the maximum channel in the Ge PP.

l max = p l max = d l max = d TM Ref. [18] l max = d l max = d l max = d GTH this work l max = d l max = d l max = d GTH Ref. [19] l max = d l max = d l max = d lower R χ (6.
Given these premises, the question arises on the hypothesized impact of this construction on the structural properties of glassy GST. To assess this point, an amorphous GST system has been generated following the same recipe adopted for BLYP-TM and taking the Ge PP with l max = d (this model is called hereafter BLYP-TM-D). The Ge electronic structure is [Ar] 3d 10 4s 2 4p 2 , thus the d channel included in the pseudopotential is empty. Results are compared to experiments and to the BLYP-TM case (l max = p) in Fig. 1.

For the BLYP-TM-D case, we find a Wright parameter R X equal to 10.95%. This value is similar to the one found with the BLYP-GTH, BOMD-PBE-TM and SGCPMD-PBE-GTH schemes. Overall, it appears that the FPMD performances in the case of glassy GST acquire a stronger quantitative character when adopting l max = p in the Ge PPs. This legitimates the choice adopted throughout the present paper.

Partial pair correlation functions (PCFs) are reported in Fig. 3, while related coordination numbers are listed in Table II. These are determined using as a cutoff distance criterion the first minimum of each g αβ (r), apart from the Sb-Te bonds, where a cutoff radius of 3.4 Å was selected on the basis of the outer edge of the Sb-Te PCF in the crystal phase. For completeness, a comparison of the partial pair correlation functions obtained via BLYP-TM and BLYP-TM-D is provided in Fig. 4. The first peak in both the Ge-Ge and Ge-Sb PCFs is assigned to the presence of Ge-Ge and Ge-Sb "wrong bonds." Following the terminology of Ref. [51], we define as wrong bonds those that do not occur in the crystal. At variance with former models obtained by typical melt-quench procedures [19,[START_REF] Akola | European Phase-Change and Ovonics Symposium[END_REF] our results show a considerable amount of "wrong bonds." In this respect, our results are in line with the estimate obtained on the basis of experimental high energy x-ray and K-edge EXAFS measurements [28]. While the PBE-TM model slightly overestimates bond lengths

B. Coordination number and analysis of local environment

The distribution of the coordination numbers for Ge, Sb, and Te are shown in Fig. 5. The local environment turns out to be rather similar for both the BLYP-TM and the PBE-TM systems, the difference being less than 5%, whereas significant differences are found for the BLYP-GTH case. Glassy Ge 2 Sb 2 Te 5 exhibits a large amount (∼80%) of fourfold Ge sites. Among these fourfold atoms, Ge-Ge n Te 4-n (n = 0,1) are the most abundant, i.e., 38.6% for BLYP-TM and 38.9% for PBE-TM. Sb is mainly threefold (∼60%) with ∼30% of fourfold configurations and Te is mainly twofold (∼60%) with ∼30% of threefold atoms. Such a topology provides a new picture of glassy GST, corresponding to the best agreement recorded so far with experimental data.

Clear cut information on the atomic structure can be obtained from the order parameter [52] 2 . Such parameter is able to quantify the different weight of distinct coordination, the sum running over the bonded atom pairs where j is the central atom forming a bonding angle θ ij k with its neighbors. For a perfect tetrahedron, the order parameter is q = 1, while q = 0 refers to a sixfold octahedral environment. Defective octahedra correspond to TABLE 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 models generated with BLYP and PBE XC functionals. α and β denote the two atomic species. T Ge is the fraction Ge atoms in tetrahedral geometry. Experimental results are from Refs. [9,10]. SGCPMD stands for "second generation Car-Parrinello molecular dynamics," while BOMD stands for "Born-Oppenheimer molecular dynamics." It should be noted that the BLYP-TM results of Ref. [50] were obtained by combining PBE pseudopotentials and BLYP exchange-correlation functionals (private communication by the authors). We have used the following cutoff distances as determined from the first minimum of each g αβ (r): r GeGe = 3.1, 2.9, 3.1 Å, r GeSb = 3. q = 5/8 for fourfold coordination and q = 7/8 for threefold pyramidal ones [17].

q = 1 -3 8 i>k ( 1 3 + cos θ ij k )
Focusing on the BLYP-TM and PBE-TM cases only, the high fraction of fourfold Ge atoms is reflected by the corresponding peak in Fig. 6. This indicates that the majority of fourfold Ge sites has a tetrahedral configuration and the broadening up to q = 0.5 is due to defective octahedral fourfold Ge. The fraction of tetrahedral Ge can be estimated by integrating the q (n c = 4) distribution in the range q = 0.8-1.0. We obtain 68.7% for BLYP-TM and 65.7% for PBE-TM with the remaining fraction pertaining to defective octahedral configurations. The broad distribution in the range q = 0-0.6 for the fivefold Ge is also a signature of defective octahedral configuration.

IV. ELECTRONIC PROPERTIES: NETWORK CHEMICAL BONDING

We seek further insight into the interplay between atomic structure and chemical bonding of GST by resorting to the Wannier formalism, based on the notion of Wannier functions and centers WFC (the W label referring to these degrees of freedom). The WFC formalism provides a powerful tool for an unambiguous characterization of the bond network with respect to the traditional analysis based uniquely on g αβ (r) [53,54]. Recently, we have employed WFC centers to quantify the different behavior of PBE and BLYP in describing the fingerprints of the bonding and lone pairs subnetworks of glassy GeTe 4 [START_REF] Bouzid | [END_REF]. We note that, other methods such as the crystal orbital overlap population method could provide useful insight into the nature of chemical bonds as it was applied in the case of amorphous GeTe [55].

In Fig. 7, we compare the pair correlation functions g TeWFC (r), g GeWFC (r), and g SbWFC (r) obtained within the BLYP-TM and PBE-TM frameworks. All three pair correlation functions show two main peaks, corresponding, the first, to the WFC lone pairs [centered at 0.56 Å (Te-WFC), 0.29 Å (Ge-WFC), and 0.39 Å (Sb-WFC)] and the second, to the bonding WFC [centered at 1.19 Å (Te-WFC), 1.42 Å (Ge-WFC), and 1.59 Å (Sb-WFC)]. In the case of glassy GeTe 4 , BLYP and PBE differed by the presence or absence of a split second peak in g TeWFC (r), denoting homopolar and heteropolar bonding. It appears that such marks are broadened and merged in the case of GST. In view of the very similar BLYP and PBE coordination environments, the main difference between the two XC schemes in terms of correlations with the Wanniers centers is found in the intensity of the peak corresponding to the WFC lone pairs. This difference simply reflects dissimilar amounts of miscoordinated atoms, at the very origin of the presence of lone pairs. The corresponding units are visualized in the right panels of Fig. 7, according to the formalism proposed in Refs. [6,54]. In the definition of M(m,n), M is the central atom (Ge, Sb, Te), m is the nearest-coordinated neighbor, and n is the number of lone pairs. The above analysis show that the better performances of BLYP in terms of comparison with experimental structural data are not reflected by any striking difference in the electron localization behaviors.

V. CONCLUSIONS

In conclusion, our BLYP-TM model of glassy GST provides a very good account of experimental structural results. We are able to recover realistic values for Ge-Te bonds, corresponding to a combined tetrahedral/octahedral network. Within our model, the tetrahedra are clearly predominant, this conclusion being in agreement with several experimental pieces of evidence and in disagreement with others [START_REF] Carria | [END_REF]6,9,10,12,14,15,20]. For instance, it was argued in Ref. [20] that Raman features arising from Ge atoms in tetrahedral coordination are absent in the experimental counterpart, calling for Raman calculations to be performed on the present configurations (see Ref. [56]). Our PBE-TM model is also very well performing, albeit at a slightly lower level in terms of agreement with the pattern of the total structure factor and the interatomic distances. Overall, the bonding features of the two schemes are found to be quite similar, as shown via an analysis of the correlations between the ionic positions and the Wannier centers positions. The predictive power of the BLYP exchange-correlation functional for disordered chalcogenide materials stems from a specific choice of the corresponding Ge PP. It remains true that the debate on the atomic structure of glassy GST awaits further experimental and theoretical contributions, better accounting, as done recently [3,57,58], for structural relaxation and nucleation processes.

  FIG. 1. Calculated and experimental neutrons (a) and x-ray (b) total structure factor S(k) of amorphous Ge 2 Sb 2 Te 5 . X-ray and neutrons experimental data are taken from Refs. [26,27], respectively. The BLYP-TM (magenta line), PBE-TM (red line), BLYP-GTH (green line), and BLYP-TM-D (brown line) are compared to experimental data (dashed blue line).

FIG. 2 .

 2 FIG. 2. Calculated and measured total neutrons structure factor S(k) of amorphous Ge 2 Sb 2 Te 5 . (a) The BLYP-GTH results (green line), SGCPMD-PBE-GTH simulations from Refs. [19] (red line) and BOMD-PBE-TM simulations from Ref. [18] (maroon line) are compared to experimental data from Ref. [27] (dashed blue line). (b) The BLYP-TM results (magenta line) is compared to experimental data (dashed blue line).

  84%) in the BLYP-TM case [Fig.2(b)] deserves to be rationalized in detail. In this context, the selection of the exchange-correlation functional cannot be invoked to account for the recorded behaviors, since BLYP-GTH, BOMD-PBE-TM, and SGCPMD-PBE-GTH lead to very close results. Therefore we analyzed the specific pseudopotential constructions adopted in the present and in previous studies.In order to understand what makes our Ge TM PP different from other PPs, we looked at the Kleinman-Baylander (KB) construction used in each work. As shown in Table

FIG. 3 .

 3 FIG. 3. Partial pair correlation functions of amorphous Ge 2 Sb 2 Te 5 for the BLYP-TM (magenta line), PBE-TM (red line), and BLYP-GTH (green line) models. and coordination numbers, the BLYP-TM one reproduces accurately the r GeGe , r GeTe , and r SeTe distances.

FIG. 4 .

 4 FIG. 4. Partial pair correlation functions Ge-Ge, Ge-Sb, and Ge-Te as obtained from BLYP-TM (magenta line) and BLYP-TM-D (orange line) models.

  3, 3.2, 3.3 Å, r GeTe = 3.2, 3.3, 3.5 Å, r SbSb = 3.2, 3.3, 3.4 Å, r SbTe = 3.4, 3.4, 3.4 Å, and r TeTe = 3.1, 3.1, 3.2 Å, for BLYP-TM, PBE-TM, and BLYP-GTH, PBE-GTH PBE-TM BLYP-TM PBE-TM HSE a:[9], b:[10] BLYP-TM PBE-TM BLYP-GTH

FIG. 5 .

 5 FIG. 5. Fractions of n-fold Ge, Sb, and Te atoms (n = 1, 2, 3, 4, 5 or 6) for the BLYP-TM (magenta bar), PBE-TM (red bar), and BLYP-GTH (green bar) models.

FIG. 6 .

 6 FIG. 6. Distribution of the local order parameter around Ge atom for the BLYP-TM and PBE-TM models. Results are given for different coordination numbers (n c ).

FIG. 7 .

 7 FIG. 7. (Left): Partial pair correlation functions g TeWFC (r) (top), g GeWFC (r) (center), and g SbWFC (r) (bottom). Blue lines refer to BLYP-TM and red lines to PBE-TM data. (Right) Structural units in GST identified with the formalism proposed by Refs. [6,54] based on M(m,n), where M is the central atom (Ge, Sb, Te), m is the nearestcoordinated neighbor, and n is the number of lone pairs. WFC centers corresponding to bonds and lone pairs are depicted by light blue and green spheres, respectively.

TABLE I .

 I Maximum angular momentum l max used in the KB construction of pseudopotentials.

		Kleinman-Baylander construction
	Simulation	Ge	Sb	Te
	TM this work

ACKNOWLEDGMENT

We acknowledge Pôle HPC Equip@Meso of the University of Strasbourg and GENCI (DARI N. x2016095071) for computer time allocations.


Hayat Zaoui

Pier Luca Palla

Fabrizio Cleri

Évelyne Martin

Evelyne Lampin
email: evelyne.lampin@univ-lille1.fr

Fourier-like conduction and finite 1D thermal conductivity in long silicon nanowires by approach-to-equilibrium molecular dynamics

Fourier-like conduction and finite one-dimensional thermal conductivity in long silicon nanowires by approach-to-equilibrium molecular dynamics

I. INTRODUCTION

The thermal conductivity of silicon nanowires (NWs) 115 to 22 nm in diameters and micrometer long has been measured more than 10 years ago. 1 At room temperature, it decreases from about 150 W.K -1 .m -1 in bulk silicon 2 to about 10 W.K -1 .m -1 in the thinnest NWs. Such a drop in thermal conductivity has opened up a wide interest for thermoelectricity, a major improvement in the figure of merit of thermoelectric converters being expected provided that a high value of electrical conductivity can be preserved. Still a controversy remains whether the heat transport behaves or not according to the Fourier law. This controversy is fed by the length dependence of the thermal conductivity. Experimentally, it is extremely challenging to measure the spatial distribution of temperature and heat flow at the nanoscale. This challenge becomes ever dramatic when very thin NWs are targeted, such as the sub-10 nm NWs used in the most advanced nanoelectronic technologies. The problem has been however adressed on other 1D systems : carbon and boron nitride nanotubes 3 and SiGe nanowires 4 . The thermal conductivity presents a dependence even for lengths L longer than the average phonon mean free path (MFP) l ph that one can extracted from the bulk thermal conductivity κ bulk using the kinetic theory (C is the heat capacity and v the sound velocity):

κ bulk = 1 3 Cvl ph (1)
i. e. around 50 nm for Si at room temperature. A length dependence of Si NWs thermal conductivity is also obtained by atomic scale simulations 8-10 using the direct method 5 give a length dependence of the thermal conductivity of silicon nanowires [8][9][10] . At short length, a first regime shows a pronounced length dependence, which is attributed to a ballistic regime for L < l ph [8][9][10] . It must be emphasized that a length dependence is also obtained using the same method when bulk Si is studied 5 . At long length, the calculations by Hu 10 saturate to a constant thermal conductivity while N. Yang 8 and X. Yang 9 obtain a divergence of κ versus L.

The length dependence of the thermal conductivity is not directly related to a violation of the Fourier law nor to a ballistic regime. Indeed, the length dependence is also observed in atomistic simulations on bulk silicon. Recently, it has been shown 6,7 that the phonon MFPs are distributed on several orders of magnitude, and reach 10 µm in bulk silicon. The average value approximated by the above kinetic theory can not therefore be a criterium of saturation of the thermal conductivity. Until the system length has not reached the maximum phonon MFP, the contributions of the latter progressively accumulate and result in a length dependence. Though, several questions remain open related to the width of the phonon MFPs distribution in Si NWs and to a possible deviation from the diffuse regime described by the Fourier law. In the present work, we adress these issues by a systematic study of the thermal conductivity dependence on the Si NWs diameter and length. The latter have been pushed up thanks to a very efficient atomistic MD method for the simulation of thermal properties that we recently developed: 11,12 the Approach-to-Equilibrium MD (AEMD). The speed up is due to the exploitation of heat transients, rather than the stationary regime common to other MD-based approaches. A length dependence is observed, as in the direct method. Its origin has been clearly explained [START_REF] Zaoui | Length dependence of thermal conductivity by approach-to-equilibrium molecular dynamics[END_REF] and its correlation with the phonon MFP distribution has been demonstrated. The methodology, recalled in the next section, is used to answer the following questions : does the temperature profile and evolution comply to the heat equation (i.e. Fourier-like) ? Is the thermal conductivity finite or not at long length ? What is the evolution of the thermal conductivity as a function of the diameter ? And finally what can we learn on the phonon MFP distribution in the nanowire ?

II. ATOMISTIC METHOD

Smooth-surface NWs are obtained by cutting a cylinder of diameter D and length L from a silicon block with diamond-lattice atomic structure (Fig. 1). The circular section is preferred to the square one since it is closer to the real shape of Si NWs usually fabricated from either the top-down or bottom-up technologies 14,[START_REF] Ernst | Proc. IEEE Int. Electron Device Meeting[END_REF] . The interactions between silicon atoms are modeled using the Tersoff interatomic potential [START_REF] Tersoff | [END_REF] . The lattice parameter is taken equal to that of bulk crystalline Si for this potential scheme, a 0 =5.45 Å at 500 K. The axis of the nanowire is along the [001] direction (noted z). The diameter D is varied from 1 nm for the thinnest wires to 14 nm for the largest. Periodic boundary conditions are applied along the z-direction.

The Approach-to-Equilibrium Molecular Dynamics method is used to create and monitor a temperature transient, by looking at the decay time of a temperature difference ∆T to obtain the NW thermal conductivity. Two Nosé-Hoover thermostats are used to split the NW into a cold block at temperature T 1 = 400 K for 0 < z < L/2, and a hot block at temperature T 2 = 600 K for L/2 < z < L. Due to the presence of the surface, the thermalisation of the hot/cold blocks could be slower compared to the bulk case but we checked that the duration of the constant {NVT} MD simulation can be maintained within the range 100 ps to 1 ns without impacting the final result. Therefore we most often take 100 ps, to further minimize the computational cost. After the thermalization phase the two thermostats are switched off, and the system is left free to reach thermal equilibrium at constant {NVE}. During the approach-to-equilibrium, the temperature difference ∆T between the two blocks is followed up. It can be seen from the semi-log plot in Fig. 2 that ∆T tends to 0 with an exponential decay, i. e. the same behaviour already observed in bulk materials 11,12 . The final temperature of the nanowire is equal to 500 K, therefore quantum effects, which are not described by classical MD simulations, should not be relevant (Debye temperature of silicon = 645 K).

A fine determination of the decay time τ is made for each decay curve. Many exponential fits are performed with an incremental initial time to extract the mean value and error, and to remove the initial fast transient corresponding to the switching from {NV T } to {NV E}.

In this way, we recover only the intrinsic thermal decay time. In addition, averaging is performed over a couple of simulations for the smallest sizes more sensitive to the particular set of initial atomic velocities determined randomly following a Maxwell-Boltzmann distribution.

III. TRANSIENT ANALYSIS AND 1D THERMAL CONDUCTIVITY

The nanowires studied by molecular dynamics have a finite section by construction.

Therefore, the thermal current does not need to be normalized by unit section to be a finite quantity as in the case of a bulk or by unit length as in the case of a plane like graphene. The 1D formulation of the Fourier law then gives:

ϕ(W) = -κ 1D (W.m.K -1 ) gradT (K) ( 2
)
where ϕ is the thermal current and κ 1D is a 1D thermal conductivity.

We now derive the 1D heat equation in the nanowire. The nanowire supercell unit of length L with periodic conditions along z that we study by molecular dynamics is schema-tized in Fig. 3. The amount of heat δQ supplied to this unit between t and t + dt is equal to:

δQ = (ϕ i -ϕ o ) dt (3)
where ϕ i and ϕ o are the inward and outward thermal currents. The change in enthalpy dH of this unit between the same time interval is equal to:

dH = ∂ ∂t L 0 ρCT dz dt (4)
with C the heat capacity (J.K -1 ) and ρ the number 1D density in m -1 . As in the derivation of the heat equation in 3D, the time derivative can be displaced inside the integral where it applies to T only:

dH = L 0 ρC ∂T ∂t dz dt (5)
On the other hand, δQ can be rewritten using an integral over z :

δQ = - L 0 ∂ϕ ∂z dz dt = - L 0 κ 1D d 2 T dz 2 dz dt (6)
According to the first thermodynamics law, the two quantities dH and δQ are equal and we obtain the 1D heat equation:

∂T ∂t = κ 1D .L 3Nk B ∂ 2 T ∂z 2 (7)
where the Dulong and Petit law has been assumed for the heat capacity. This equation can be solved as in the bulk case 17 and the solution is a Fourier series with a dominant contribution :

T (t, z) = 1 2 (T 1 + T 2 ) + 2 π (T 2 -T 1 ) sin 2π L z × e -t/τ (8)
where

τ = 3Nk B L 4π 2 1 κ 1D (9)
The temperature difference during the approach-to-equilibrium is a sum of decaying exponentials:

∆T (t) = ∞ m=0 8(T 2 -T 1 ) (2m + 1) 2 π 2 e -(2m+1) 2 t/τ (10)
This expression tends to a single exponential decay for long times that corresponds to the major contribution for m = 0. The decay time extracted by the procedure presented in the previous part corresponds to this contribution, and can therefore be used to obtain the 1D thermal conductivity following:

κ 1D = 3Nk B L 4π 2 1 τ (11)
The NW 1D thermal conductivities are plotted versus length in Fig. 4, for diameters ranging from 1 to 14 nm. The speed up provided by the exploitation of the transient regime allowed us to study diameter and lengths that far exceed previous simulations [8][9][10] .

The number of atoms in the biggest systems (D=14 nm and L = 1.2 µm) exceeds 9 millions.

The thermal conductivities strongly depend on the NW length. A length dependence has also been observed when we apply the same methodology to bulk systems and was attributed to a progressive accumulation of the contributions of the largest phonon MFPs when the length of the simulation supercell increases. [START_REF] Zaoui | Length dependence of thermal conductivity by approach-to-equilibrium molecular dynamics[END_REF] The length dependence in the NWs and the consequences for the phonon MFPs will be discussed in the following sections.

IV. CONDUCTION REGIME IN THE NANOWIRES

In the range of lengths where the NWs thermal conductivity varies, an open question is whether the Fourier regime is valid, and therefore whether the temperature profile along z and the temporal evolution of the temperature difference between the hot and cold blocks still agree with Eqs. 8 and 10. We have investigated this point in the case of the thicker NW, because the temperature fluctuations will be the smallest among the familly of NWs studied due to an average on a larger number of atoms in a slice. We have looked into details to three supercell lengths : a long length (600 nm) in the plateau of the κ(L) curve, and two small lengths in the region of strong variation of κ : 75 and 150 nm. For these 3 lengths, the temperature difference behaves exponentially at long time as evidenced by the red lines superimposed on the MD evolution in Fig. 5. However, a small deviation from the exponential decay exists at the beginning of the approach-to-equilibrium, and can be seen in the figure for the longest wire before 50 ps. As discussed in the previous part, it is due to the switching from {NVT} to {NVE}. We have followed up the temperature profile in intervals far away from this initial transient, at times where only the Fourier regime should be accounted for. The intervals are indicated by arrows in Fig. 5 for each size. In these intervals, the temperature difference is around 100 K. The temperature profiles averaged on these intervals are presented in Fig. 6. The temperature profiles are presented on the same z scale, they are repeated according to the periodic boundary conditions for the two smaller lengths. The MD results in black are sine function of the form T (z) = T m + ∆T sin(2πz/L) as evidenced by the superimposing of the red curves obtained by adjusting only T m and ∆T .

The temperature profiles are therefore of the form obtained in Eq. 8 by solving the heat equation, even for the smallest length in the domain where the dependence of the thermal conductivity with the length is significant. In conclusion, both the temporal evolution of the temperature difference of the hot and cold blocks, and the temperature profile along the nanowire length do agree with the heat equation. We are therefore in a Fourier regime even if the thermal conductivity is not that of bulk silicon, and does vary with the supercell length.

V. THERMAL CONDUCTIVITY AT LONG LENGTH

A. Does the conductivity diverge ?

In the litterature, the thermal conductivity of the nanowires computed by MD is found to saturate at long length by Hu et al. 10 while it diverges in the studies by N. Yang et al 8 and X. Yang et al 9 . In these three works, the direct method is used to obtain the thermal conductivity of square-section nanowires. The same interatomic potential by Stillinger and Weber 27 is used, and the calculations are performed at room temperature, and up to 1000 K in the works of N. Yang and X. Yang. The thermal conductivities that we obtained by AEMD do not present any divergence when the NW length increases. In particular, Fig. 7 focuses on the D=2 nm nanowire studied in this work, one of our nanowire with the highest aspect ratio, up to 250. For lengths greater than 200-250 nm, the thermal conductivity remains constant. The error bar increases for the two longest NWs with L=400 and 500 nm due to a progressive self-heating at long time that prevents the whole temporal evolution of the temperature difference to be valuable for the determination of the decay time. This self-heating is reduced when the initial temperature difference is reduced for example to 100 K. 18 It could be due to a long-time surface relaxation. The nanowire studied in Fig. 7 has a cross section comparable to the NWs studied by N. Yang and X. Yang. X. Yang obtained a thermal conductivity at L=400 nm twofold the thermal conductivity at L=200 nm, while N.

Yang observed a twofold increase of the thermal conductivity from L=200 and L=1000 nm.

It has been checked that the temperature effects are not at the origin of such a difference, since our calculations also saturate at 300 and 1000 K for the same diameter. The origin of the divergence is unclear, and its very existence appears quite doubtful. It must be noticed however that is possible to extract finite thermal conductivities from experiments. 1

B. Comparison with the litterature

The length dependence is an effect of the present MD method, and since our calculations converge to a finite thermal conductivity at long length, the relevant value is the thermal conductivity that corresponds to a nanowire of a given diameter but infinite length. These values are presented in Fig. 8. As a matter of comparison with bulk thermal conductivities or bulk-like determination of nanowire thermal conductivities, it is possitble to obtain a pseudo 3D-thermal conductivity κ by dividing the 1D thermal conductivity by the nanowire section S:

κ 3D (W.m -1 .K -1 ) = κ 1D (W.m.K -1 ) S(m 2 ) ( 12
)
The variation range is smaller for κ 3D than for κ 1D since the smaller the section, the smaller the 1D thermal conductivity, as can be seen in Fig. 8 where the 3D thermal conductivities obtained from our 1D values are plotted in blue. It may explain, at least in part, the quasi absence of variation from D=2 nm to 15 nm in the EMD results. 28,29 At D=1 nm however, the increase of thermal conductivity is not obtained in our approach. In these NWs, the number of atoms by plane is only equal to 4, or in other words the core of the NWs is nearly nonexistent and the inconsistencies in the system are not relevant. The 3D thermal conductivities obtained from our 1D values range from κ 3D = 3.5 W.m -1 .K -1 for D=1nm and 28 W.m -1 .K -1 for D=14nm. The reduction compared to a silicon bulk described in the same conditions (AEMD, Tersoff potential, 500K, i. e. κ bulk =145±2 W.m -1 .K -1 ) is ≈40

for D=1nm and ≈5 for D=14nm. It is not possible to compare these values to experiments.

Indeed although typical NW diameters of 10-15 nm are used in NW-FETs, 14,[19][20][START_REF] Hiramoto | Proc. 10th IEEE ICSICT[END_REF][START_REF] Koyama | Proc. of the European solid-State Device Research Conference[END_REF] and that ultra-thin NWs with diameters in the [3-5] nm-range [START_REF] Goldberger | [END_REF]24 were successfully fabricated in the most advanced technologies, the thermal conductivity has been measured only for NW diameters greater than 22 nm. 1 Concerning the modeling, Zou and Balandin 25 solve the Boltzman transport equation by accounting for the confinement effects on the phonon dispersion curve.

They determine the dependence of the thermal conductivity on the nature of the reflection at the surface, from specular to diffuse. In molecular dynamics, both confinement and surface effects are implicitly accounted for. Our values are considerably smaller than their result in the specular case, that would a priori be more adapted to the smooth nanowires studied in the present work. However the maximum diameter we were able to study is smaller than the case they studied (D=20nm), and a decrease of the thermal conductivity is expected when the diameter decreases. An additional comparison element could be obtained thanks to the determination by molecular dynamics of the impact of the surface roughness, but is outside the frame of the present work. Concerning atomistic calculations, a large part of the results published in the litterature have been performed on square-section NWs. Differences are expected between circular and square sections due to the presence of edges. Moreover, it is not clear whether the results should be compared at a given cross section or at a given characteristic length (diameter or side). A strict comparison will therefore be restricted to calculations of NWs with a circular shape. The work of Ponomareva et al. 26 is dedicated to circular NWs oriented along the [111] direction. Although they use a different interatomic potential, the Stillinger-Weber 27 , known to give higher thermal conductivity for bulk Si 13 , their value of 17 W.m -1 .K -1 for D=4.2 nm gives a 1D thermal conductivity of 2×10 -14

W.m -1 .K -1 is in agreement with our results (Fig. 8). The thermal conductivity reduction of 3-4 obtained by He and Galli 29 by EMD and Tersoff potential in a cylindrical NW is obtained at room temperature. We have applied the ratio to the bulk thermal conductivity at 500 K and converted the value to a 1D thermal conductivity (κ 1D =7×10 -14 W.m -1 .K -1 ).

This point for D=15 nm is plotted in Fig. 8 and is in agreement with our calculations.

VI. PHONON MEAN FREE PATHS IN THE NANOWIRES

We have recently demonstrated [START_REF] Zaoui | Length dependence of thermal conductivity by approach-to-equilibrium molecular dynamics[END_REF] that the length dependence observed in the AEMD method is due to a cutoff at L of the phonon MFPs distribution. When increasing L, the thermal conductivity increases until a threshold length (L th ) is reached. Above L th , the thermal conductivity is constant and results from the whole distribution of phonons MFPs in the supercell unit. The value of L th corresponds to the upper limit of the phonon MFP distribution in the material Λ max . The curves κ 1D = f (L) plotted in Fig. 4 are used to obtain the value of Λ max as a function of the NW diameter (Fig. 9). The maximum phonon MFPs are in the range [100-800] nm for the diameters under study. These values are samller than in bulk silicon : using the same methodology, we were not able to reach the saturation of the κ(L) curve even for the largest supercell at L=1.2 µm. Other calculations in the litterature 6,7 show that the mean free path distribution extends up to 10 µm in bulk silicon.

The MFPs, although drastically lower than in bulk silicon, still remain large, hundred of times the NW diameter. This could mean that a way to improve the dissipation of selfheating in NW-based devices would be to preserve these long phonon MFPs for example thanks to a surface state as smooth as possible. In order to quantify the contribution of these long MFPs to the thermal conductivity, we present in Fig. 10 with the two latter diameters, could benefit of an improvement in the heat dissipation if the surface states is optimised, and even if the length of the NWs themselves will probably be smaller than the maximum phonon MFPs, they are generally connected to larger blocks acting for example as source/drain where these phonon large MFPs could exist.

VII. CONCLUSION

In conclusion, we systematically calculated the thermal conductivity of Si nanowires with diameters from 1 to 14 nm and lengths up to 1.2 µm by approach-to-equilibrium molecular dynamics. The temperature transient and spatial distribution comply with the timedependent solution of the heat equation, whenever the length is equal to a few tenths of nms or several hundreds of nms. Therefore a deviation from the Fourier regime is not evidenced.

A 1D thermal conductivity is obtained from the temperature transient. 26 In green the calculation from He et al. 29 The two latter, obtained in W.m -1 .K -1 by the authors, have been converted to 1D thermal conductivities using the cross section.

FiguresFIG. 1 :

 1 Figures

FIG. 2 :FIG. 3 :FIG. 4 :FIG. 5 :FIG. 6 :FIG. 7 :FIG. 8 :

 2345678 FIG. 2: Temperature difference betwen the hot and cold blocks during the approach-to-equilibrium. D=10 nm and L = 200 nm.

FIG. 9 :FIG. 10 :

 910 FIG.9: Maximum phonon mean free path versus NW diameter.

  an histogram of the contribution to the total conductivity of the first, second, third and last quarter of the MFPs distribution. For the 3 diameters, the distribution is approximately the same : 60 to 70 % of the thermal conductivity is provided by phonon MFPs distributed from 0 to Λ max /4. Another 25 % is provided by MFPs between Λ max /4 and 3.Λ max /4. The longest MFPs, from 3.Λ max /4 to Λ max contributes to ≈ 5% of the total thermal conductivity. For the thinnest NWs (D = 3nm), 30 % of the thermal conductivity is provided by phonon MFPs in the range [75:300] nm. This range shifts to larger values when the diameter increases : [150:500] nm for D=10 nm and [200:700] nm for D=14 nm. In conclusion, the NW-FETs acutally developed

  The length dependence of the thermal conductivity is analysed. At small length, the thermal conductivity increases when the length increases. At long length, the thermal conductivity saturates. The threshold corresponds to the maximum value of the phonon MFPs distribution. It remains as large as hundreds of nanometers in the smooth NWs studied in this work. The low-length regime corresponds to the accumulation of the phonon MFPs contributions, as observed in bulk silicon. At long length, all the MFPs are included and the thermal conductivity is fixed at a given diameter. A divergence of the thermal conductivity is not observed. dynamics, in Amorphisation and thermal properties of semiconductors: From bulk to nanostructures, ed. K. Termentzidis, Pan Stanford (in press).

VIII. ACKNOWLEDGEMENTS

This work was funded by the French ANR via the project n. ANR-13-NANO-0009 "NOO-DLES" and by Equipex EXCELSIOR (http://www.excelsior-ncc.eu). The computations were performed by using resources from GENCI (Grand Equipement National de Calcul Intensif) (Grant No.[2016]097600, [2015]096939, [2014]096939). Evelyne Lampin acknowledges Christophe Delerue for discussions on the 1D electical and thermal conductivities.