Paul Eyméoud
email: paul.eymeoud@univ-amu.fr

Fabienne Ribeiro

Rémy Besson

Guy Tréglia

How to take into account local concentration in Ising-based Monte-Carlo: illustration with zirconium hydrides

Keywords:

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

Introduction

Many interatomic energetic models such as Embedded-Atom Method [START_REF] Daw | Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals[END_REF][START_REF] Daw | Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals[END_REF] (EAM), Cluster Expansion Model [START_REF] Sanchez | Generalized cluster description of multicomponent systems[END_REF] (CEM), or its simpler so-called Ising form [START_REF] Ising | Beitrag zur theorie des ferromagnetismus[END_REF], generally depend on the concentration of atomic species in the considered alloy or intermetallic. When that kind of concentrationdependent interatomic potential is used to implement a Monte-Carlo (MC) scheme, one can choose from two possible representations: (i) for all sites of the simulation box, affect the same concentration, computed by averaging compositions on the whole simulation box (the global concentration approach), (ii) for each atomic site of the simulation box, affect a specific local concentration, depending the atomic neighboring of this site (the local concentration approach).

When interatomic potentials are slowly varying with concentration, the global concentration approach (i), simpler, is generally sufficient to reproduce the correct physics for the considered system [START_REF] Mottet | Theoretical investigation of chemical and morphological ordering in pdccu1-c clusters[END_REF][START_REF] Lopes | Ordering and surface segregation in co1-cptc nanoparticles: A theoretical study from surface alloys to nanoalloys[END_REF][START_REF] Ludwig | Indium-gallium segregation in cuinxga1-xse2: An ab initio-based monte carlo study[END_REF][START_REF] Lavrentiev | Monte carlo study of thermodynamic properties and clustering in the bcc fe-cr system[END_REF]. However, in many situations, the concentration dependency of interatomic potential can be quite important. In that case, the local concentration (ii) is often required to correctly reproduce the physical behavior of the considered system. For instance, Shmakov et al. [START_REF] Shmakov | Decomposition kinetics in fe-cu dilute alloys. monte carlo simulation using concentration-dependent interactions[END_REF] have shown that the correct description of Fe-Cu alloys decomposition kinetics in CEM-based Monte-Carlo requires to take into account local concentration dependency of interactions. Similarly, Levesque et al. [START_REF] Levesque | Simple concentration-dependent pair interaction model for large-scale simulations of fe-cr alloys[END_REF] (resp. Stukowski, Caro et al. [START_REF] Stukowski | Efficient implementation of the concentration-dependent embedded atom method for molecular-dynamics and monte-carlo simulations[END_REF][START_REF] Caro | Classical many-body potential for concentrated alloys and the inversion of order in iron-chromium alloys[END_REF]) had to develop a local-concentration dependent pairwise interaction model (resp. EAM) implementing their Monte-Carlo approach, in order to correctly reproduce order vs segregation tendency and short-range order parameter trend with increasing at%Cr.

This question of taking or not local concentration dependency to reproduce the correct thermodynamics, is all the more important that in the case of Ising models, concen-tration dependency of pairwise interactions can sometimes be quite strong, with sign changes [START_REF] Ducastelle | Electronic structure and ordering[END_REF]. Furthermore, using clustering developments Gonis et al. [START_REF] Gonis | Use of the ising model in the study of substitutional alloys[END_REF] have demonstrated that concentration-independent interactions are rigorously valid only for strictly finite systems, such as a nanoparticle, whereas treatment of infinite systems, such as bulk alloys, requires an explicit concentration dependency of interactions. In the continuity of previously developed local concentration methodologies for interatomic potential MC studies of substitutional binary alloys [START_REF] Shmakov | Decomposition kinetics in fe-cu dilute alloys. monte carlo simulation using concentration-dependent interactions[END_REF][START_REF] Levesque | Simple concentration-dependent pair interaction model for large-scale simulations of fe-cr alloys[END_REF][START_REF] Pasichna | Alternative algorithms for simultaneous modeling of ordering and intermediate compound growth during reactive diffusion[END_REF], we propose here a detailed local concentration methodology for Isingbased MC studies of intermetallics. The procedure will be developed through the example of zirconium hydrides interstitial ordering, studied by a Tight-Binding Ising Model [START_REF] Tréglia | Segregation and ordering at surfaces of transition metal alloys: the tight-binding ising model[END_REF] (TBIM).

This article is organized as follows. First, we will briefly present the system chosen as a support for our study, namely zirconium hydrides, and the energetic model used for implementing our thermostatistical approach, TBIM. We will show the necessity, for this particular intermetallic system treated with this specific energetic model, of introducing the local concentration dependency of pairwise interactions in canonical Monte-Carlo Metropolis algorithm [START_REF] Metropolis | The monte carlo method[END_REF][START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], in order to correctly reproduce the ordering tendencies. This will lead us to develop a methodology for the consideration of local concentration in Monte-Carlo, and we will apply it to the construction of a phase diagram of hydrogen-vacancy ordering on interstitial tetrahedral sublattice of face-centered cubic Zr-H.

At that point, we can notice that detailed derivation of the Zr-H TBIM (recalled in section 2.2), and a preliminary Monte-Carlo study, have been already performed in Ref. [START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF]. However, the present investigation is clearly distinct from this previous work. In ref. [START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF], we failed at directly implementing the concentration-dependent TBIM in our MC with the global concentration approach, because of unphysical demixing phenomena. To overcome the problem, we implemented our global concentration approach MC with an alternative Ising model, fitted on DFT calculations by Connolly-

Type Stoichiometry

Bravais type for H occupancy of tetrahedral interstitial sites [START_REF] Suman | Hydrogen in zircaloy: Mechanism and its impact[END_REF] Zr cage-lattice [START_REF] Sidhu | Neutron and x-ray diffraction studies of nonstoichiometric metal hydrides[END_REF] FCC H atoms with no particular order [START_REF] Kolesnikov | Neutron scattering studies of ordered gamma-zrd[END_REF][START_REF] Bowman | Proton nmr line shapes in zrhx[END_REF][START_REF] Petrunin | Investigation of phase equilibria in zirconium deuterides[END_REF], potentially forming ordered patterns at low temperature [START_REF] Topchyan | Low-temperature phase transitions in the zirconiumhydrogen system[END_REF][START_REF] Bydlinskaya | New low-temperature phase transition in zirconium hydride[END_REF] ZrH 1.66-2.0 FCT H atoms nearly saturating tetrahedral sites [START_REF] Sidhu | Neutron and x-ray diffraction studies of nonstoichiometric metal hydrides[END_REF] Table 1 Structural characteristics of the three main zirconium hydrides.

X H 1 H 2 H 3 H 4 H 5 H 6 H 7 H 8 A 0.5 1 1 0 0 0 0 0 0 B 0.5 1 0 0 1 0 0 0 0 C 0.5 1 0 0 0 0 0 0 1 A 1.0 1 1 1 1 0 0 0 0 B 1.0 1 1 1 0 1 0 0 0 C 1.0 1 1 1 0 0 1 0 0 D 1.0 1 0 0 1 1 0 0 1 E 1.0 1 1 1 0 0 0 0 1 F 1.0 1 0 0 1 0 1 1 0 A 1.5 0 0 1 1 1 1 1 1 B 1.5 0 1 1 0 1 1 1 1 C 1.5 0 1 1 1 1 1 1 0 Table 2
Definition of a set of X ordered structures with respective ZrH □ 2-stoichiometries, using occupation factors of tetrahedral interstitial sites ( H equals 1 if site is occupied by an H atom, 0 otherwise). Elementary cell drawn with Ovito [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito-the open visualization tool[END_REF].

Williams inversion scheme [START_REF] Connolly | Density-functional theory applied to phase transformations in transition-metal alloys[END_REF]. In the present work, we manage to employ our strongly concentration-dependent TBIM in MC simulations, thanks to a local concentration methodology. This local concentration MC based on TBIM allows us to build a Zr-H schematic phase diagram (section 4), similar to the one obtained on Figure 9 of ref. [START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF] by global concentration MC based on DFT-fitted Ising model.

Position of the problem for Zr-H system

System under study

We will here focus on the case of Zr-H system, an intermetallic of particular interest in nuclear safety [START_REF] Suman | Hydrogen in zircaloy: Mechanism and its impact[END_REF][START_REF] Motta | Hydride formation in zirconium alloys[END_REF].

This binary system can form several kinds of sub-stoichiometric compounds, with a face-centered cubic (FCC) or tetragonalcentered cubic (FCT) zirconium lattice, and hydrogen atoms located on tetrahedral interstitial sites. Characteristics of these intermetallic hydrides are detailed in Table 1.

In our simulations, we will consider crystallographic variants with FCC Zr lattice at cell parameter fixed at 4.82Å (equilibrium parameter of fluorite-like ZrH 2 structure determined by first-principle calculation [START_REF] Dufresne | Tight-binding n-moment potential for zirconium hydride atomistic modeling[END_REF]), and various hydrogen concentrations and occupations. These variants of stoichiometry ZrH □ 2-, where represents the H/Zr ratio and □ an empty tetrahedral interstitial site, are described on Table 2. In our approach, this ZrH □ 2-system will be treated as a three-components alloys: zirconium Zr on FCC cage-lattice, hydrogen H on tetrahedral site, and vacancy □ on tetrahedral site.

Energetic model: Tight-Binding Ising Model (TBIM)

Here is a recap of procedure and results developed in Ref. [START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF][START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF], which will be necessary for what follows.

Starting point: the Tight-Binding (TB) Hamiltonian

The starting-point Hamiltonian model is the following Tight-Binding (TB) Hamiltonian, expressed on orthogonal atomic orbital basis {| , ⟩} ( , = , , and = Zr, H, □ respectively index site, atomic orbital and atomic type) :

= ∑ ∑ ∑ | , ⟩ ⟨ , | + ∑ , ≠ ∑ , ∑ , | , ⟩ ⟨ , | (1)
where denotes energy level of orbital of atomic type located on site , and hopping integrals from orbital of atomic type located on site to orbital of atomic type located on site (here, are assumed to be vanishing if and are not first-neighbors). The factor represents the atomic occupation : it equals 1 if site is occupied by atomic type , and zero otherwise.

This TB Hamiltonian relies on an spd parameterization splitting total energy into a sum of a band term (obtained by integration of density of states), and a repulsive Born-Mayer term [START_REF] Dufresne | Tight-binding n-moment potential for zirconium hydride atomistic modeling[END_REF][START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF] (assumed to be fixed at given hydrogen concentration and lattice geometry): tot = band + rep . Only the band term will be considered in what follows. Numerical values of energy levels and hopping integrals are given in Appendix A (see [START_REF] Dufresne | Tight-binding n-moment potential for zirconium hydride atomistic modeling[END_REF][START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF] for detailed fitting procedure).

Furthermore, computing partial electronic filling on , , orbitals, we find a good reproduction of charge transfer (see [START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF], Table 4.3), coherent with H -Zr electronegativity difference: partial filling at Fermi level near 1.2 (resp., lower than 4.0) for H (resp., for Zr).

Derivation procedure: the Generalized Perturbation Method

The energetic model we will use here to implement our Monte-Carlo approach is the Tight-Binding Ising Model (TBIM). Within this model, the band energy of an ordered variant X of Table 2 is expressed in the following way:

band (X ) = band (Ω ) + ord (X ) (2)
The first term of right member of equation ( 2), band (Ω ), represents the band energy of H/□ interstitial disordered state Ω . The second term, called ordering energy, is the little part of the energy characterizing the specific ordered state X , and is defined by:

ord (X ) = 1 2 ∑ , H -⟨ H ⟩ H -⟨ H ⟩ ( ) (3)
Within eq. ( 3), differences H -⟨ H ⟩ corresponds to the concentration fluctuation on site with respect to the disordered state, and ( ) denotes pairwise interaction between two hydrogen ℎ neighbors.

These effective pairwise interactions between hydrogen atoms ( ) appearing in equation ( 3) are derived from TB Hamiltonian (1) using the Generalized Perturbation Method [START_REF] Tréglia | Generalised perturbation theory in disordered transition metal alloys: application to the self-consistent calculation of ordering energies[END_REF][START_REF] Landesman | Electronic structure and pairwise interactions in substoichiometric transition metal carbides and nitrides[END_REF] (GPM) based on Coherent Potential Approximation [START_REF] Soven | Coherent-potential model of substitutional disordered alloys[END_REF][START_REF] Faulkner | Electronic states of substoichiometric compounds and application to palladium hydride[END_REF]

(CPA).

This technique falls into two parts. First, one needs to represent the H/□ interstitial disorder on the tetrahedral sublattice for a given composition , by resolving the following self-consistent equation (this is what we call "CPA"):

⋅ H -Σ( ) 1 -( H -Σ( )) ̄ ( ) -(2 -) ⋅ 1 ̄ ( ) = 0 (4)
where ̄ ( ) denotes the Green's function of the disordered effective medium, and Σ a complex potential representing the effective energy level of tetrahedral interstitial site. The second step (called "GPM") consists in taking the previously computed disordered state as a reference, and computing H-H pairwise interactions by perturbative calculation through the following analytic formula [START_REF] Landesman | Electronic structure and pairwise interactions in substoichiometric transition metal carbides and nitrides[END_REF] (for two H atoms located on sites 1 and 2 ):

( ) = - 2 Im ∫ F -∞ 1 1 -( H -Σ( ) ̄ ( )) ̄ ( ) 2 ⋅ 1 ⃗ ⃗ ′ ⋅ ∑ ⃗ , ⃗ ′ Tr HH 1 2 ( , ⃗ ) HH 2 1 ( , ⃗ ′ ) ⋅ ( ⃗ -⃗ ′ )⋅⃗ 1 2 ⋅ (5)
wherein ⃗ 1 2 denotes the distance vector between 1 and 2 sites, HH 1 2 ( , ⃗ ) the projection of Green's function ̄ ( ) on sites 1 , 2 at point ⃗ of first Brillouin zone, ⃗ the number of ⃗ points, and F the energy at Fermi level. Formula ( 5) can be generalized for triplets, quadruplets, etc, corresponding to higher order additional terms in perturbative development (3). In practice, we have proven [START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF][START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF] that for the Zr-H system, these multiplet interactions are negligible with respect to pairwise interactions.

Result: a concentration-dependent hydrogen-hydrogen pairwise interaction model

In Ref. [START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF][START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF], we have applied the methodology described in section 2.2.2 to compute the H-H pairwise interactions on a large range of concentration from = 0.0 to = 2.0 (see Table 4 in Appendix). The procedure led to three major conclusions: (i) H-H pairwise interactions are negligible beyond fourth neighbors, (ii) orders of magnitudes

(| 3 | >> | 1 |, | 2 |, | 4 |
) point out that ordering energies are mostly driven by third-neighbors H-H pairwise interactions, (iii) signs of pairwise interactions ( 1 > 0, 2 < 0, 3 < 0, 4 < 0) indicate that an hydrogen atom tends to be preferentially surrounded by vacancy at first neighborhood and hydrogen at second, third and fourth neighborhood. On top of that, we have noticed that pairwise interactions were presenting a substantial dependency in hydrogen concentration.

In order to explicit this concentration dependency of pairwise interactions, we have performed here a second-order polynomial fitting on data from Table 4, leading to the following analytic expressions in meV ( 3 corresponds to an average on the two kinds of third neighborhoods):

1 = 7.2 ⋅ -0.2 2 = 4.4 ⋅ 2 -15.4 ⋅ -6.1 3 = 4.4 ⋅ 2 + 5.7 ⋅ -70.4 (6)
4 = 3.9 ⋅ 2 + 1.5 ⋅ -16.8 Importance of correctly implementing this concentration dependency in calculations of ordering energy and canonical Monte-Carlo algorithms will be developed next.

Reproduction of ordering processes using

TBIM: order, disorder, segregation

Methodology: benchmark TBIM/DFT

In order to validate the ability of the TBIM constructed in section 2 to reproduce H/□ interstitial ordering processes on tetrahedral interstitial sub-lattice, we have performed a benchmark with Density Functional Theory [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistend equations including exchange and correlation effects[END_REF] (DFT) results. This ones have been performed at fixed geometry 2, and blue points to the disordered state Ω (resp. computed using CPA/SQS technique in TBIM/DFT approach). Red triangles identify the demixed state computed using eq. ( 7) from Zr and ZrH 2 energies, and red squares (resp. circles) correspond to the results obtained using global (resp. local) concentration formulation of ordering energy described in section 2.3.3.

(CFC Zr cell with lattice parameter = 4.82Å, fixed volume of cell, form of cell, and internal coordinates), using the Vienna Ab initio Simulation Package [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF], with 350 eV energy cut-off, ultrasoft pseudo-potentials [START_REF] Vanderbilt | Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[END_REF][START_REF] Kresse | Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements[END_REF], and 13 × 13 × 13 Γ-centered k-points grid. In order to represent the disordered state in DFT, we have used one 32-Zr atoms Special Quasirandom Structure [START_REF] Zunger | Special quasirandom structures[END_REF] (SQS) per concentration, built with the mcsqs modulus [START_REF] Van De Walle | Efficient stochastic generation of special quasirandom structures[END_REF] of ATAT code [START_REF] Van De Walle | The alloy theoretic automatic toolkit: A user guide[END_REF], taking into account clusters containing atomic pairs up to fourth neighboring length H-H pairs. Energy sequences of ordered variants of Table previously computed by both TBIM and DFT approach in ref. [START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF], are presented on Figure 1. On these graphics we have also added the novel results of disordered phase Ω (resp. computed using CPA/SQS technique in TBIM/DFT approach) and demixed state in order to evaluate the ordering tendency.

Ability of TBIM to reproduce ordering sequences, stability of order vs disorder

Considering Figure 1, we can draw two preliminary major conclusions. First, the procedure validates the ability of the TBIM to reproduce H/□ interstitial ordering processes, since it correctly reproduces the energy sequences of ordered variants obtained by DFT calculations (green points on graphics). Secondly, the ZrH □ 2-system tends to organize into interstitial ordered variants, rather than interstitial disorder, since for each concentration the energy of the disordered phase (blue cross on graphics) is higher than the energy of the most stable ordered variant: resp. C 0.5 , D 1.0 (isomorphic to the experimentally-observed phase [START_REF] Sidhu | Neutron and x-ray diffraction studies of nonstoichiometric metal hydrides[END_REF][START_REF] Kolesnikov | Neutron scattering studies of ordered gamma-zrd[END_REF]), C 1.5 for = 0.5, 1.0, 1.5.

For ZrH □ 2-, we are expecting three possible behav-iors for H/□ interstitial occupation of tetrahedral sub-lattice: (i) disorder (random occupation of H), (ii) order (stability of one or several ordered variants), (iii) phase separation (segregation of H on the one hand, and □ on the other hand).

We have shown the preferential stability of interstitial ordering with respect to interstitial disorder: it thus remains the question of the phase separation tendency.

Stability of order vs segregation: importance of local concentration

The energy of the demixed state at concentration can be obtained by computing the respective energy (Zr) and (ZrH 2 ) of two independent structures Zr and ZrH 2 , then applying the following formula:

( ) = 2 ⋅ (ZrH 2 ) + (1 - 2 ) ⋅ (Zr) (7)
In DFT, terms (Zr) and (ZrH 2 ) are computed on Zr and ZrH 2 cells. In TBIM, one applies formula [START_REF] Ludwig | Indium-gallium segregation in cuinxga1-xse2: An ab initio-based monte carlo study[END_REF] in cases = 0.0 and = 2.0, wherein the ordering energy (defined by eq. ( 3) is vanishing, to compute (Zr) and (ZrH 2 ). Both approaches (red triangles on graphics of Fig. 1) reveal that the energy of the demixed state is higher than the one of the most stable ordered variant (resp. C 0.5 , D 1.0 , C 1.5 at concentrations = 0.5, 1.0, 1.5). Thus, the ZrH □ 2-system tends to organize into interstitial ordered variants, rather than demix. This energetic comparison between segregation and ordering has been based on calculations on separated Zr and ZrH 2 structures. In practice, applying a TBIM-based Monte-Carlo on a supercell, demixing phenomena will occur in the following way: a part of interstitial sites of the simulation box saturated in H, and the other part unoccupied. To compute an average ordering energy on the simulation box, two formulations are possible.

The first one, called "global concentration approach", neglects the local concentration dependency of pairwise potential. In other words, in the hydrogen-rich region (ZrH 2like environment) of the demixed state , pairwise potential take the value ( ). With this formulation, in eq. ( 3) concentration corresponds to the ratio of number of H atoms in simulation box, on all interstitial sites of simulation box.

The second one, called "local concentration approach", takes into account local concentration dependency of pairwise potential. In other words, in the hydrogen-rich region (ZrH 2 -like environment) of the demixed state , pairwise potential take the value ( = 2.0) (corresponding to a local environment with local concentration = 2.0). With this formulation, in eq. ( 3) concentration is fixed to the value = 2.0 (local ZrH 2 environment).

Both approaches are numerically applied by computing ordering energy of demixed phase using eq. ( 3) (with or = 2.0 for global or local concentration approach), then subtracting the energy of the most stable ordered variant (resp. C 0.5 , D 1.0 , C 1.5 for = 0.5, 1.0, 1.5). One can see that global concentration formulation of ordering energy (red squares on Fig. 1, left) leads to a unphysical segregation tendency, whereas local concentration formulation of ordering energy (red circles on Fig. 1, left) is coherent with previous results (segregation disfavored with respect to ordering).

This observation enhances the importance of taking into account local concentration in TBIM-based studies of the ZrH □ 2-system, to reproduce the ordering tendencies, which correspond to the experimental stabilities of ordered phases from Table 1. We will see in the following that this last remark, due to the substantial dependency of pairwise potentials in concentration (noticed in section 2.2.3), will take a particular importance in our Monte-Carlo study.

Local concentration model for Monte-Carlo

Preliminaries

We will use here a canonical Monte-Carlo Metropolis algorithm (fixed number of atoms in the simulation box). Each simulation box contains 8000 interstitial tetrahedral sites, and takes into account periodic boundary conditions. The box geometry is fixed (rigid lattice), with lattice parameter of simple cubic (SC) interstitial lattice fixed to ∕2 with = 4.82Å. For each fixed temperature, we have performed 8 ⋅ 10 6 hydrogen-vacancy exchanges in simulation box.

Insufficiency of global concentration approach in Monte-Carlo, leading to phase separation at 0K

The global concentration in hydrogen on the simulation box, is defined as the ratio of the number of sites of the simulation box occupied by H atoms, on the total number of sites box = 8000 of the simulation box. It can be expressed using the following formula:

= 1 box ⋅ box ∑ =1 H ( 8
)
Before iterations: After 8 ⋅ 10 6 iterations: interstitial disorder H/□ demixtion wherein H denotes the occupation number of interstitial site (equals to 1 if site is occupied by an H atom, and 0 otherwise).

In our Monte-Carlo simulations, the global concentration approach consists in assigning the value ( ) (with defined by eq. ( 8)) to a pairwise interaction between two -th neighbors of the simulation box occupied by H atoms. With this formulation, the value of pairwise interaction does not depend on local hydrogen environment on the pair: for instance, a couple of two hydrogen atoms -th neighbors will have the same value of pairwise interaction, whether it is surrounded by hydrogen (local environment ZrH 2 ) or vacancies (local environment Zr).

Implementation of this global concentration Monte-Carlo approach at low temperature leads to a general tendency of the system to hydrogen segregation (see Figure 2 for illustration), whatever the imposed concentration, and even though the dimensions of simulation box are increased (double dimension in each direction , , ). This is in contradiction with experimental observations ( and stabilized phases of Table 1) and DFT results from section 2.3. Such a problem is caused by the necessity of correctly taking into account concentration dependency of pairwise potentials in order to avoid a general segregation tendency, which was previously evoked in section 2.3.3. For that reason, we will implement the local concentration dependency in Monte-Carlo, and detail the methodology.

Definition of the local concentration model

In section 2.2.3, we have fitted the analytic dependency of pairwise potentials in concentration (eq. ( 6)).

We will now define several particular sets of interstitial sites, represented on Figure 3.

The first ones are ( ): the site for = 0, and the sets of sites which are -th neighbors of on SC interstitial sublattice for > 0. Respectively, the cardinal of sets 0 ( ), 1 ( ), 2 ( ), 3 ( ), 4 ( ) equals to 1, 6, 12, 8, 6.

The second one is ̃ ( ), corresponding to the union of site and its sites first, second, third and fourth neighbors (Card( ̃ ( )) = 33):

̃ ( ) = 4 ⋃ =0 ( )
Considering the cluster ̃ ( ), one can define the local concentration as the ratio between the number of sites occupied by H atoms within this cluster, on the number of total sites of this cluster:

loc ( ) = ∑ ∈ ̃ ( ) H Card( ̃ ( )) = 4 ∑ =0 ∑ ∈ ( ) H 4 ∑ =0 Card( ( )) (9)
One can link local concentration (defined by eq. ( 9)) to global concentration (defined by eq. ( 8)) through the following formula, whose proof in given in Appendix C:

= 1 box ⋅ box ∑ =1 loc ( ) (10)
Within this formalism, one can define a local ordering energy for each site , depending on local concentration on this site:

loc ( ) = 1 2 ⋅ H ⋅ 4 ∑ =1 ∑ ∈ ( ) H ⋅ ( loc ( )) (11)
The global ordering energy of the simulation box can then be expressed as the arithmetical average on all local ordering energies of the simulation box:

ord = 1 box ⋅ box ∑ =1 loc ( ) (12)

Comments on the local concentration model

Equations ( 9), ( 10), ( 11) and ( 12), are defining our local concentration model. This last one is preserving process reversibility during an occupation exchange between two sites and of the simulation box with different initial atomic occupations ( H ≠ H ). In other words, if one performs the inverse exchange process (go back to initial occupation numbers), local concentrations take their values prior exchange: consequently, it will be the same for local ordering energies, and therefore for global ordering energy.

Moreover, in order to reduce CPU time, global ordering energy will be computed using formula [START_REF] Caro | Classical many-body potential for concentrated alloys and the inversion of order in iron-chromium alloys[END_REF] only once, before starting Metropolis algorithm. Then, at each proposition of exchange occupation between two sites and (with different initial occupations), we will add the energy variation Δ with respect to initial state. This quantity Δ will involve the sums of local ordering energies affected by the exchange, that is to say, the sum of ordering energies of sites located in neighboring environment of sites and :

Δ = ∑ ∈ ̃ ( )∪ ̃ ( ) loc ( ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ after exchange - ∑ ∈ ̃ ( )∪ ̃ ( ) loc ( ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ before exchange (13)

Reproduction of ordering tendencies at 0K using local concentration model

Contrary to the global concentration approach used in section 3.2, which has led to a general segregation tendency, our local concentration approach defined in section 3. Moreover, at the early stages of the iterations at 0K, some twinned domains can be observed: for instance, at = 1.0, germinating D 1.0 structures with different grain orientations. However, they reorganize in a single monovariant after a sufficient number of iterations (8 ⋅ 10 6 in our approach), leading to the quasi-perfect ordered structures presented on Fig. 4.

The stable variants D 1.0 , ZrH 2 , C 1.5 , C 0.5 can be respectively identified as phase (hydrogen occupation along (110) planes [START_REF] Sidhu | Neutron and x-ray diffraction studies of nonstoichiometric metal hydrides[END_REF][START_REF] Kolesnikov | Neutron scattering studies of ordered gamma-zrd[END_REF], see Table 1), phase (H saturating tetrahedral sites [START_REF] Sidhu | Neutron and x-ray diffraction studies of nonstoichiometric metal hydrides[END_REF], see Table 1), phase (low-temperature ordered variant with ZrH 1.5 stoichiometry [START_REF] Topchyan | Low-temperature phase transitions in the zirconiumhydrogen system[END_REF][START_REF] Bydlinskaya | New low-temperature phase transition in zirconium hydride[END_REF], see Table 1), phase (applying an HCP from FCC structural change on this phase with Zr 2 H stoichiometry [START_REF] Zhao | Identification d'une nouvelle phase d'hydrure de zirconium et modélisation à l'échelle mésoscopique de sa précipitation[END_REF]).

Application to the thermostatistical exploration of Zr-H phase diagram

We can use the local concentration model defined in section 3.3 in order to build a phase diagram of H/vacancy chemical order on interstitial tetrahedral sublattice, by canonical Monte-Carlo.

Selecting several concentrations between 0.0 and 2.0 (resp. = 0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75), we have performed a temperature loop for each concentration, starting from a randomly-distributed H/vacancy repartition at 0K. Formation of ordered phases are identified by following the evolution of both long-range order parameters, defined in Appendix D, and Warren-Cowley short-range order parameters, defined as the average number of first, second, third and fourth hydrogen neighbors of an hydrogen atom in the simulation box moy ( ∈ {1, 2, 3, 4}). For instance, one can observe on Figure 5 the order-disorder phase transition (of second-order type) for D 1.0 and C 1.5 ordered variants, confirming stability of phase (experimentally evidenced by [START_REF] Lanzani | Comments on the stability of zirconium hydride phases in zircaloy[END_REF][START_REF] Sidhu | Neutron and x-ray diffraction studies of nonstoichiometric metal hydrides[END_REF][START_REF] Kolesnikov | Neutron scattering studies of ordered gamma-zrd[END_REF], see Table 1) and low-temperature interstitial ordering of phase (experimentally evidenced by [START_REF] Topchyan | Low-temperature phase transitions in the zirconiumhydrogen system[END_REF][START_REF] Bydlinskaya | New low-temperature phase transition in zirconium hydride[END_REF], see Table 1).

The path lead us to identify several coexistence domains of ordered variants at low temperature: Zr+C 0. respectively around 900K, 950K, and 600K for , and phase.

At that time, writing [START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF] we were not able to use our TBIM to build the phase diagram by "classical" global concentration Monte-Carlo, since this approach was leading to a generalized hydrogen segregation tendency. To overcome this difficulty, we used an Ising model fitted on DFT calculations by Connolly-Williams inversion scheme [START_REF] Connolly | Density-functional theory applied to phase transformations in transition-metal alloys[END_REF].

With the new local concentration methodology presented in the present article, we have managed to correctly take into account the strong-concentration dependency of TBIM in Monte-Carlo, and stabilize, at low temperature, the same kinds of experimentally-detected ordered variants than the CEM-based global concentration Monte-Carlo of Ref. [START_REF] Eyméoud | Tight-binding modeling of interstitial ordering processes in metals: Application to zirconium hydrides[END_REF].

Conclusion and perspectives

We have first presented the system chosen as a support for our study, zirconium hydrides (section 2.1), and the energetic model chosen to implement the canonical Monte-Carlo, the Tight-Binding Ising Model (section 2.2). We have shown the importance, with those specific system and energetic model, to take into account the local concentration dependency for H-H pairwise interactions, in order to avoid a non-physical general segregation tendency (sections 2.3.3 and 3.2). We have then detailed the implementation of local concentration dependency in Monte-Carlo (sections 3.3 and 3.4). This local concentration Monte-Carlo has been then employed to build a schematic diagram of hydrogen-vacancy ordering on interstitial tetrahedral sublattice of Zr-H (section 4).

To sum up, the main interest of this work is to give a ready-to-use Monte-Carlo methodology for treating local concentration, which can be applied for any alloy or intermetallic whose ordering energy is described by an Ising model with strongly concentration-dependent pairwise interactions (which is often the case for substitutional binary alloys [START_REF] Ducastelle | Electronic structure and ordering[END_REF], for instance). The path can be also enlarged for more complex interatomic potentials. Particularly, several CEM [START_REF] Holliger | Hexagonal-based ordered phases in h-zr[END_REF][START_REF] Holliger | Reciprocal-space cluster expansions for complex alloys with long-range interactions[END_REF][START_REF] Besson | Ab initio thermodynamics of fcc h-zr and formation of hydrides[END_REF], derived by Connolly-Williams inversion scheme by fits on DFT results, were developed to study interstitial ordering of zirconium hydrides by global concentration MC. It would be interesting to introduce the local concentration in these MC investigations, and evaluate its impact on phases stability.

Beyond the case of MC bulk studies of binary materials with strong concentration dependency of interatomic potentials (aforesaid intermetallic hydrides [START_REF] Holliger | Hexagonal-based ordered phases in h-zr[END_REF][START_REF] Holliger | Reciprocal-space cluster expansions for complex alloys with long-range interactions[END_REF][START_REF] Besson | Ab initio thermodynamics of fcc h-zr and formation of hydrides[END_REF] and substitutional alloys [START_REF] Shmakov | Decomposition kinetics in fe-cu dilute alloys. monte carlo simulation using concentration-dependent interactions[END_REF][START_REF] Levesque | Simple concentration-dependent pair interaction model for large-scale simulations of fe-cr alloys[END_REF]), the methodology developed here can be broadened to more complex MC investigations. For instance, it can be employed to study the segregation phenomena around dislocations, implying important variations of local concentration: C around steels dislocations [START_REF] Lüthi | Ab initio thermodynamics of carbon segregation on dislocation cores in bcc iron[END_REF], P and Cr in -iron dislocation core [START_REF] Medouni | Role of dislocation elastic field on impurity segre-gation in fe-based alloys[END_REF], etc. It can also be useful for investigations on high-entropy alloys, wherein local concentration fluctuations could be quite significant [START_REF] Rao | Atomistic simulations of dislocations in a model bcc multicomponent concentrated solid solution alloy[END_REF]. Pairwise interactions values between first, second, third and fourth order neighboring hydrogen atoms, with hydrogen concentration , extracted from [START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF].

A. Numerical implementation of TB Hamiltonian

Table 3 gives the numerical parameters for TB Hamiltonian defined by eq. [START_REF] Daw | Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals[END_REF].

Detailed fitting procedure is given in Ref. [START_REF] Dufresne | Tight-binding n-moment potential for zirconium hydride atomistic modeling[END_REF][START_REF] Eyméoud | Modélisation atomistique de la fragilisation des gainages combustibles nucléaires par les hydrures : caractérisation de l'ordre chimique interstitiel des hydrures de zirconium à l'aide d'un modèle d'Ising effectif dérivé des liaisons fortes[END_REF].

B. Data for TBIM fitting

The data employed for the fitting of formula ( 6) is given in Table 4.

C. Relation between global and local concentration: mathematical proof

Here, for all subset of sites of the simulation box, the notation ℭ box ( ) designs the complement of in the simulation box.

Using definition [START_REF] Shmakov | Decomposition kinetics in fe-cu dilute alloys. monte carlo simulation using concentration-dependent interactions[END_REF] of local concentration, one can write: Let a site of the simulation box, this site has 32 neighbors, which are elements of the set ̃ ( ) ∩ ℭ box ({ }).

1 box ⋅ box ∑ =1 loc ( ) = 1 box ⋅ 1 33 ⋅ box ∑ =1 ∑ ∈ ̃ ( ) H ⏟⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞ ⏟ = ( 14
)
Reciprocally, there are 32 sites , different from , which contain in their set ̃ ( ).

Therefore, among the 33 ⋅ box terms of double sum , there are 33 terms taking the value H .

Since this argument is valid for each sites of the simulation box ( ∈ {1, ⋯ , box }), one can rearrange these terms to write:

= box ∑ =1 33 ⋅ H (15)
Injecting [START_REF] Pasichna | Alternative algorithms for simultaneous modeling of ordering and intermediate compound growth during reactive diffusion[END_REF] in equality ( 14) one gets:

1 box ⋅ box ∑ =1 loc ( ) = 1 box ⋅ box ∑ =1 H ( 16
)
By identification with definition (8) of global concentration , one then gets the following equation:

1 box ⋅ box ∑ =1 loc ( ) = (17)
corresponding to formula [START_REF] Levesque | Simple concentration-dependent pair interaction model for large-scale simulations of fe-cr alloys[END_REF].

D. Definition of long-range order parameters

Order parameter for D 1.0 (resp. C 1.5 ) variant can be defined by dividing interstitial tetrahedral sublattice within two (resp. four) sublattices, as represented on Figure 7, then defining linear combinations of concentrations on these sublattices:

(D 1.0 ) = - (C 1.5 ) = 1 3 ( 2 + 3 + 4 ) -1
As concerns C 0.5 variant, one uses the definition of (C 1.5 ) with permutations H ⟷ □.

Figure 1 :

 1 Figure 1: Sequences of formation energies with respect to the most stable ordered state (resp. Zr, C 0.5 , D 1.0 , C 1.5 , ZrH 2 at concentrations = 0.0, 0.5, 1.0, 1.5, 2.0): comparison between TBIM and DFT results. Green points correspond to ordered variants of Table2, and blue points to the disordered state Ω (resp. computed using CPA/SQS technique in TBIM/DFT approach). Red triangles identify the demixed state computed using eq. (7) from Zr and ZrH 2 energies, and red squares (resp. circles) correspond to the results obtained using global (resp. local) concentration formulation of ordering energy described in section 2.3.3.

Figure 2 :

 2 Figure 2: Global concentration approach at = 1.0 and 0K: atomistic view of the simulation box (interstitial sublattice only: green for H atom, red for vacancy) before and after convergence.

Figure 3 :

 3 Figure 3: Schematic representation of ensemblist description used on the SC tetrahedral interstitial sublattice.

= 0. 25 : 2 Figure 4 :

 2524 Figure 4: Simulation results at 0K with the local concentration approach: atomistic view of the simulation box (interstitial sublattice only: green for H atom, red for vacancy).

1 Figure 5 :C 1 . 5 moy

 1515 Figure 5: Order-disorder phase transition for D 1.0 and C 1.5 variants, identified by plotting long-range and short-range order parameters with respect to temperature. In short-range approach, one uses the average number of -th H neighbors of an H atom in the simulation box, in comparison with the known values for ordered variants D 1.0 moy , C 1.5 moy (continuous lines) and disordered phase Ω moy = ∕2 ⋅ ZrH 2moy (dotted lines). In long-range approach, one uses parameters defined in Appendix D (only the heating part of temperature loop is presented in this case).

  5 , C 0.5 +D 1.0 , D 1.0 + C 1.5 , C 1.5 + ZrH 2 for respective concentration ranges ∈ [0.00; 0.50], [0.50; 1.00], [1.00; 1.50], [1.50; 2.00]. With the correspondences established in section 3.5, these domains respectively correspond to phase coexistence Zr + , + , + , + . Such considerations are summarized on the schematic phase diagram of Figure 6. This last one is coherent with the diagram presented on Figure 9 of Ref. [19]: stabilization of + , + , + domains, and order-disorder temperature

Figure 6 :

 6 Figure 6: Phase diagram of H-vacancy interstitial ordering on the tetrahedral interstitial sublattice, resulting from our local concentration Monte-Carlo. Experimental data from [47].

3

 3 Energy levels (in eV) and Slater-Koster parameters ( ) = 0 ij exp[-( -)∕ )] (in eV) for hopping integrals. is the interatomic distance between first neighbors ( ZrZr = 3.20Å, and ZrH = 2.09Å), ZrZr = 2.40, and ZrH = 1.8.

Figure 7 :

 7 Figure 7: Representation of specific sublattices used to define long-range order parameters.

Table 4

 4

	4 (meV)

1 (meV) 2 (meV) 3 (meV)

Acknowledgments

The authors thank ANR SIZHYP for financial support. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.

Data availability statement

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

CRediT authorship contribution statement

Paul Eyméoud: Conceptualization, Methodology, Software, Writing -Original Draft, Writing -Review and Editing. Fabienne Ribeiro: Supervision. Rémy Besson: Software. Guy Tréglia: Supervision.