M Danoix

Helena Gouné

Philippe Zapolsky

S Chentouf
email: sara.chentouf@gmail.com

S Cazottes

F Danoix

M Goune

H Zapolsky

P Maugis

Effect of interstitial carbon distribution and nickel substitution on the tetragonality of martensite: a first-principles study

come

Effect of interstitial carbon distribution and nickel substitution on the tetragonality of martensite: A first-principles study

S. Chentouf, Sophie Cazottes, F.

Introduction

The relation between the lattice parameters and the carbon content of ferrous martensite has been the subject of many experimental studies in the last decades. From the review of experimental measurements, Cheng et al. [1] have established the linear dependences of the bct iron-carbon martensite a and c parameters as functions of the number of interstitials per iron atom, 𝑥 # : 𝑎 = 0.28664 -0.027 ± 0.001 𝑥 # 𝑛𝑚 Eq.1 𝑐 = 0.28664 + 0.243 ± 0.003 𝑥 # 𝑛𝑚 Eq.2

From these equations, the tetragonality ratio is:

𝑐/𝑎 = 1 + 0.95 ± 0.01 𝑥 # Eq.3 where 0.28664 nm is the adopted lattice parameter of pure bcc-iron at 300 K.

The above equations match well many experimental data. However, some differences in tetragonality are observed for steels having the same carbon content [2]. It was shown that these differences are related to both interstitial ordering [3][4] and alloying [2,[5][6].

It is well accepted that tetragonality of martensite formed by quenching from high-temperature cubic face centered austenite phase is a direct consequence of carbon atoms preferential location at one of the three possible sub-lattices of octahedral interstices. In 1948, Zener [7] has used the elastic interactions to explain the temperature-dependent transition between ordered and disordered interstitial carbon distributions in martensite. Later, Khachaturyan [START_REF] Khachaturyan | Theory of structural transformation in solids[END_REF] proposed a detailed description of expected strain induced ordering of solutes using the Microscopic Elasticity Theory (TEM). He demonstrated that the strongly repulsive chemical interactions and elastic interactions between carbon atoms must be taken into account in order to predict the order-disorder transition. The calculations of Khachaturyan have been recently reviewed by Udyanski [START_REF] Udyansky | [END_REF] et al. using an ab initio EAM potential for Fe-C system and molecular dynamic simulations to estimate the chemical interactions between carbon atoms.

It was also shown that martensite tetragonality depends on substitutional content. It has been observed that freshly formed martensite steels containing Mn, low-Ni, Rh and Cr have an abnormally low tetragonality [10,[11][12]. On heating these steels from room temperature the c parameter increases and the a parameter decreases [13][14][15]. Meanwhile, abnormally high tetragonality was found in freshly formed Al and high-Ni martensite [2,5]. Heating freshly formed Al and high-Ni martensite led to a decrease of c and increase of a parameters. More recently, by using ab initio calculations Al-Zoubi et al. [16] found that Al and Ni addition enhances significantly the c/a ratio of ordered Fe 16 C 1 .

The present study was carried out to investigate (i) the effect of interstitial carbon atom content and distribution and (ii) the effect of the substitutional Ni atoms on martensite tetragonality. The variation of the a and c lattice parameters is investigated at different carbon contents, ranging from small to high concentrations. The typical carbon content in supersaturated ferrous martensite is about 1 wt.% C (i.e. 4.48 at.%). However, a martensite quenched from austenite without tempering, if held at temperatures close to room temperature, will decompose to form a modulated structure of carbon-rich and carbon-lean areas. TEM and one-dimensional tomography revealed that in Fe-Ni-C alloys the nanoscale carbon-rich modulated structure results from spinodal decomposition [17][START_REF] Olson | Martensite. A Tribute to Moriss Cohen[END_REF]. Then, the carbon content, the structure and the stability of the carbon-rich regions have been the subjects of numerous experimental studies. It was suggested that these zones consist of an ordered phase. Because of the high carbon concentration of these regions (about 11.1 at.%), it is assumed that carbon atoms ordering in one type of octahedral site might be important. Thus, three types of long-range order distributions can be considered, from the least to the most ordered: (i) carbon randomly distributed on the three octahedral sub-lattices. That is the case of highly supersaturated ferrite noted a-Fe 16 C 2, (ii) carbon ordered on one octahedral sub-lattice and randomly distributed on its octahedral sites. That is the case of Zener ordered martensite noted a Z -Fe 16 C 2 and (iii) carbon ordered onto one octahedral sub-lattice with occupation of only some of its sites leading to the stoichiometric structure a²-Fe 16 C 2 . Taylor et al. [17] have found that, for carbon compositions close to 11 at.% and for low temperatures, the stoichiometric structure type (iii) has the lowest free energy among the three types. The partially ordered a Z -Fe 16 C 2 structure (type (ii)) would be the next most stable while the least stable would correspond to a-Fe 16 C 2 (type (i)). These results have been confirmed by Molecular Dynamic (MD) calculations performed by Sinclair et al [START_REF] Sinclair | [END_REF].

In the first part of this paper, the effect of interstitial carbon content (up to x c =12.5 at.%) and ordering on the structural parameters of martensite is investigated. Besides, to get information about the spinodal decomposition, the mixing energies are calculated for different carbon concentrations by taking as references the stoichiometric structure a²-Fe 16 C 2 and pure bcc-iron. The Special Quasirandom Structure method has been applied to generate different possible carbon distributions on the octahedral sites within the bcc-iron lattice. Special attention has been paid to ordering of the structure a²-Fe 16 C 2 . To discuss the differences between the mixing energies of supercells with different distributions of carbon, the carbon-carbon interaction energies are calculated. In the second part, the effect of 25 at.% of Ni substitution on the lattice parameters of martensite is studied. The calculations of equilibrium lattice parameters and energies have been carried out using the Density Functional Theory (DFT).

II. Calculation details

The first-principles calculations are performed using the Projector Augmented Wave (PAW) [20] method as is implemented in the VASP code [21][22]. The calculations employ the Generalized Gradient Approximation (GGA) [23] for exchange-correlation energy. The calculations are done with the energy cutoff of 400 eV. k-points meshes are created by Monkhorst-Pack scheme [24]. 4x4x4 k-meshes are used for all cells with 128 atoms. The energy convergence criterion of the electronic self-consistency is chosen as 10 -5 eV for all calculations. The reciprocal space energy integration is performed by the Methfessel-Paxton [25] technique for structure relaxation. All calculations are performed within the spin polarized approximation.

III. Results and discussion

III.1. Ordering of carbon atoms in Fe 16 C 2

To study the properties of Fe 16 C 2 compound, the three types of long-range ordering are considered, corresponding to the following structures: (i) a-Fe 16 C 2 structure where carbon atoms are randomly distributed on the three octahedral sub-lattices, (ii) a Z -Fe 16 C 2 structure where carbon atoms are randomly distributed on one octahedral sub-lattice and (iii) a²-Fe 16 C 2 structure where carbon atoms are ordered and occupy specific sites on one octahedral sublattice. The unit-cells of these three structures are shown in Fig. 1.

First, the total energies per iron atom for these three structures have been calculated. For the two structures a Z -Fe 16 C 2 and a-Fe 16 C 2 , three SQS supercells (with 128 Fe atoms + 16 C atoms) have been considered with different carbon distributions on their respective interstitial sublattices. The results of these calculations are summarised in Table 1. Our results are in good agreement with those obtained from molecular dynamics calculations [START_REF] Sinclair | [END_REF]. The ordered a²-Fe 16 C 2 structure is the most stable one among the three considered structures. The partially ordered a Z -Fe 16 C 2 structure is the second most stable while the fully disordered a-Fe 16 C 2 is the least stable. The ordering energy from a to a² is 24.9 meV/Fe atom.

The relaxed a²-Fe 16 C 2 structure is perfectly tetragonal with a ratio c/a=1.107. The calculated lattice parameters using the PAW [20] potential and PBE [23] approximation for the exchange correlations function reproduce well the previous theoretical data obtained using the USPP-GGA methods [START_REF] Sinclair | [END_REF]. For the case of a Z -Fe 16 C 2 and a-Fe 16 C 2 structures, the average c/a ratios are 1.10±0.01 and 1.03±0.02, respectively. The c/a uncertainty is defined as the standard deviation of the three SQS supercells. Since, in a-Fe 16 C 2 structure the carbon atoms are randomly distributed in the three octahedral sub-lattices, the calculated tetragonality, as expected, is close to 1. On the other hand, the carbon distributions on one octahedral sub-lattice lead to a c/a ratio of around 1.1 for both a Z -Fe 16 C 2 and a²-Fe 16 C 2 structures.

III.2. Carbon-carbon interactions in martensite

The carbon-carbon interactions have been calculated in order to discuss the differences between the mixing energies of supercells having different distributions of carbon atoms for a given concentration. In this section, the carbon-carbon interactions from the first nearest neighbours (pair length a Fe /2) to the nineteen nearest neighbours (pair length 2a Fe ) in the bcc-Fe are presented and discussed.

For designation, the first 19 positions of carbon octahedral sites as function of pair length are shown in Figure 2. The carbon-carbon interaction energies calculated using supercells containing 128 atoms of iron and two atoms of carbon, are summarised in Table 2. Fig. 3 presents the comparison of our results with the theoretical results obtained by Domain et al. [26]. The interaction energy is defined as:

𝑉 𝐴 : , 𝐴 < = 𝐸 𝐴 : + 𝐴 < + 𝐸 >?@ -𝐸 𝐴 : + 𝐸 𝐴 < Eq. 4

where 𝐸 𝐴 : + 𝐴 < , 𝐸 >?@ and 𝐸 𝐴 A are the energies obtained when only atomic positions have been allowed to relax. 𝐸 𝐴 : + 𝐴 < is the total energy of a supercell containing two carbon atoms. 𝐸 >?@ and 𝐸 𝐴 A are the total energies of the supercells without carbon and with one carbon atom in the octahedral sub-lattice, respectively. With this definition, a positive interaction energy means repulsion while a negative interaction energy denotes attraction. All supercells have the same size, i.e. equivalent to the size of the relaxed supercell of bcc-iron. Note that 𝐸 𝐴 : = 𝐸 𝐴 < for bcc lattice.

Our calculations (Table 2) show that most interactions are repulsive. They predict that the interactions become slightly attractive for configuration 19 (400) while the previous DFT results of Domain et al. [26] predict attractive interactions for configurations 6 (102), 8 (220) and 15 (222). Nevertheless, our results seem to be in agreement with the calculations of Ruban [27] and Hatcher et al. [28]. The DFT results by Ruban et al. [27] and Hatcher et al. [28] can be found in graphical form in the original papers. We suggest that the discrepancies in interaction energies are directly related to the size of cell used for calculations and the details of the electronic calculations. On the other hand, from Figure 3, it can be seen that the interactions are large at the first shells and decrease with distance with an exception of shells 4 (020) and 18 (040) when one or two iron atoms are located along [010] direction between the two carbon atoms. Ruban [27] has found that strain-induced interactions become dominant just at the point where carbon atoms start attracting each other, i.e. from about the 14 (311) interaction. However, in our calculations we did not separate the chemical and elastic interactions. (USPP-GGA)

V(1-2) 1 1/2 ½ [001] 1.9537 V(1-2) 2 2/2 ½ [011] 0.8298 0.65 V(1-1) 3 3/2 ½ [111] 0.2669 0.09 V(1-1) 4 1 ½ [020] 1.964 1.67 V(1-1) 5 1 ½ [200] 0.015 0.09 V(1-2) 6 5/2 ½ [102] 0.046 -0.13 V(1-2) 7 6/2 ½ [211] 0.073 V(1-1) 8 2 ½ [220] 0.0379 -0.14 V(1-1) 9 2 ½ [202] 0.1557 V(1-2) 3/2 ½ [300] 0.0645 V(1-2) 3/2 ½ [122] 0.1722 V(1-2) 10/2 ½ [031] 0.0712 V(1-1) 11/2 ½ [131] 0.3035 V(1-1) 11/2 ½ [311] 0.0154 V(1-1) 3 ½ [222] 0.0306 -0.16 V(1-2) 13/2 ½ [302] 0.0015 V(1-2) 14/2 ½ [132] 0.0055 V(1-1) 2 ½ [040] 0.8286 V(1-1) 2

III.3. Effect of carbon concentration on the lattice parameter of Fe 16 C 2X

To study the structural properties of ferrous martensite, the carbon atoms are distributed randomly on a single octahedral sub-lattice with occupation of specific sites only: precisely those sites that correspond to carbon sites in the a²-Fe 16 C 2 stoichiometric structure with 12.5 at.% of carbon. Here the martensite is approximated by the phase Fe 16 C 2X . The parameter X (0<X<1) is the occupancy factor of the carbon sublattice, i.e. the ratio of occupied interstitial sites by carbon atoms to the total number of interstitial sites in the interstitial sublattice of structure a². This gives us the possibility to treat the coherent spinodal decomposition between Fe (X=0) and a²-Fe 16 C 2 (X=1). The calculations have been performed with different supercells containing 128 iron atoms. Nine concentrations (0£x c £12.5 at.%) are considered in total with 0, 2, 4, 6, 8, 10, 12, 14 and 16 interstitial atoms of carbon. For each concentration, three different supercells are used and an average of the structural parameters is calculated.

As described above, the a²-Fe 16 C 2 structure (with 16 carbon atoms and 128 iron atoms) is perfectly tetragonal. However, the structures with lower concentrations of carbon are slightly orthorhombic with very close lattice parameters a and b. In these cases, the average (a+b)/2 is calculated to represent the lattice parameter, noted a FeC , for the considered carbon concentration.

In order to scale our results in the same way as the data of Cheng et al. [1] and the molecular dynamics results of Becquart et al. [29], the c FeC /a Fe and a FeC /a Fe ratios are presented in Figure 4. Our calculations have been done at 0 K. In this case the lattice parameter of iron is a Fe =2.833 Å. It should be noted that both the experimental data of Cheng et al. and the MD results of Becquart et al. [29], given for comparison, are obtained at room temperature (300K). At this temperature, the lattice parameters of pure iron are 2.866 Å and 2.855 Å, respectively. These values are slightly higher than the one obtained in the present work. It is well known that the ab initio calculations at 0 K underestimate the lattice parameter of bcc-iron. It can be seen in Fig. 4 that our calculations remarkably reproduce the linear dependence of lattice parameters with carbon content up to the high concentration of x C =12.5 at% corresponding to Fe 16 C 2 . Quantitatively, our results are in good agreement with the c FeC /a Fe ratio given by Cheng et al. [1] from the experimental data up to 8 at.%. This indicates that our calculations reproduce the expansion of the lattice parameter c of martensite as a function of carbon concentration. However, for the case of a FeC /a Fe , from about 2.5 at.% of carbon our results diverge slightly, by 1%, from both the experimental data and theoretical values of Becquart et al. [29].

III.4. Mixing energies of Fe 16 C 2X

The mixing energies of different structures of Fe 16 C 2X system as function of carbon concentration are calculated with supercells containing 128 atoms of iron. The mixing energy per iron atom, as reference to pure iron and a²-Fe 16 C 2 stoichiometric structure, is defined as: Where, 𝐸 𝐹𝑒 :<G 𝐶 I is the total energy of a cell containing 128 atoms of iron and n atoms of carbon (n=16X). 𝐸 𝐹𝑒 :<G is the calculated total energy of bcc-iron and 𝐸 𝐹𝑒 :<G 𝐶 :J is the total energy of the ordered a²-Fe 16 C 2 structure.

∆𝐸 CAD =
The calculated mixing energies at nine different concentrations (0£x c £12.5 at.%) are presented in Figure 5. Similarly to the lattice parameters, for each concentration, calculations of ∆𝐸 CAD have been performed for three different SQS supercells with different carbon spatial distributions on the interstitial sub-lattice. For each concentration of carbon, the average of the three different mixing energies is also calculated and presented in Figure 5. Remembering that X = 8x C is the occupancy factor of the a'' interstitial sub-lattice (0<X<1), the average energy has been fitted according to the sub-regular solution model [36] in the form of a polynomial: Lattice parameters (lattice units)

Carbon concentration Xc (at.%) aFeC[Present work] cFeC [Present work]

aFe C [Becquart] cFeC [Becquart] aFeC [Cheng] cFeC [Cheng] ∆𝐸 CAD = 𝜔 N?# 𝑋 1 -𝑋 < Eq.6

with w FeC = 43 meV. As can be seen in Figure 5, the mixing energies are positive for all carbon concentrations. The curve of the average ∆𝐸 CAD is not symmetric around the medium concentration x C =6.25 at.% (i.e with 8 interstitial carbon atoms): it reaches a maximum for a carbon concentration around 4.7 at. % (i.e. with 6 interstitial carbon atoms).

Among the three different mixing energies for each concentration, some are grouped very close together and others are scattered. These differences are likely to be related to the carbon distributions in the interstitial sub-lattice and therefore to the carbon-carbon interactions. To clarify this point, we have chosen to analyse the carbon atom distribution in the three different SQS supercells for the carbon concentration of 4.68 at. %. For this concentration, the differences between the mixing energies are significant (from 3.13 to 9.86 meV/Fe atom).

We have identified all carbon-carbon pairs in the three different SQS supercells. Table 3 summarises the number of different carbon-carbon pairs as a function of their lengths. The carbon-carbon interaction energies have been calculated for the first nearest (1NN) and second nearest neighbour (2NN) pairs in the a²-Fe 16 C 2 structure. The details of the calculations and the results of the interaction energies are presented and discussed in section III.2. In the case of the a²-Fe 16 C 2 structure, the 1NN and the 2NN correspond to the configurations 15 and 19 with lengths 3 a Fe and 2a Fe respectively (see Figure 2).

The two values of the interaction energies for 1NN and 2NN carbon atoms are +0.025 and -0.028 eV, respectively. Therefore, the interactions are repulsive between the first nearest neighbours and attractive between 2NN. Supercells '3' with the low mixing energy has eight 1NN pairs, while the two other supercells have nine 1NN pairs each. Thus, supercell '3' possesses the least number of repulsive interactions and this is why its mixing energy is lower Here, supercell '2' has additional attractive interaction energies between the 2NN than the supercell '1'. Consequently, the mixing energy of supercell '2' is lower than that of structure '1'. Note that, the number of 1NN and 2NN pairs is only one part of the explanation of the energy differences between the three SQS structures.

III.5 Effect of carbon concentration on the lattice parameters of (Fe,Ni) 16 C 2X

In the following sections, we investigate the role of Ni substitution on the properties of Fe 16 C 2X .

As in the case of Fe 16 C 2X , nine carbon concentrations are considered for the calculations of the lattice parameters of (Fe,Ni) 16 C 2X . Calculations have been performed using disordered structures by considering simultaneously (i) the substitutional disorder of nickel atoms on the iron sub-lattice and (ii) the interstitial disorder of carbon atoms on the interstitial sub-lattice.

Here also, the positions of interstitial carbon sites correspond to those of the a²-Fe 16 C 2 structure.

In order to scale the lattice parameters of (Fe,Ni) 16 C 2X systems to those of Fe 16 C 2X , we have presented the c FeNiC /a Fe and a FeNiC /a Fe ratios in Figure 6. Note that the resulting lattice parameter of carbon-free SQS-Fe25at.%Ni structure for one bcc unit cell is a FeNi = 2.854 Å [30]. This lattice parameter is calculated using the same methods in the present study.

Figure 6 shows that the lattice parameters of (Fe,Ni) 16 C 2X are linear functions of the carbon concentration. Additionally, the expansion and contraction of the lattice parameters are greater in (Fe,Ni) 16 C 2X than in Fe 16 C 2X . For example, at x C = 12.5% the longitudinal deformations 𝜀 = 𝑐 -𝑎 N? 𝑎 N? induced by carbon insertion are e FeNiC = 0.145 and e FeC = 0.106. We show in the following that this behaviour can be explained in the frame of the elasticity theory of point defects. Indeed, the observed linear dependence of the lattice parameters vs the carbon concentration can be written 𝑎 𝑥 Q = 𝑎 R 1 + 𝛿 T 𝑥 # and 𝑐 𝑥 Q = 𝑎 R 1 + 𝛿 Q 𝑥 # where d a and d c are the lattice expansion concentration coefficients. The lattice parameter a 0 for carbon-free alloys being respectively a Fe and a FeNi for Ni-free and 25at.%Ni systems (see Table 4). For lattice parameter c, 25 at.% addition of nickel results in an increase of carbon-induced deformation up to 37%: e FeNiC / e FeC = 1.37. This finding is in qualitative agreement with the DFT results of Al-Zoubi et al. [16]: the authors have found that 5 at.% of nickel enhances significantly the c/a ratio of Fe 16 C 1 . According to Bacon [32], Cochardt et al. [33] and Clouet et al. [34], the quadratic deformation induced by one carbon atom on a z-octahedral site in bcc lattice can be calculated from the dipole tensor P ij associated to the point defect. According to the tetragonal symmetry of the octahedral interstitial sites, the dipole tensor is diagonal and is characterized by only two components: P a and P c related respectively to the lateral and longitudinal deformations induced by the carbon atom. The components of the dipole tensor verify the following relations:

𝑃 T = : < 𝑎 R
V 𝐶 :: + 𝐶 :< 𝛿 T + 𝐶 :< 𝛿 Q Eq.7

𝑃 Q = : < 𝑎 R V 2𝐶 :< 𝛿 T + 𝐶 :: 𝛿 Q Eq.8
where C ij are the elastic constants of the carbon-free matrix and d a , c are the lattice expansion concentration coefficients defined above. Applying these equations to bcc-Fe yields the values P a = 8.03 eV and P c = 15.4 eV. On the other hand, for SQS structure Fe25at.%Ni, P a = 4.48 eV and P c = 9.65 eV. It is clear that both dipole components and elasticity constants are weaker in the case of Ni substitution. A rough estimate of the consequences on the longitudinal deformation can be derived as follows. From the inversion of Eq. 7 and 8 we obtain:

𝛿 Q = < T W X Y<# Z[ \ ] ^ (# ZZ ^#Z[ )\ _ (# ZZ Y# Z[ )(# ZZ ^<# Z[ )
Eq.9

Noting that the component 𝑃 T has a minor influence on d c compared to 𝑃 Q , it appears that d c scales approximately as the ratio P c /C¢, where C¢=(C 11 -C 12 )/2 is the shear modulus. From Table 4, we have (P c /C¢) FeNiC /(P c /C¢) FeC =1.39, which is close to the value of 1.37 observed on e. We conclude that the higher c/a ratios in the nickel-bearing system is due the higher compliance of the latter due to nickel substitution.

Based on the results of neutron diffraction experiments [35], Kurdjumov et al. [31] have deduced that during the martensitic transformation of FeNiC system with high Ni content, all carbon atoms occupy the sites of a single octahedral sub-lattice and consequently induce the maximum tetragonality for a given carbon content. While, for the case of other freshly formed martensite i.e with low Ni content and Ni-free, some of carbon atoms (about 20 %) are distributed on the two other octahedral sub-lattices. In our calculations of Fe 16 C 2X and (Fe,Ni) 16 C 2X systems , 100% of the carbon atoms have been located on the z-octahedral lattice. Nonetheless, a high effect of nickel substitution on tetragonality is observed. We conclude that the influence of the nickel-induced softening on tetragonality cannot be discarded in this respect. where, 𝐸 𝐹𝑒 `J𝑁𝑖 V< 𝐶 I is the total energy of a cell containing 96 atoms of iron, 32 substitutional atoms of nickel and n interstitial atoms of carbon, and 𝐸 𝐹𝑒 `J𝑁𝑖 V< is the total energy of the SQS-Fe25at.%Ni structure. 𝐸 𝐹𝑒 `J𝑁𝑖 V< 𝐶 :J is the total energy average of two structures with disorder on the substitutional lattice only, the interstitial carbon sublattice being totally filled.

The calculated mixing energies for nine different concentrations are presented in Figure 7.

Similarly to the case of lattice parameters, for each concentration, two different SQS-(Fe,Ni) 16 C 2X supercells, with different substitutional nickel and interstitial carbon distributions, have been taken into account for the calculations of ∆𝐸 CAD . Additionally, for each carbon concentration, the average of the two different mixing energies is calculated and presented in Figure 7. The average energy of mixing is fitted with an equation of type Eq. 6 with w FeNiC = 75 meV.

As shown in Figure 7, all average mixing energies are positive. Similarly to the Fe 16 C 2X system, the curve of the average ∆𝐸 CAD is not symmetric around the medium concentration x C =6. 25 In order to quantify the influence of 25at.%Ni substitution on the driving force for unmixing, the free energy of mixing DF mix = DE mix -TDS mix is calculated. Under the hypothesis of random distribution of the carbon atoms on the interstitial sites of a''-Fe 16 C 2 , the entropy of mixing per metal atom writes ∆𝑆 CAD = -: G 𝑘 e 𝑋𝑙𝑛𝑋 + 1 -𝑋 𝑙𝑛 1 -𝑋 , Eq.11

where X is the occupancy factor of the interstitial sublattice, as defined above. The free energies of mixing for both systems are presented in Figure 8. The temperature T=300K has been chosen for the purpose of illustration. Free energy of mixing (meV/metal atom)

Carbon concentration Xc (at.%)

FeC

FeNiC

The vanishing slopes of DE mix functions in the vicinity of the stoichiometric compound Fe 16 C 2 (i.e. for 𝑥 # close to 12.5 at%) (Figures 6 and7) result in minimal and negative values of DF mix in this region. From the common tangent construction, a two-phase equilibrium appears at T=300K between a very dilute phase (𝑥 # ≈ 0) and an under-stoichiometric phase (𝑥 # <12.5%), for both systems. In addition, from the shape of the curves, it can be concluded that the driving force for unmixing of martensite is higher with 25%Ni alloying than without.

IV. Conclusion

Using first-principles calculations based on the PAW potential and PBE approximation for the exchange correlation, the effect of carbon distribution and nickel substitution on the (Fe,Ni) 16 C 2X system has been investigated. Three structures have been studied: a²-Fe 16 C 2, a Z -Fe 16 C 2 and a-Fe 16 C 2 . Our results show that the stoichiometric ordered a² structure is the most stable among the three. The lattice parameters of both Fe 16 C 2X and (Fe,Ni) 16 C 2X systems vary linearly with carbon concentration, up to the high value of 𝑥 # =12.5 at.%. We find that nickel substitution enhances the expansion of parameter c of the Fe 16 C 2X system by 37%. This is related to the fact that the disordered Fe25at.%Ni alloy is more compliant than pure bcc-iron. The mixing energies of Fe 16 C 2X and the (Fe,Ni) 16 C 2X systems are positive for all carbon concentrations. The highest mixing energy value for (Fe,Ni) 16 C 2X (~11 meV/metal atom) is greater than that of Fe 16 C 2X (~6 meV/Fe atom). This indicates that unmixing of martensite by spinodal decomposition is favoured by nickel substitution. It also suggests that the higher Zener ordering observed in nickel-bearing iron alloys is due to the softening effect of nickel substitution.

Figure 1 .

 1 Figure 1. (a) The a-Fe 16 C 2 unit cell with the positions of carbon sites on the three (x, y and z) octahedral sub-lattices. (b) The a Z -Fe 16 C 2 unit cell. The carbon atoms occupy the z-octahedral interstitial sub-lattice, with a lattice parameter close to a Fe . (c) The a²-Fe 16 C 2 unit cell with the positions of carbon sites on the z-octahedral interstitial sub-lattice. The lattice parameter is close to 2a Fe . Iron atoms are in grey, carbon sites are coloured in accordance to their sublattice (blue, green and red for x, y and z respectively).

Figure 2 .

 2 Figure 2. The first nineteen positions of the carbon octahedral interstitial sites. Iron atoms in grey, carbon sites coloured in accordance to their sub-lattice.

Figure 3 .

 3 Figure 3. Calculated total carbon-carbon interaction energies in (eV). The ab initio results of Domain et al. [26] are shown by red spheres.

Figure 4 .

 4 Figure 4. Lattice parameters of Fe 16 C 2X system (lattice units) as a function of carbon concentration. For comparison, the results of Cheng et al. [1] and Becquart et al. [29] are presented. The dotted lines correspond to the linear interpolations.

Figure 5 .

 5 Figure 5. Mixing energy (in meV/Fe atom) of Fe 16 C 2X structures as a function of carbon concentration.

Figure 6 .

 6 Figure 6. Lattice parameters of the (Fe,Ni) 16 C 2X and Fe 16 C 2X systems (in units of a Fe ) as functions of carbon concentration.III.6. Mixing energies of (Fe,Ni) 16 C 2XThe mixing energies of different (Fe,Ni) 16 C 2X structures as function of carbon concentration have been calculated with supercells containing 128 metallic atoms (with 96 Fe+32 substitutional Ni). In these calculations, the mixing energy per atom was calculated using the following expression:

Table 1 .

 1 The

	Energy (meV/Fe atom)	-24.6	n/a	-9.9	0
	Lattice parameter, a (Å)	5.64	5.71	5.64±0.02 5.77±0.05
	Lattice parameter, c (Å)	6.25	6.31	6.22±0.04 5.94±0.08
	Tetragonality, c/a	1.107	1.105	1.10±0.01 1.03±0.02

equilibrium total energies per iron atom (meV/Fe atom), by taking as reference the energy of a-Fe 16 C 2, and lattice parameters of the a²-Fe 16 C 2 , a Z -Fe 16 C 2 and a-Fe 16 C 2 structures from the ab initio calculations. The lattice parameters correspond to the unit cells constructed with 2x2x2 bcc iron. Property a²-Fe 16 C 2 a²-Fe 16 C 2 [19] a Z -Fe 16 C 2 a-Fe 16 C 2

Table 2 .

 2 The calculated carbon-carbon interaction energies (in eV) at the same V(1-1) and different V(1-2) octahedral sub-lattices. The vectors of the nineteen neighbours as a function of pair length are correspondingly presented.

	Type of interaction	Distance of	Vectors This work-128 atoms	Domain et al. [26]
		separation	(PAW-PBE)
		(units of a Fe )

  other supercells. Let's compare now the number of second nearest neighbours of supercells '1' and '2'. Although they have the same number of 1NN, the first supercell has six 2NN and the second has eight.

	than for the
		14.00
	Mixing energy (meV/metal atom)	2.00 4.00 6.00 8.00 10.00 12.00						Model SQS Average
		0.00
		0.00	2.00	4.00	6.00	8.00	10.00	12.00	14.00
				Carbon concentration Xc (at.%)

Table 3 .

 3 Number and length of carbon-carbon pairs in the three SQS supercells of Fe 128 C 6 (i.e with carbon concentration x c =4.68 at.%). The mixing energies are also given to identify the different supercells.

	Length	Supercell '1'	Supercell '2'	Supercell '3'
			Mixing energy (meV/Fe atom)
		+ 9.86	+ 6.26	+ 3.13
			Number of C-C pairs in Fe 128 C 6
	3 a Fe	9	9	8
	2 a Fe	6	8	6

Table 4 .

 4 Parameters used in the elasticity model of point defects. The elastic constants C ij 's and shear modulus C¢ (in GPa) are extracted from[30].

		a 0 (Å)	d a	d c	C 11	C 12	C¢
	bcc-Fe	2.833	-0.025	0.84	267	147	60
	SQS-Fe25%Ni	2.854	-0.25	1.07	186	132	27

  8. Free energies of mixing (in meV/metal atom) for Fe 16 C 2X and (Fe,Ni) 16 C 2X systems as a function of carbon concentration.

		14.00
	Mixing energy (meV/metal atom)	0.00 2.00 4.00 6.00 8.00 10.00 12.00	0.00	2.00	4.00	6.00	8.00	10.00		12.00 Model Average SQS	14.00
		-2.00
					Carbon concentration Xc (at.%)
		10
		8
		6
		4
		2
		0
		0	2	4	6	8	10	12	14
		-2

Acknowledgments

This study has been performed in the framework of the ANR-13-BS08-0014 SPIDERMAN project funded by the Agence Nationale de la Recherche (ANR). The calculations were carried out partially on the mesocentre of the Aix-Marseille University (project No. 16b045), and partially on the clusters of IDRIS (project No. 100524).


Ege Dundar

Nicolas Chanut

Pascal Boulet

Philip L Llewellyn

Bogdan Kuchta

Modeling of adsorption of CO 2 in the deformed pores of MIL-53(Al)

Carbon dioxide is a greenhouse gas implicated in global warming [1] but also used in existing hydrogen production processes [2]. Therefore, with a future world economy based on hydrogen as an energy vector, coupled with recent commitments to the Kyoto Agreement whereby efforts must be made to reduce CO 2 emissions, there is a need for materials to efficiently capture and store CO 2 [3]. The MIL series of hybrid porous frameworks are candidates for the purpose of CO 2 adsorption. In particular, the breathing MIL-53 structures are of special interest because their adsorption properties can be modulated by applying external pressure to modify pore geometry [3, 4].

Molecular simulations were performed to predict CO 2 adsorption in flexible metal-organic frameworks. For this purpose, a generic force field was fitted to our experimental data to describe the non-bonded (electrostatic and Van der Waals) interactions between the CO 2 molecules and the large pore (lp) and narrow pore (np) forms of the MIL-53(Al) framework.

With this validated force field, it will be possible to predict CO 2 uptake and enthalpy of adsorption at various applied external pressures that will modify the structure's pore configuration and allow us to have more control over the adsorption/desorption process . This will be done in future work by creating rigid intermediate structure forms between the lp and np shapes, and simulating adsorption of CO 2 in each of them.

In parallel of our modeling strategy, we have devised an experimental device in which it is possible to modulate the porosity of the MIL-53 material with the aid of external mechanical pressure. Equally, calorimetry experiments have been carried out on the unconstrained material to directly measure the adsorption enthalpies. This set of experiments allows a comparison of the isotherms on the material constrained to various mechanical pressures, and the enthalpies on the non-constrained material, with the results obtained by molecular simulation. This joint experimental-computational approach allows for us to further understand this flexible class of materials and is a way forward to understanding the possibility to optimize pore size for specific gas separations.

1.0

Introduction

With recent global commitments to reduce greenhouse emissions, there is an increasing need for high quality materials capable of capturing and storing harmful gases like CO 2 1 .

Therefore, great attention has lately been given to the development and elaboration of porous materials such as metal-organic frameworks (MOFs). These are versatile, chemically and thermally robust materials containing reactive sites and large pores with high surface areas 2,3 .

Because of these important features, MOFs are able to efficiently capture and release guest molecules and thus find use in various industrial applications like gas storage/separation, catalysis, and biomedicine 3 . Some MOFs have been known to show selective flexibility during gas adsorption via unusual breathing 4 and gate opening process 5 properties. The MIL-53 (Cr, Al) (Matériaux de l'Institut Lavoisier) structures synthesized by the group of Férey are among the most interesting structures displaying these features because of their high uptake capacity for gases like CO 2 and CH 4 6 . Concerning CO 2 particularly, an unusual microscopic adsorption behavior in the aluminum and chromium-containing MIL-53 structures has been reported, based on microcalorimetric measurements. This behavior was associated with a breathing mechanism of the framework where the structure is able to reversibly change its configuration from a narrow pore (np) form to a large pore (lp) form 6 . Although the signature of these structural changes can be detected experimentally, the microscopic mechanism involved in the transitions is too fast to be well understood.

To date, computational methods have played a significant role in elucidating the mechanism of adsorption at an atomistic level and in providing insight on the microscopic arrangement of the adsorbates in the flexible framework pores. Ramsahye et al 1 The goal of this work was to verify whether a validated force field found in the literature can effectively predict our own experimental adsorption isotherms and enthalpy of adsorption measurements for the same gas, same adsorbent, in the same thermodynamic conditions. In the eventuality that it can, our goal is to create rigid intermediate structures that would produce isotherms with uptakes located between those of the lp and np forms to simulate various external pressures applied to the material. However, if the predictions are not satisfactory, then the goal is to fine-tune some of the force field parameters to improve the capture of our experimental data's behavior.

Computational methodology

Charge distribution

The initial atomic coordinates for the large pore (lp) and narrow pore (np) forms of the MIL-53(Al) framework were taken directly from the structure's crystallographic information files (cif) obtained by X-ray diffraction. Table I shows the crystallographic cell parameters for the structure and the references from which they were obtained. The obtained charges are presented in Table II. The calculations were performed using the PW91 GGA density functional 12

Interatomic interactions

The total energy (E) is expressed as the sum of the interaction energies between the adsorbate and the framework (E AF ) and that between the adsorbate molecules (E AA ) 14 :

𝐸 = 𝐸 !" + 𝐸 !! (1)
E AF and E AA are sums of pairwise additive potentials of the form:

𝑒 !" = 4𝜀 !" ! !" ! !" !" - ! !" ! !" ! + ! ! ! ! ! !" (2)
where e ij is the sum of the non-bonded energetic interactions for E AF and E AA . The first term in Eq (2) is the usual repulsion-dispersion Lennard-Jones (LJ) potential with ε ij and σ ij representing the potential well and the distance at which the potential between particles i and j is zero respectively, and r ij is the distance between the particles. The second term is the Coulombic contribution term accounting for electrostatic interactions between the particles, with q i and q j being the point charges for particles i and j respectively. The Lorentz-Berthelot mixing rules were applied in order to calculate the parameters for each pair of atoms:

𝜎 !" = ! ! !! ! ! (3)
𝜀 !" = 𝜀 ! 𝜀 ! (4)
The LJ potential term in Eq (2) can also be expressed in the following form, more appropriate for calculations:

𝑒(𝐿𝐽) !" = 𝐷0 !" !! !" ! !" !" - !! !" ! !" ! (5)
where R0 ij = σ ij *2 1/6 represents the distance at which the potential well is at its maximum absolute value of depth and D0 ij is the potential well in units of kCal mol -1 . Table III presents the LJ interaction parameters used in this work. Parameters for the CO 2 -MOF interactions were calculated with Eqs (3) and ( 4) from the LJ parameters for each atom type used by

Hamon et al 15 . In their work, the parameters describing the atoms in the organic part of the structure were taken from the Dreiding force field 16 and those describing the atoms in the inorganic part come from the Universal Force Field (UFF) 17 . Parameters for the CO 2 -CO 2 interactions were again calculated with the Lorentz-Berthelot mixing rules with LJ parameters for each atom taken from Yung and Harris 13 . To our knowledge, this force field is the most accurate for the simulation of CO 2 adsorption on the MIL-53 family of frameworks and was frequently and successfully used to explain the breathing behavior of these flexible frameworks 1,2,10,15 .

Grand Canonical Monte Carlo (GCMC) simulations

Absolute adsorption isotherms were computed using the GCMC algorithm implemented in Accelrys' Sorption module (ref Accelrys), which allows exchange, displacement (translation and rotation), regrowth and conformer moves 18 . These simulations consist of evaluating, for a given volume and temperature, the average number of adsorbate molecules whose chemical potential is equal to that of the bulk fluid. All simulations in this work were run at 303 K, using one unit cell for the np and lp forms of the MIL-53(Al) framework. Characteristically, 1

x 10 7 Monte Carlo (MC) steps were used for each run, the frameworks of both structures being considered rigid. The electrostatic interactions were taken into account by using the Ewald summation and the short range interactions were limited to a cut-off distance of 12 Å 19 .

Results and discussion

Fig. 2 shows the GCMC simulations for CO 2 uptake and enthalpy of adsorption in the MIL-53(Al) framework at 303 K. The LJ parameter values used for these simulations are taken from tables III and IV. In Fig. 2 a, although the equilibration values for the isotherms seem to be in the right range of values (8 and 2 mmol g -1 for the lp and np forms respectively) the gasframework interaction appears to be too strong, resulting in a steep uptake slope for both forms at low fugacity. As stated by Ramsahye et al 1 and mentioned above, between 0.25 and 6 bar the np form is the most energetically preferred form for CO 2 adsorption on MIL-53(Al).

Our experimental values show that for a fugacity of 0.25 bar, the CO 2 uptake is 0.5 mmol g -1 whereas the simulations show an almost fourfold adsorbed value of 1.89 mmol g -1 for the np form (Fig. 2 c). Also, for both forms of simulated isotherms, saturation seems to occur at low fugacity values while experimentally CO 2 molecules are still being adsorbed in the framework's pores at fugacities well over 10 bar. demonstrates the gap between the two. In the light of these results, we conclude that in our case the initial parameter values found in the literature do not produce an adequate prediction of CO 2 adsorption in the MIL-53(Al) framework. The reasons for this are not quite clear, as the same software and same simulation conditions as in 1,2,10,15 were followed and the structures were taken from the same references 8,9 . Differences between the experimental data found in the literature and those obtained in this work can however be attributed to density variations between batches, quantities used and accumulated experimental error with increasing pressure. Since it appears from Fig. 2 b that the enthalpy of adsorption for both lp and np forms of the structure is too high, we have tried to manually alter the R0 ij and D0 ij parameters for the CO 2 -CO 2 interactions, without touching the CO 2 -MOF interaction parameters. Furthermore, we have chosen not to modify the electrostatic interactions since the charges for the structure were obtained from our own DFT calculations and Mulliken charge partitioning method, and not taken from the literature. Electrostatic charges for the c_co2, o_co2 atoms were left unmodified as well. This strategy is rather unconventional but it was our thinking that the various modifications of the R0 ij and D0 ij parameters for the CO 2 -CO 2 interactions would give some valuable information about the extent of their sole influence on the adsorption isotherms and the enthalpy of adsorption. In order to get a sense of the dependence of the isotherms and the enthalpy of adsorption on the interaction between the CO 2 molecules, we first varied the D0 ij parameters by adding or subtracting 80% of their value while keeping R0 ij constant. Then we kept the D0 ij parameters at their original values and varied the R0 ij parameters by 20% and 80% of their value. The parameter values for variations in D0 ij ,and R0 ij are shown in Tables V and VI respectively, while Figures 3 and4 show variations in CO 2 uptake and enthalpy of adsorption due to changes in D0 ij ,and R0 ij respectively.

Influence of D0

Figure 3a shows us that the variation of the energetic D0 parameter produces what we would intuitively expect: by reducing this parameter value, the interactions between the CO 2 molecules decrease and result in a lower uptake than what was obtained with the original D0 value (Table IV). This behavior is much more evident for the LP form than for the NP form due to the steric effect of the latter. Indeed, in the NP configuration there is less space for the molecules to move around so the impact of reduction of the energetic interaction between adsorbate molecules is less striking. By increasing the D0 parameter the expected increase in adsorbed molecules occurs as well as a more pronounced slope at low pressure, suggesting high adsorbate-adsorbent interaction. This is probably due to the higher number of CO 2 molecules being clustered together, taking up more space and yielding higher enthalpy of adsorption. As was the case for the decrease in D0, the effect of the increase of this parameter is not obvious for the NP form of the structure because of the lack of space. These characteristics of CO 2 uptake variations due to changes in the D0 parameter can also be seen on Figure 3b where the enthalpy of adsorption is plotted as a function of fugacity. It is noticed that for the three parameter values of the LP form, the enthalpy of adsorption spans a much larger energy spectrum than that for the NP form. For the NP form, the three parameter values do not have a big influence on the enthalpy of adsorption, meaning that whether we increase or decrease D0 by 80% of its value or we leave it as it is, the energy released by the adsorbate due to adsorption on the pore walls remains arguably the same.

Influence of R0

Figure 4 illustrates the influence of the R0 parameter on CO 2 uptake and enthalpy of adsorption. It can be seen from Figure 4a that as the equilibrium distance parameter R0 decreases, uptake is increased because of the reduction of the repulsion and increase in the dispersion contributions of the LJ potential, thus allowing CO 2 molecules to come closer together. Consequently, increasing R0 increases the repulsion and decreases the dispersion contributions of the LJ potential and keeps clusters of CO 2 molecules from forming, lowering the amount of adsorbed molecules in the pore. This behavior is detected for both LP and NP forms, the latter being less pronounced than the former, although decreasing R0 by 80% of its value seems to have a more significant impact on the uptake than increasing it. An interesting phenomenon is observed when we look at the enthalpy of adsorption shown on Figure 4b: for both forms, the energy released due to adsorption is higher for both decreased and increased parameter values. Also, in both cases the enthalpy of increased R0 is smaller than the decreased one. Thus, for a case where the energy between the adsorbate molecules is too high, a decrease in both R0 ij and D0 ij parameters contributes in obtaining a better fit to the experimental data.

After many modifications of the c_co2 -c_co2 and o_co2 -o_co2 interactions (the third interaction parameters c_co2 -o_co2 are calculated using the Lorentz-Berthelot mixing rule once values for the first two are set), we obtained the optimal parameters shown in Table IX.

The isotherms and enthalpies of adsorption with the modified CO 2 -CO 2 parameters are shown in Fig. 5.

Conclusions

A generic force field with parameters widely used in literature was implemented to predict experimental adsorption data for CO 2 on MIL-53(Al) framework's open (lp) and closed (np) forms. However, the accuracy of the predictions was not satisfying with regard to our experimental isotherms and the calorimetrically measured enthalpy of adsorption. The main problems in the predictions with the original parameters were the large increase in enthalpy of adsorption with increasing pressure, and failure to capture the isotherms' shapes. Therefore, attempts were undertaken to study the susceptibility of the calculated properties to parameters' variations, and to fine-tune the LJ potential parameters, equilibrium distance R0, and energetic parameter D0, for the CO 2 molecule's carbon and oxygen atoms. First, a sensitivity analysis regarding the CO 2 LJ parameters D0 and R0 was carried out. Four variations of 10% of the original R0 and D0 parameter values were tested, while keeping every other parameter in the force field constant: 10% increases in D0 and R0, and 10% decreases in D0 and R0. This simple analysis allowed us to quantify to what extent, and in which "direction," a 10% variation in the parameter values would modify the CO 2 uptake and the energy generated by the adsorbed molecules. It followed from the analysis that a slight decrease in each parameter (5% and 7% in D0 and R0, respectively) improved the force field's prediction for the CO 2 experimental data.

It is important to note that our exploration of the parameter variations' impact is purely numerical. The low variation from the original parameters (10%) tested in this work was part of the goal of keeping the changes physically sound. However, despite this precaution, the parameters that best predicted the experimental adsorption data in our work (5% decrease in D0 and 7% decrease in R0) were not able to reproduce the CO 2 VLE curve accurately. It was found that variations of 5% in the D0 parameter, and around 2% in the R0 parameter should not be exceeded. Manual iterations in order to obtain an optimal combination of D0 and R0 parameter values where the adsorption experimental data would be better predicted while keeping the nature of the CO 2 molecule physically sound are possible but can prove laborious.

The discussed results confirm our hypothesis that adsorption of CO 2 in MIL-53 MOF structure is very sensitive to structural parameters. Moreover, unless the desired changes from the employed force field's original parameters are relatively small, manual fine-tuning of LJ parameters should be done with extreme caution, and, very importantly, validated for physical soundness.

One possible factor that could be important for the quality of the force field is that of multibody effects. Indeed their influence was not verified for CO 2 as the initial force field does not take them into account in the pairwise parameterization process. Therefore, it is in fact very likely that including a many-body contribution would explain some of the discrepancies, and result in both an appropriate VLE curve and a better uptake/enthalpy of adsorption prediction.

As suggested in recent work, taking into account many-body interactions can result in sharp changes in the grand-canonical partition functions of single-component systems and nanoconfined fluids, considerably deviating adsorption properties like isotherms and heats of adsorption [24, 25]. Although not in the present paper because of the necessity to modify the whole CO 2 model, which will require a non-negligible amount of time, it is definitely something that we should look into in the very near future, and is an excellent path to follow to answer some of the questions that have arisen in this work.

  have focused on the theoretical study of CO 2 adsorption in MIL-53 (Al) using Grand Canonical Monte Carlo (GCMC) simulations based on atomistic potential parameters and a partial charge model. By comparing the simulated isotherms and differential enthalpies of adsorption with experimental data they were able to observe that at pressures below 0.1 bar and above 6 bar, the framework takes the lp form while between 0.25 bar and 6 bar the np form is present. In their case the framework was considered rigid with the constituent atoms positioned at the coordinates reported by crystallography. In order to properly model the breathing of the MIL-53(Cr) structure, Salles et al 2 extended their simulations to encompass molecular dynamics. The force field thus had to be modified and included forces like bond potential, bending potential and torsion potential, in addition to the non-bonded Lennard-Jones and electrostatic interactions. By doing this the group successfully captured the two-step structural switching, and correctly reproduced the empty large and narrow-pore structures as well as the large magnitude of the breathing in presence of CO 2 . More recently, Garcia-Pérez et al have reproduced adsorption and simulated diffusion of CO 2 and CH 4 in the flexible NH 2 -MIL-53(Al) while using a rigid force field and a minimum number of experiments7 . Their twostep process consisted of an optimization of the interactions between the CO 2 molecules and the bridging hydroxyl groups of the framework, followed by the use of a linear combination of the np and lp forms based on experimental results.

Figure 1 .

 1 Figure 1. Labels of the force field types corresponding to each atom of the MIL-53(Al) structure and the CO 2 molecule.

Fig. 2 b

 2 shows the simulated enthalpy of adsorption for both forms of the MIL-53(Al) structure. It can be seen that although the np form's values display a relatively flat behavior and are in the higher than expected 40 kJ mol -1 range 1,19 , the lp form's increase from around 20 kJ mol -1 to 35 kJ mol -1 denotes uncharacteristically strong CO 2 -CO 2 interactions and excessively high energy values. Fig. 2 d compares our experimental data for enthalpy of adsorption to the simulated values and clearly

Figure 2 .

 2 Figure 2. GCMC simulations of (a) CO 2 adsorption isotherms (b) enthalpy of adsorption on MIL-53(Al) at 303 K using CO 2 -CO 2 interaction parameters from Ref 13. Figures (c) and (d) compare the simulated combinations of lp and np forms to experimental data.

Figure 4 .

 4 Figure 4. GCMC simulations of (a) CO 2 adsorption isotherms and (b) enthalpy of adsorption on MIL-53(Al) at 303 K using variations of R0 ij parameters for CO 2 -CO 2 interactions.

Figure 3 .

 3 Figure 3. GCMC simulations of (a) CO 2 adsorption isotherms and (b) enthalpy of adsorption on MIL-53(Al) at 303 K using variations of LJ parameters for CO 2 -CO 2 interactions.

Figure 5 .

 5 Figure 5. GCMC simulations of (a) CO 2 adsorption isotherms (b) enthalpy of adsorption on MIL-53(Al) at 303 K using our most accurate manually fitted CO 2 -CO 2 interaction parameters. (c) and (d) compare the simulated combinations of lp and np forms to our experimental data.

Figure 6 .

 6 Figure 6. Comparison between low pressure GCMC simulations of CO 2 adsorption isotherms and enthalpy of adsorption on lp and np forms of the MIL-53(Al) structure at 303 K. (a) and (b) are the initial CO 2 -CO 2 interaction parameters found in the literature and (c)and (d) are the manually fitted parameters we have used in this work.

Table I :

 I Crystallographic

	c	12.813	6.805
	α	90.000	90.000
	β	90.000	105.85
	γ	90.000	90.000
	Space group	Imma (no 74)	C 2/c (no 15)
	Since the positions of the H atoms are not detectable by experimental techniques, they were
	added to the framework's organic group using Accelrys Materials Studio 8.0's Visualizer
	software [ref Accelrys] 10 . The µ 2 positions' H atoms were added by editing the structures' cif
	files. The periodic models were then geometry optimized with Accelrys' Dmol 3 module [ref
	Accelrys], using the crystallographic coordinates of the atoms as starting configurations. With
	the same module, the partial charges for the porous structures were then extracted using
	density functional theory (DFT) calculations and the Mulliken charge partitioning method 11
	data for MIL-53(Al) structure used in this paper's calculations
	Crystallographic	MIL-53(Al)
	data
		large pore 8	narrow pore 9
	a	6.6085	19.716
	b	16.675	8.310

Table II :

 II Atomic partial charges carried by the different atoms in the open and closed forms of MIL-53(Al), and in the CO 2 molecule in electron units. The framework and CO 2 molecule labels are described in Figure 1.

	Atoms	h_c	cg1	cg2	c_c	o_c	µ2_o	h_o	Al
	MIL-53(Al) lp	0.139	-0.07 -0.083 0.584 -0.557 -0.712 0.299	1.392
	MIL-53(Al) np	0.134 -0.071 -0.077 0.566 -0.562 -0.741 0.292	1.478
			Atoms		c_co2	o_co2
			CO 2		0.6512 13	-0.3256 13

Table III :

 III Repulsion-dispersion Lennard-Jones parameters for each interacting pair of the adsorbentadsorbate interaction. The parameters were taken from Ref.15.

	Adsorbent-adsorbate pair D0 ij , kCal mol -1	R0 ij , Å
	c_c -c_co2	0.073056	3.5004
	c_c -o_co2	0.12374	3.6553
	cg1 -c_co2	0.073056	3.5004
	cg1 -o_co2	0.12374	3.6553
	cg2 -c_co2	0.073056	3.5004
	cg2 -o_co2	0.12374	3.6553
	h_c -o_co2	0.049286	3.2849
	h_o -c_co2	0	1.5473
	h_o -o_co2	0	1.7022
	µ 2 _o -c_co2	0.057756	3.2984
	µ 2 _o -o_co2	0.097822	3.4533
	o_c -c_co2	0.057756	3.2984
	o_c -o_co2	0.097822	3.4533

Table IV :

 IV Repulsion-dispersion Lennard-Jones parameters for each interacting pair of the adsorbateadsorbate interaction. The parameters were taken from Ref.13.

	o_co2 -o_co2	0.16002	3.4044
	Adsorbate-adsorbate pair	D0 ij , kCal mol -1	R0 ij , Å
	c_co2 -c_co2	0.055782	3.0946
	c_co2 -o_co2	0.094478	3.2495

Table V :

 V Variation of the repulsion-dispersion Lennard-Jones parameter D0 ij for each interacting pair of the adsorbate-adsorbate interaction. GCMC simulations of (a) CO 2 adsorption isotherms and (b) enthalpy of adsorption on MIL-53(Al) at 303 K using variations of D0 ij parameters for CO 2 -CO 2 interactions.

	Figure 3.
	Adsorbate-adsorbate pair	-80 % * D0 ij	+80 % * D0 ij
	c_co2 -c_co2	0.011156	0.10041
	c_co2 -o_co2	0.018896	0.17006
	o_co2 -o_co2	0.032004	0.28803

Table VI :

 VI Variation of the repulsion-dispersion Lennard-Jones parameter R0 ij for each interacting pair of the adsorbate-adsorbate interaction.

	Adsorbate-adsorbate pair	-80 % * R0 ij	-20 % * R0 ij	+20 % * R0 ij +80 % * R0 ij
	c_co2 -c_co2	0.61893	2.4757	3.7135	5.5703
	c_co2 -o_co2	0.64991	2.5996	3.8994	5.8491
	o_co2 -o_co2	0.68088	2.7235	4.0853	6.1279

Table VII :

 VII Effects of the variation of the repulsion-dispersion Lennard-Jones parameters D0 ij on the

	CO 2 uptake and enthalpy of adsorption
	Variable	-80 % * R0 ij	+80 % * R0 ij
	CO 2 loading, mmol g -1	0.61893	5.5703
	Q ads , kJ mol -1	0.64991	5.8491

Table VIII :

 VIII Effects of the variation of the repulsion-dispersion Lennard-Jones parameters R0 ij on the

	CO 2 uptake and enthalpy of adsorption
	Variable	-80 % * R0 ij	-20 % * R0 ij	+20 % * R0 ij +80 % * R0 ij
	CO 2 loading, mmol g -1	0.61893	2.4757	3.7135	5.5703
	Q ads , kJ mol -1	0.64991	2.5996	3.8994	5.8491

Table IX :

 IX Manually modified repulsion-dispersion Lennard-Jones parameters for each interacting pair of the adsorbate-adsorbate interaction.

	Adsorbate-adsorbate	D0 ij , kCal mol -1	R0 ij ,Å
	pair	(this work)	(this work)
	c_co2 -c_co2	0.01	2.80
	c_co2 -o_co2	0.01732	2.90
	o_co2 -o_co2	0.03	3.00


H Balout

P Boulet

M.-C Record

Strain-induced electronic band convergence: effect on the Seebeck coefficient of Mg 2 Si for thermoelectric applications

Keywords: Thermoelectricity, isotropic constraints, silicides, band structures, Density-functional theory, Boltzmann transport theory

The present theoretical study, performed using density-functional theory and Boltzmann transport theory formalisms, shows that, under 2.246% isotropic tensile strain, the two energy-lowest conduction bands of Mg 2 Si overlap. The two, threefold degenerated orbitals become a unique, sixfold degenerated orbital. It is demonstrated that, such degeneracy implies an increase of the Seebeck coefficient, of the electrical conductivity, of the power factor, and in fine of the figure of merit.

Introduction

The rarefaction of fossil fuels and their environmental impact have made the concerns about energy grown fast over the past decades. New ways of energy production, such as batteries, photovoltaic and thermoelectricity, are intensively explored. The principles of thermoelectricity are based on the direct conversion between electrical and thermal energies [START_REF]Thermoelectrics Handbook: From Macro to Nano[END_REF]. The performance of a thermoelectric material is measured by three parameters, namely the Seebeck coefficient S, the electrical conductivity σ and the thermal conductivity κ; these parameters are connected by a dimensionless factor, the so-called figure of merit zT , with zT = T S 2 σ /κ, where T is the absolute temperature [START_REF] Ioffe | Semiconductor Thermoelements, and Thermoelectric Cooling[END_REF]. According to this equation, it appears that a way to reach high zT is that the material shows a large Seebeck coefficient, a good electrical conductivity and, at the same time, a low thermal conductivity. However the improvement of the heat-to-electricity conversion by materials is made difficult due to the inverse relationships between the electronic thermal conductivity and the electrical conductivity (Wiedemann-Franz law) and between the Seebeck coefficient and the electrical conductivity (Boltzmann's transport theory).

Hence, the best materials found so far result from the best compromise between the transport properties.

The improvement of the power factor S 2 σ , and hence of zT , could be obtained by separately controling the electrical conductivity and the Seebeck coefficient, which can be achieved by band structure engineering. Band engineering consists in electronic band manipulations so as to create density of states (DOS) resonant levels [START_REF] Heremans | [END_REF] or electronic band degeneracy (see for instance [4][5][6][7]). The DOS effective mass is then increased which in turn increases the Seebeck coefficient of the material [START_REF] Heremans | [END_REF]8,9]. Chasmar and Stratton [10] showed that, when assuming that in doped materials the dopants only affect the chemical potential level, the best figure of merit value can be achieved by increasing the β quality factor defined as:

β = hk 2 B 3π 2 N v C L m * I Ξ 2 κ L T (1)
With h and k B the Planck and Boltzmann constants, respectively, and N v , C L , m * I , Ξ , and κ L the number of degenerate valleys for the band, the average longitudinal acoustic modulus, the inertial effective mass, the de-Strain-induced electronic band convergence: effect on the Seebeck coefficient of Mg 2 Si for thermoelectric applications 3 formation energy potential coefficient, and the lattice thermal conductivity, respectively. Hence, the strategy is to maximize the band degeneracy while minimizing the effective mass and the deformation potential coefficient if acoustic phonon scattering dominates, which is known to be the case in the best, known thermoelectric material. In previous theoretical studies dealing with both compressive and tensile isotropic strains applied to Mg 2 Si the authors have shown that for Mg 2 Si convergence of the two energy-lowest conduction valleys occurs by applying an isotropic tensile strain on the crystal structure [11,12]. Since the energy difference between these two valleys is small (∆ E = 173 meV) we can expect that for a reasonable amplitude of the strain they will become degenerated. In Mg 2 Si these orbitals are already threefold degenerated, hence the valley under constrain will be sixfold degenerated. This should enhance the Seebeck coefficient. The purpose of the present work is to demonstrate this feature. This investigation has been performed by computational approach. The paper is divided into two sections in which the calculations procedure and the results are presented, and a conclusion summarizes this work.

Calculations Procedure

Magnesium silicide Mg 2 Si crystallizes in a face-centered cubic lattice with a lattice parameter a equals to 6.35 Å [START_REF] Villars | Pearson's Crystal Data: Crystal Structure Database for Inorganic Compounds[END_REF]. This arrangement possesses the anti-fluorite structure. In the primitive cell, the silicon atom occupies the 4a (0, 0, 0) site and the magnesium atoms occupy the 8c (0.25, 0.25, 0.25) sites. The Fm3m space group fixes the fractional coordinates of all atoms. Each Mg is surrounded by a tetrahedron of Si atoms, and each Si is surrounded by a cube of Mg atoms.

The primitive cell has been adopted to simulate isotropic strains that are modelled by changing the lattice parameter a 0 . Throughout the paper the isotropic strain is given in units of the relative change of the lattice parameter as ∆ a/a 0 = a/a 0 -1, where a 0 is the relaxed cell parameter. Since only a tensile strain tends to make the two sets of low-lying orbitals of the conduction band converge [12], in the present work ∆ a/a 0 is positive.

All the structural optimizations have been performed by means of first-principle density-functional theory (DFT) calculations [START_REF] Hohenberg | [END_REF]15] as implemented in the Quantum Espresso code [16]. The self-consistent total energy calculations have been carried out using the plane-wave and ultrasoft pseudopotential approaches [17].

The generalized gradient approximation (GGA) Perdew-Wang91 functional has been used for the electronic exchange-correlation interactions treatment [START_REF] Perdew | Electronic Structure of Solids '91[END_REF]. A cutoff energy of 410 eV has been used throughout this Ry/Bohr.

From these well converged results, the Seebeck coefficient has been obtained by solving the linearized Boltzmann equation within the constant relaxation time approximation as implemented in the BoltzTraP code [START_REF] Madsen | [END_REF]. It is customary to assume a constant relaxation time (τ) value for evaluating transport properties since they involve the derivative of the Fermi function, which is significant only in a small energy window around the Fermi level and τ usually does not change much over this interval. In Addition, it allows for the direct calculation of the Seebeck coefficient as a function of temperature and doping level with no adjustable parameters.

Results

Determination of the optimal tensile strain

In order to determine the optimized cell parameter a, the Mg 2 Si structure has been fully relaxed. The optimized cell parameter a equals 6.37 Å, which shows a reasonably agreement with the experimental value (6.35 Å [START_REF] Villars | Pearson's Crystal Data: Crystal Structure Database for Inorganic Compounds[END_REF]).

Figure 1 shows the electronic band structures of the bulk Mg 2 Si. The path is plotted along the high symmetry k-points W -Γ -X -W . The Fermi level is marked by a dotted line and placed at E = 0 eV. Mg 2 Si is a semiconductor, its electronic band structure presents an indirect band gap, where the valence band maximum (VBM) and the conduction band minimum (CBM) are located at the Γ and X k-points, respectively. The band gap energy along the Γ -X line amounts to about 256 meV, which is underestimated by about 500 meV with respect to experimental data [20] but in agreement with previous theoretical calculations (e.g., [21,22]). The band structure shows also that the direct band gap is located at Γ with an energy of about 1830 meV. In addition, the two energy-lowest conduction orbitals present a threefold degeneracy and are separated by an energy of 173 meV at the X k-point. These two sets of orbitals tends to converge towards each other, as described in the next paragraph, when the compound undergoes tensile strain. Figure 2 presents the evolution of the energy of the two lowest conduction orbitals at the X k-point with respect to the isotropic tensile strain. It can be seen that the isotropic tensile strain leads to a band crossing resulting in a sixfold degenerated conduction orbital for the ∆ a/a 0 values of 2.246 %.

In the next sections the calculation of the properties of the constrained material are performed using this value of the isotropic tensile strain.

3.2 Band structure and transport properties under constraint

Band structure

Figure 3 shows the electronic band structure of strained Mg 2 Si with ∆ a/a 0 = 2.246 %. Under this constraint the conduction bands especially at the X k-point are modified compared with those of the strain-free Mg 2 Si (Figure 1). The changes in the conduction bands cause an increase of the indirect band gap up to 327 meV and a decrease of the direct band gap down to 1444 meV.

The electronic density of states (DOS) of both the energy-highest valence band (VB) and energy-lowest conduction band (CB) for unstrained and 2.246%-strained Mg 2 Si is presented in Figure 4(a). The comparison of the DOS between the strained and the unstrained compound shows that the strain does not noticeably affect the shape of the valence band near the Fermi level E f as the strained and unstrained Mg 2 Si VB curves are superimposed. They depart from each other only at -0.5 eV below E f . Overall the DOS of the valence band is slightly shrinked. From the figure 4(b) it can be seen that the convergence of the conduction bands results from both the upshift and downshift of the first lowest CBM (CBM1) and second lowest CBM (CBM2), respectively.

Seebeck Coefficient

This section presents the influence of the convergence of the conduction bands, caused by 2.246 % isotropic tensile strain, on the Seebeck coefficient of the Mg 2 Si compound. Since the relaxation time (τ) has been assumed to be constant, the Seebeck coefficient has been directly calculated from the density of states of the systems. This allows for a reliable comparison between the strained and strain-free Mg 2 Si bulk property. The Seebeck coefficient (S) is first given at predefined temperatures and various electron (n-type material) and hole (p-type material) doping levels. As a second step, the evolution of S with respect to temperature and at fixed doping levels (10 18 and 1.2×10 20 cm -3 ) is presented. These concentrations have been chosen from the fact that the Seebeck coefficient of strain-free Mg 2 Si exhibits its highest values at 300 K for 10 18 cm -3 (electrons and holes) and at 900 K for about 1.2×10 20 cm -3 (electrons and holes). [12] Fig. 3 Electronic band structure (left) and total density of states (right) of the strained bulk Mg 2 Si under 2.246% tensile strain.

The evolution of S at the three defined temperatures of 300, 600 and 900 K is presented in Figure 5, where the left and right panels refer to the hole and electron dopings, respectively. Irrespective of the type of doping the maximum of S decreases with increasing temperature and shifts toward high doping level. Indeed, at 300 K bulk Mg 2 Si shows a maximum S value of about 500 and -407 µV.K -1 at ∼10 18 cm -3 for hole and electron doping levels, respectively. As the temperature increases up to 900 K these maxima decrease down to about 224 and -225 µV.K -1 at ∼1.2 ×10 20 cm -3 for hole and electron doping levels, respectively.

Due to the marginal deformation on the valence bands, the Seebeck coefficient for hole doping (p) seems to be almost the same for the strained and unstrained Mg 2 Si, especially for p above 2×10 18 , 10 19 and 10 20 h.cm -3 for 300, 600 and 900 K, respectively (Figure 5 ac). We note that, at low temperature (300 K) and hole doping (p < 2 × 10 18 h.cm -3 ) the strained structure shows Seebeck values higher than the bulk Mg 2 Si, and amounts to 552 µV.K -1 . This increase could be due to the band gap enlarging. Since the convergence of the bands occurs in the conduction bands, a good enhancement of the Seebeck coefficient of the n-type material is expected. This is indeed what is observed in Figures 5 df ; Irrespective of the temperature and for n < 2 × 10 20 e.cm -3 the strained structure shows the highest values of S. In order to highlight further this feature, the improvement percentage of S max under isotropic tensile strain for nand p-type materials is plotted as a function of temperature in Figure 6. The largest improvement is in favor of the n-type material for which it reaches its maximum of 39% at 300 K.

As the temperature increases the improvement decreases down to below 5 % at high temperature. As can be seen in Figure 7, this feature is ascribed to the position of the chemical potential with respect to the bottom of the conduction band. At 300 K the chemical potential is 50 meV and 80 meV below the bottom of the CB for unstrained and strained Mg 2 Si, respectively, hence the improvement in S is large. At high temperature, the chemical potential for unstrained structure is inside the CB (40 meV above the bottom of the CB), which explains the low S value. Under isotropic constraint the gap increases, hence the chemical potential is shifted inside the band gap at 30 meV below the bottom of the CB. Therefore, S is improved, though less than at low temperature for which the distance to the bottom of the CB is larger (80 meV).

Figure 8 presents the evolution of the Seebeck coefficient as a function of temperature at constant hole and electron carrier concentrations of 10 18 cm -3 and 1.2 ×10 20 cm -3 , where the left and right panels refer to the hole and electron doping, respectively. For the strain-free structure, at low temperatures and low doping levels, the Seebeck coefficient reaches values of about 540 µV.K -1 for hole doping while for electron doping it reaches about -425 µV.K -1 (Figure 8 a,c). However, for temperatures above 400 K, due to the bipolar intrinsic regime, the Seebeck coefficient rapidly drops and levels off to a small negative value (-50 µV.K -1 ). For the heavy doping level, the maximum of S decreases compared with the low doping case and shifts toward the high temperature domain around 800 K (Figure 8 b, d).

Electrical conductivity

The evolution of the electrical conductivity with the doping level for the unstrained and 2.246%-strained structure is shown in Figure 9. Basically, there seems to be no change in the electrical conductivity whether Mg 2 Si undergoes constrains or not. In fact, a slight increase of σ of at most 17% (at 900 K and electron doping of 1.2x10 20 cm -3 ) is observed from the calculated data. As to the variation of σ with the doping level, we observe a constant increase, typical for semi-conductors and heavily doped semi-conductors.

Effective masses

The band effective mass (BEM), which is proportional to the inverse of the second derivative of the band energy (∝ 1/(dE 2 /d 2 k)), is known to play a crucial role in the effectiveness of materials for electronic transport properties. A large BEM tends to lower the carriers mobility, hence decreasing the electrical conductivity and the thermoelectric efficacy via the quality factor β . By constrast, large BEM contributes to large Seebeck coefficient. In Figure 10 is depicted conduction band structure around the X k-point and its successive derivatives for the unstrained and 2.246%-strained structure. It can be first noticed that, for unstrained Mg 2 Si both the CBM1 and CBM2 conduction bands have their minimum slightly shifted away from the X k-point by 0.01 Å -1 in the direction of Γ . Second, the mass of CBM1 is larger than that of CBM2 at X k-point (bottom panels in Figure 10). Upon strain, the CBM2 minimum remains slightly shifted from X whereas that of CBM1 is located at the X k-point.

The change in CBM1 and CBM2 mass induced by the strain is opposite for each band: the mass of CBM1 increases by about 65% whereas that of CBM2 decreases by about 20% (Figure 10). Overall, the DOS effective mass, given by m * = N contributes to the improved Seebeck coefficient when Mg 2 Si is subjected to the isotropic tensile strain (Figure 8 (c) and (d)).

Power factor

The power factor has been calculated from the results presented above for the electrical conductivity and the Seebeck coefficient. The improvement of the Seebeck coefficient is responsible for the increase in the power factor. This is particularly visible at low temperature (300 K) for electron doping for which, a 50% gain is observed (Figure 11(d)). As the temperature increases the gain becomes less (Figure 11(e)-(f)). In addition, under tensile strain the thermoelectric operating temperature of the material is downshifted. This type of behaviour has already been observed for Mg 2 Si under biaxial strains [23]. For hole doping, since the effective mass of the valence band is not modified under strain, the Seebeck coefficient remains unmodified, and so is the power factor, irrespective of the temperature (Figure 11(a)-(c)). Hence, we expect an improvement of the thermoelectric response of Mg 2 Si only in the case of electron doping.

Fig. 9 Electrical conductivity of the strain-free and 2.246%-strained Mg 2 Si as a function of hole and electron doping levels.

Figure of merit

To estimate the figure of merit zT , the knowledge of the thermal conductivity κ is indispensible. The thermal conductivity is contribute to by a lattice and an electronic parts. In semi-conductors it is generally accepted that the electronic contribution is minor compared to the lattice one, at most one-third of the total κ value at low and intermediate temperature (see e.g. Ref. [24]). The thermal conductivity has not been simulated in this work, hence we rely on experimental data already published (see for instance [24][25][26]), although these data have been determined for unstrained materials. Roughly, κ varies between 4 and 10 W m -1 K -1 [24] depending on the temperature and carrier concentrations. Another approximation has to be made in regards to the relaxation time τ of the electrons that intervenes into the calculated electrical conductivity. Recalling that the constant relaxation time approximation of the Boltzmann theory leaves τ undefined, the relaxation time can be determined by comparing the experimental electrical conductivity and the calculated one. According to Akasaka et al. [24] the electrical conductivity amounts to 1.05x10 5 S m -1 at 600 K for n=1.38x10 20 e.cm -3 . In the same conditions the calculated σ /τ value is 0.127x10 20 S m -1 s -1 for strain-free Mg 2 Si. Hence τ equals 8.265x10 -15 s. Assuming an electron-phonon T -dependence of σ and a n-dependence of τ the expression of τ is of the form [27]:

τ = τ 0 T -1 n -1/3 (2)
Using the above conditions for T and n (600 K and 1.38x10 20 e.cm -3 ) we obtain τ 0 = 2.563 10 -5 K cm s by inverting the previous equation. Because we do not know the experimental electrical conductivity for the tensilely strained structure, we cannot calculate τ 0 for strained Mg 2 Si. Hence, in the following we use the same relation to determine τ for the strain-free and strained conditions. It is consequently possible to determine the thermoelectric data σ and zT . Table 1 gathers all the data calculated from the above relations. As expected, the electrical conductivity increases with temperature whereas the Seebeck coefficient decreases. Overall, the figure of merit increases too, and is larger for strained than for unstrained Mg 2 Si at a given temperature. The improvement is particularly noticeable at 900 K where the increase is of about 25 %.

This behaviour is mainly caused by both the larger absolute value of S and the larger value of the electrical conductivity for the strained structure.

As compared with the experimental data [24] our calculated values for zT are smaller at 300 K and in good agreement at 600 K and 900 K (Table 1). The higher the temperature the smaller the difference. The deviation is probably caused by the various approximations used in our calculations.

Conclusion

The conduction bands overlap at the X k-point of Mg 2 Si has been achieved through a tensile isotropic strain of 2.246% applied to the structure. Based on DFT and Boltzmann transport theory formalisms a detailed analysis of the electronic properties of the strain-free and strained structure has been done. Under strain an enlargement of about 71 meV of the indirect band gap Γ -X has been observed. The curvature of the valence band is barely modified by the constraint whereas that of the conduction band is notably affected. Upon strain the effective mass of the electrons is increased by about 80%, which contributes to the improvement of the Seebeck coefficient as compared with the unstrained case. The enhancement of the Seebeck coefficient reaches up to 40 % at low temperature for the n-type material. The figure of merit zT has been calculated with the help of the experimental thermal conductivity and the theoretical electrical conductivity that has been calculated from the relaxation time estimated at various temperatures and electron concentrations. It has been observed that, at a given temperature, upon strain the electrical conductivity increases. Overall, the increase of both σ and S leads to a significant improvement of the figure of merit at all temperatures.

  work. The k-point selection is based on the Monkhorst-Pack scheme; the k-point mesh used to sample the Brillouin zone has been set to 20 × 20 × 20 for the cell parameter determination and 80 × 80 × 80 for the subsequent electron density calculations. The lattice constant of the bulk Mg 2 Si was optimized until the total energy has converged to at least 10 -9 Ry and the forces borne by the atoms have reached a value smaller than 10 -5

Fig. 1

 1 Fig. 1 Electronic band structure (left) and total density of states (right) of bulk Mg 2 Si.

Fig. 2

 2 Fig. 2 Calculated energies of the first and the second minimum of conduction bands at X-point as a function of the isotropic tensile strain.

Fig. 4

 4 Fig. 4 Density of states of the strain-free and 2.246%-strained Mg 2 Si.

Fig. 5

 5 Fig.5Seebeck coefficient of the strain-free and 2.246%-strained Mg 2 Si as a function of the electron and hole doping level at three defined temperatures of 300, 600 and 900 K.

Fig. 6

 6 Fig. 6 Improvement of the Seebeck coefficient with respect to temperature under isotropic tensile strain.

Fig. 7

 7 Fig. 7 Energy of the chemical potential in the band structure at various hole and electron concentrations and temperatures.

  b , where N v the pocket degeneracy and m * b the average BEM, is increased, which

Fig. 8

 8 Fig.8Seebeck coefficient of the strain-free and 2.246%-strained Mg 2 Si as a function of temperature for two electron and hole doping levels of 10 18 and 1.2×10 20 cm -3 .

Fig. 10

 10 Fig. 10 Conduction band structure (in eV, top panels), first derivative (in eV Å, middle panels) and second derivative (in eV Å 2 , bottom panels) for strain-free and 2.246%-strained Mg 2 Si. The two, lowest conduction bands (CBM1 and CBM2) are depicted at the vicinity of the X-point and the band structure is plotted along the Γ (left) and W (right) paths. The units along the abscissa are in Å -1 .

Fig. 11

 11 Fig.11Power factor (S 2 σ ) with respect to electron and hole doping levels for strain-free and 2.246%-strained Mg 2 Si at 300 K, 600 K and 900K.

Table 1

 1 Transport properties for strain-free and 2.246%-strained Mg 2 Si.

			unstrained			strained
	T (K)	300	600	900	300	600	900
	κ (W K -1 m -1 ) *	9.5	5.7	4	9.5	5.7	4
	n (cm -3 )	10 18	4 × 10 19	1.2 × 10 20	6 × 10 17	4 × 10 19	1.2 × 10 20
	τ (10 -14 s -1 ) * *	8.543	1.249	0.5773	10.128	1.249	0.5773
	σ (10 3 S m -1 ) * *	5.597	44.742	66.579	6.935	54.835	78.102
	S (µV K -1 ) * *	-388	-258	-225	-558	-299	-242
	zT * *	0.027	0.32	0.76	0.068	0.52	1.03
	zT *	0.05	0.25	0.7	-	-	-

* Experiment from

[24]

. * * Calculated, this work.

Ackowledgements

This work was granted access to the HPC resources of CINES (Centre Informatique National de l'Enseignement Supérieur, Montpellier, France) under Allocation C2015086881 made by GENCI (Grand Equipment National de Calcul Intensif). This work was also supported by the computing facilities of the CRCMM (Centre Régional de Compétences en Modélisation Moléculaire de Marseille) and by the mésocentre d'Aix-Marseille Université (Project 13b020). The Fondation EADS is acknowledged for financial support of H.B.'s PhD thesis.