Modeling of Phase Equilibria in Ni-H: Bridging the Atomistic with the Continuum Scale

: In this paper, we present a model which allows bridging the atomistic description of two-phase systems to the continuum level, using Ni-H as a model system. Considering conﬁgurational entropy, an attractive hydrogen–hydrogen interaction, mechanical deformations and interfacial effects, we obtained a fully quantitative agreement in the chemical potential, without the need for any additional adjustable parameter. We ﬁnd that nonlinear elastic effects are crucial for a complete understanding of constant volume phase coexistence, and predict the phase diagram with and without elastic effects.


Introduction
The understanding of phase equilibria and transitions is essential for the comprehension of many processes in nature as well for a theory guided search for novel materials with superior properties. The modeling of two-phase coexistence in binary alloys is based on the common tangent construction: In equilibrium, the chemical potentials of each species in the two phases of interest have to be equal. As shown in pioneering work by several authors [1][2][3], elastic effects can significantly modify the phase coexistence behavior, especially in solid phases, e.g., due to density differences.
For a true multi-scale modeling of phase equilibria and transitions in complex materials, ranging from macroscopic dimensions down to the nanoscale, an efficient and accurate matching between the atomistic simulations and formal thermodynamic and continuum concepts is critical. However, despite significant progress, there is still a substantial gap between these two levels. The purpose of the present article is therefore to seamlessly connect the atomistic and continuum scale, illustrated for the Ni-H system. This system is used as a prototype in atomistic simulations for the understanding of hydrogen embrittlement phenomena [4] and rechargeable batteries [5][6][7]. Mechanical deformations play a crucial role, and thus a careful thermodynamic inspection is essential for a thorough understanding. In previous work [8], we connected ab initio modeling of coherent phase equilibria in the presence of lattice strains and interfacial proximity effects to continuum descriptions. Here, in contrast, the focus is on a matching between Monte Carlo (MC) simulations using empirical potentials and bulk effects on the continuum level. Additionally, we consider interfacial energy contributions, which were neglected in [8].
Experimentally, it is known that metallic nickel absorbs very limited amounts of hydrogen under ordinary conditions, but is known to form a nearly stoichiometric hydride NiH at high hydrogen fugacities. Fukai et al. performed in situ measurements of the lattice parameter of Ni-H alloys over a wide range of hydrogen partial pressures. This way they constructed a pressure-composition-temperature diagram, encompassing the region of two-phase coexistence [9,10].
One of the central outcomes of the present study is the importance on nonlinear elastic effects in constant volume simulations, as large stresses can arise [11]. Moreover, our approach also provides important fundamental insights into the theory of phase equilibria in coherent solid-state systems, as it elucidates quantitatively the role of different energetic and entropic contributions.
The present work is complementary to atomistic modeling [12] and cluster expansion approaches [13,14], which aim at a description on the atomic level. Cluster expansion approaches are based on an energy description of an Ising-like model of alloys and are powerful methods for predicting phase diagrams and ordering phenomena; long-ranged elastic effects, however, are difficult to capture in such approaches. Therefore, in the present work, we focused on the bridging between a hybrid molecular static and Monte Carlo approach for finding the equilibrium configurations on the one hand and long-wavelength continuum theories on the other hand. This link across the scales is needed to capture in a closed description both microscopic effects related to H-H interactions as well as long-ranged elastic effects due to phase separation. The derived scale bridging energy functional can then serve as a direct input e.g., to phase field simulations. The resulting quantitative scale-bridging is the major outcome of this paper.
Step by step, we introduce the various contributions on the continuum scale to match the atomistically determined chemical potential of hydrogen, which is a key quantity in the description of phase coexistence.

Methods
The starting point of our study is the determination of the equilibrium spatial distribution of the interstitial H atoms in the metallic matrix, employing grand canonical Monte-Carlo simulations and molecular statics. The simulations are performed within a periodic 10 × 10 × 10 fcc supercell containing 4000 Ni atoms at constant volume (corresponding to the equilibrium volume of pure Ni) and temperature. H atoms are introduced on the octahedral interstitial sites according to the Metropolis algorithm [15]. The energies of the trial steps are determined by relaxing the ionic degrees of freedom-thus fully capturing elastic interactions-using the LAMMPS simulation package [16,17] in conjunction with a modified version [4] of the Ni-H EAM potential by Angelo et al. [18,19]. The effects of temperature are thus only considered to account for the configurational entropy contributions to the free energy, while vibrational contributions are disregarded. This description has previously been used to model hydrogen mediated dislocation-dislocation interaction [4].
The continuum modeling involves the construction of a free energy functional, which is derived step by step in the following section. This partly requires solving continuum mechanics problems, for which we used in particular finite element modeling using ABAQUS.

Monte Carlo Modeling
In the grand canonical simulations, the chemical potential is prescribed, and the number of hydrogen atoms can vary. Depending on the H chemical potential dilute or condensed H distributions are obtained. The condensed hydride precipitates remain coherent and adopt characteristic shapes depending on the bulk H concentration, as shown in Figure 1.
A peculiarity, which motives the present scale bridging investigations, is the functional form of the hydrogen chemical potential µ H as function of the average H concentration, see dots in Figure 2, which will be discussed in detail below. Classically, we would expect (in a strain free situation) the chemical potential to be constant in the entire two-phase region as a consequence of the common tangent construction. Obviously, this is not the case here, and it is one of the goals of this paper to shed light on this effect.

Free Energy Formulation
Our starting point for transferring the atomic scale behavior to the continuum level is the free energy for the single phase material, with the number of hydrogen atoms N H and the elastic free energy F el , the configurational free energy F c and the H-H interaction F H-H . µ 0 is the solvation energy needed to insert an isolated hydrogen atom into the (empty) matrix, in contrast to the aforementioned chemical potential µ H , which also includes mutual interactions, elastic effects, etc.

Solvation Energy and H-H Interaction
For the formation of a hydride phase, it is essential to include a description of the lattice mediated interaction between the hydrogen atoms [4]. We use the Margules model [5], F H-H /N Ni = −αc 2 /2 + βc 3 /3, where N Ni is the number of nickel atoms and c = N H /N Ni the (homogeneous) hydrogen concentration; it is normalized to 1 for a crystal where all interstitial octahedral sites are occupied. The parameters are determined from the Monte Carlo calculations by averaging over several random mixtures of fully relaxed homogeneous solutions. We point out that, for this matching, we only consider homogeneous "lattice gas" configurations and not phase separated states, which involve additionally (coherency) elastic and interfacial effects, and which will be parametrized below. These reference calculations are performed under constant pressure conditions, P = 0, to suppress (external) stress effects. The resulting curve F(c) is shown in Figure 3. It is fitted by a polynomial F f it /N Ni = a 0 + a 1 c + a 2 c 2 + a 3 c 3 . Identification with the above parameters therefore uniquely gives µ 0 = a 1 , α = −2a 2 and β = 3a 3 (a 0 is an irrelevant zero-point energy). Hence, we obtain µ 0 = −2.148 eV relative to monatomic hydrogen in vacuum, α = 0.751 eV and β = 0.355 eV. These numbers differ slightly from those given in Ref. [5] due to a different pair potential cutoff to stabilize the hydride, as this leads to a positive value of the elastic constant C 44 , see Ref. [4] for details. The positive sign of α favors an increase of the local hydrogen concentration, thus leading to the hydride formation.  The chemical potential is defined in the usual way, and is split into additive contributions analogous to the free energy in Equation (1). We therefore get for the chemical potential contribution by the H-H interaction µ H-H = −αc + βc 2 . It is shown in Figure 4 for both the present fitting parameters and the ones given in Ref. [5], exhibiting only small differences between the predictions.

Configurational Entropy
The configurational free energy F c stems from the different possibilities to occupy the interstitial sites with hydrogen (or the hydride with vacancies), and is therefore given by under the assumption of equal occupation probability for all octahedral sites [20]. In the low concentration regime, it gives the dominant nontrivial contribution, as elastic and interaction effects are negligible there, see Figure 5. We note that the configurational term is the only temperature dependent one in the present description. To verify this dependence we also computed in the Monte Carlo simulations cases with different temperatures, and the results are shown in Figure 6, showing very good agreement.

Maxwell Construction
Thus far, the system is described for a single phase state only, and elastic effects are suppressed by a free volume expansion. Due to homogeneity, elastic stresses as well as interfacial contributions do not appear, and the description is therefore complete on this level.
According to the usual thermodynamic picture, a single phase state is unstable in regions where the slope of the chemical potential is negative, µ H (c) < 0, and phase separation should occur there. As can be seen in Figure 6, this is the case for a wide concentration regime for low temperatures, whereas for high temperatures the hydrogen solubility limit is significantly larger.    Phase coexistence on this level is described by Maxwell's equal area rule, which states that phase separation sets in at the intersection points of a horizontal line with the S-shaped van der Waals loop, cutting it into two equal areas above and below this Maxwell line. This is shown in Figure 2 for T = 300 K. From this, we get the constant chemical potential in the two-phase region, µ M = −2.405 eV for T = 300 K. In equilibrium, the system therefore enters the two-phase region at the first intersection of the horizontal Maxwell line with the S-shaped single phase curve. Apparently, this happens already for very low concentrations in the ppm range. Nevertheless, on the following ascending part of the curve µ H (c), the system is still metastable. We point out that, according to the principles of the two-phase construction, the dual phase hydrogen chemical potential should be constant and equal to µ M until the phase transformation is completed. Obviously, this is not the case for the atomistic data with the fixed volume constraint.

Elastic Effects
Before entering into a detailed analysis of elastic effects, let us briefly discuss them in relation to the H-H interaction term. The latter is sometimes also introduced as long-range elastic interaction, and we want to stress that we do not double count effects here. The H-H interaction term expresses that an H atom locally deforms the lattice, and therefore makes the placement of another H atom in the vicinity energetically favorable. In a mean field context, this is expressed by the term ∼−c 2 in F H-H . If the site occupancy gets higher and the octahedral sites more and more filled, the energy for placing more H atoms into the lattice increases again, as expressed by the +c 3 term. Although this consideration makes reference to elastic deformations, this is distinct from the elastic effects considered in the following. To make this point more explicit, we consider the following two situations, for simplicity both with homogeneous hydrogen distributions, i.e., spatially constant concentration c: First, for fixed pressure P = 0, we increase the hydrogen concentration. This leads to a widening of the lattice, as discussed in the following and shown in Figure 7. However, since the system remains stress free on a mesoscopic or macroscopic level, σ ij = 0, the elastic energy (see Equation (5)) in the continuum mechanics sense remains zero. Nevertheless, the total energy of the system is changed due to the H-H interaction, and this is expressed through F H-H , which is independent of the elastic stress state.
Second, for fixed concentration (and total number of hydrogen atoms), we change the external pressure of the system. Then, the H-H interaction term does not change, but the mesoscopic elastic energy does.
The difference between the H-H interaction term and the explicit elastic term becomes most prominent in two-phase situations. Then, the hydride has a larger equilibrium lattice constant and therefore distorts the surrounding, coherently connected matrix phase. This leads to long-ranged elastic deformations, and the range of these interactions of the order of the precipitate size, which is captured by the mesoscopic elastic energy term introduced below, but not by the H-H interaction term. For further discussion of this separation of microscopic and mesoscopic elastic contributions, we refer to [8]. There, it was shown that the elastic energy due to microscopic deformations around individual H atoms and long ranged strains resulting from the formation of mesoscopic precipitates decompose additively into microscopic and mesoscopic contributions without the appearance of a cross term.
The large deviations between the atomistic data and the continuum model clearly indicate that elastic and interfacial effects are critical and cannot be neglected. We remind that the basis for the Maxwell construction-or here equivalently the common tangent construction-is based on the assumption that the two-phase energy in a macroscopic system decomposes additively into the contributions of the two phases, which do not influence each other. This means, that they are described by free energy densities, which depend only on the local concentration. This condition, however, is violated in the presence of elastic effects. Here, e.g., a density variation in one phase also influences the other, as the global pressure in the system changes.
In the following, we discuss the modeling mainly from the continuum perspective. The technical aspects of the parameter matching between the scales are given in the appendices.
We focus first on the low-temperature regime, where the hydrogen solubility limit is low (see Figure 13), and higher temperatures will be discussed later. In the present regime, the phases are almost stoichiometric, and the chemical potential of the hydrogen in the two-phase region is modified. From a minimization of the total free energy in the two-phase region, we obtain the modified chemical potential (see Appendix A) In contrast to the Maxwell term µ M , the elastic and interfacial contributions, µ el and µ s , respectively, do depend on the concentration, and therefore the chemical potential is no longer constant in the two-phase region, as observed in the Monte Carlo simulations (see Figure 2).

Linear Elastic Effects
The pure nickel and the hydride exhibit a substantial lattice mismatch, and therefore, in the two-phase region, elastic stresses arise. The linear elastic free energy density per unit volume is given by where ij = (∂ i u j + ∂ j u i )/2 is the strain derived from the displacements u i , and the eigenstrain 0 ij = ε 0 δ ij represents the isotropic volume expansion of the hydride [21]. It is concentration dependent, and we extracted it again from T = 0 K simulations of homogeneous states with P = 0; it is well 0189, going beyond a linear dependence in the spirit of Vegard's law. The eigenstrain is directly related to the concentration dependent lattice constants with a Ni = 3.520 and a Ni−H = 3.738 via ε 0 (c) = [a(c) − a Ni ]/a Ni , since we use here the relaxed, hydrogen-free nickel as reference configuration (see Figure 7).  The elastic constants for the pure phases are listed in Table 1. We note that all values are given with respect to Ni as reference configuration; see Appendix B for a discussion of this issue. The determination of the elastic constants from the Monte Carlo simulations is described in detail in Appendix C. Table 1. Elastic constants of pure nickel and the fully saturated hydride, as obtained from the EAM potential (see Appendix C for details). We note that the relaxed configuration of Ni is used as reference state. To better understand the role of the elastic effects, it is instructive to inspect an analytical isotropic linear elastic model, where we assume for simplicity the elastic constants to be equal in both phases. We use a spherical sample of radius R consisting of nickel, with a hydride inclusion of radius R H in its center. The radial displacement component is continuous at the coherent interface and vanishes at the outer boundary due to the volume constraint. At low temperatures, the concentration is in good approximation c = 0 in Ni and c = 1 in the hydride. This gives (see Appendix D) where N 0 = 4 is the number of octahedral sites per unit cell, ∆ε 0 = ε 0 (1) − ε 0 (0), and λ and G are Lamé coefficient and shear modulus respectively. In this approximation, the elastic contribution in Equation (6) is a linear function of the concentration, and it depends quadratically on the eigenstrain difference ∆ε 0 . Since the material is relaxed for c = 0, the elastic effects therefore give only an additive constant for c = 0, which is related to the elastic hysteresis due to the internal stresses that arise upon precipitate formation. For the isotropic elastic constants we use the values for nickel and λ = C 12 and G = (C 11 − C 12 )/2 for the solid curve µ M + µ lin.el in Figure 8.
Isotropic, Ni constants Isotropic Cubic Nonlinear elastic MC Figure 8. Elastic contribution to the chemical potential in various approximations: The solid curve expresses the isotropic prediction by Equation (6), using equal elastic constants for both phases (we use the values of pure Ni). The long-dashed curve is the same isotropic model, but this time using the different elastic constants of both phases. This prediction is almost identical to a numerically calculated expression based on cubic elasticity, which assumes a spherical inclusion in a cubic box with periodic boundary conditions (short-dashed curve). In contrast to all these linear elastic approximations, the dotted line also considers elastic nonlinearities. Its prediction is very close to the Monte Carlo data.
The calculation can also be done using the elastic constants of the individual phases, leading to the long-dashed curve shown in the same graph (see Appendix D). A numerical finite element solution with ABAQUS for spherical precipitates in a cubic box with fixed volume, taking into account the full cubic elasticity and periodic boundary conditions, leads to the short-dashed curve shown in Figure 8. It is very close to the analytical results. Since we still observe a significant discrepancy to the Monte Carlo data, we conclude that the consideration of linear elasticity is not sufficient to explain the slope of the chemical potential in the two-phase region.

Nonlinear Elastic Effects
The reason for the discrepancy between the atomistic data and the continuum modeling involving linear elasticity is the appearance of large compressive stresses for higher hydrogen concentrations (at fixed volume), and therefore the linear elastic approximation breaks down. Instead, we have to take into account nonlinear elastic effects; they include nonlinearities in the stress-strain relationship, leading effectively to expressions C ijkl ({ mn }) in the elastic free energy density. Additionally, geometrical nonlinearities appear due to the volume deformation and the nonlinear strain tensor [11], We have estimated that these effects are small in comparison to the change in the constitutive law, and therefore effectively captured them in the modified elastic constants. However, it turns out to be important to carefully use consistently the same reference configuration for the Lagrangian formulation of elasticity for both phases. We note that the chemical potential contains now also a contribution due to the strain dependence of the elastic constants. We mention in passing that we do not observe plastic deformations in the atomistic simulations. From a careful analysis of the Monte Carlo simulations we find that the nonlinear corrections depend only on the volume change and not on shear effects, thus the elastic constants depend only on the trace of the (local) strain tensor, which simplifies the expressions. We therefore write them as series expansion for all elastic constants or combinations of them. From the atomistic simulations we find the coefficients of the nonlinearities for C 11 − C 12 , C 44 and the bulk modulus K, see Table 2 and Appendix C. From them, all elastic constants can be expressed using the relation 3K = C 11 + 2C 12 . We note that for the present geometry bulk compression is the dominant effect, and the material becomes stiffer in this regime of negative strain (compression). Including the nonlinear contributions the continuum chemical potential shows now a very satisfactory agreement with the Monte Carlo data, see the dotted line in Figure 8.

Interfacial Effects
In a final step, we take into account interfacial effects. As already visible from the agreement of the Monte Carlo data and the nonlinear elastic curve, their contribution is obviously small. In isotropic approximation, the surface energy for a spherical inclusion is given by F s = 4πγR 2 H with the surface energy γ. From planar interface calculations, we obtain γ = 0.11 Jm −2 for (100) interfaces and γ = 0.17 Jm −2 for (111) interfaces. From the fitting of two-phase data, we independently get an excellent agreement using γ = 0.17 Jm −2 for the spherical inclusion.
The magnitude of interfacial terms should be compared to bulk contributions to see their influence and to get an impression on the dimension of critical nuclei. For an order of magnitude estimate, we compare elastic bulk contributions in the approximation in Equation (6) to interfacial contributions for spherical precipitates for small concentrations in Equation (A42), where the influence of interfacial terms is largest. For c 1 the elastic and interfacial contributions are comparable for which is the case for about 23 H atoms. This indicates that interfacial contributions are indeed small for the present scenario in comparison to the elastic bulk terms.
For a linear isotropic material with equal elastic constants in both phases and a purely dilatational eigenstrain, the elastic energy does not depend on the shape of the inclusion but only on its volume fraction [22]. Although the conditions for this rigorous statement are not exactly fulfilled here, we find only a small dependence of the elastic energy on the precipitate shape in the Monte Carlo data, as exemplarily shown for spherical hydride inclusions in a Nickel matrix and vice versa in Figure 9. Consequently, the shape of the precipitate is determined by the interfacial energy terms alone, see Appendix E for explicit calculations. This means, that for each given average concentration c, which is related to a precipitate size, the equilibrium shape with the minimum interfacial energy appears. For isotropic surface energy we expect spherical inclusions for c < 4π/81 ≈ 0.16, followed by tubes up to c = 1/π ≈ 0.32, and slabs thereafter. The situation is symmetric for c > 0.5 with then the hydride being the matrix phase. Thus we expect nickel tubes in the range 1 − 1/π < c < 1 − 4π/81 and spheres beyond this concentration. In our Monte Carlo calculations, we indeed find all these structures in the correct ordering (see Figures 1 and 10), but the transition points differ slightly due to anisotropy and nonlinear elastic effects. Based on the observed structures and the calculated interfacial energy we added this contribution to the nonlinear elastic chemical potential in Figure 10. The transitions between different precipitate shapes explain the small discontinuities in the chemical potential, as computed from the Monte Carlo simulations. Apparently, with all the aforementioned effects taken together, we get an excellent agreement with the Monte Carlo data in the entire two-phase region.  Finally, we briefly comment on the role of surface stress. In isotropic approximation, the surface energy for a spherical inclusion is given by F s = 4πγ 0 R 2 H with the surface energy γ 0 . The surface stress contribution is F β = 4πR 2 H β( θθ + φφ ) in spherical coordinates, with the scalar surface stress coefficient β [23,24]. We obtain-in the same approximation as for Equation (6)-due to the additional pressure difference between the precipitate and the matrix ∆σ rr = 2β/R H also a contribution from the elastic bulk energy ∆F el [25] where we have neglected a small destabilizing term proportional to β 2 for |β| 3K∆ε 0 R H (we have estimated from the atomistic data β to be of the order −1 Jm −2 , thus this condition is fulfilled), and also assumed c 1, since spherical precipitates appear for small concentrations only, as discussed below. Consequently, these terms have the same radius scaling as the surface energy term ∼R 2 H , and despite their partial origin from a bulk energy, they appear effectively as surface term; therefore, we treat them via a renormalized interfacial energy density γ.
With this, the parametrization of the continuum model is completed, and it is used in the following for predicting hydrogen partial pressures, volume relaxation effects, and for the prediction of phase diagrams. Nonlinear elastic MC Nonlinear elastic + surface Dilute limits Figure 10. The chemical potential of hydrogen at T = 300 K. The dots are the data from the Monte Carlo simulations. The chemical potential including nonlinear elastic effects is very close to the Monte Carlo data, and, together with interfacial effects, the agreement is excellent. The contribution from the interfacial terms is largest for low concentrations, as there the precipitates are small and therefore have the largest surface-to-volume ratio. The dotted curves near c = 0 and c = 1 are the analytical predictions for the single phase dilute limits. The H 2 partial pressure is based on an ideal gas model, which is not realistic for higher pressures and only serves for illustrational purposes.

High Concentration Limit
We can also correctly describe the limit of high hydrogen concentrations, where almost all octahedral sites are filled with H. For 1 − c 1, the system is again in a single phase state, with a dilute distribution of vacancies, thus the steep divergency of the chemical potential is due to the configurational degrees of freedom. Nevertheless, the central difference to the dilute limit c 1 is that the material is here severely under (nonlinear) stress (as we focus on fixed volume situations, where pure Ni is stress free, hence stresses are maximized for c = 1) , thus elastic effects play a major role here. In contrast to the phase separated region we need precise knowledge of the concentration dependence of the elastic constants, their nonlinearities and the eigenstrain. Consideration of all these effects leads to the good agreement between the Monte Carlo data and the continuum picture, as shown in Figure 10.

Conversion to Partial Pressures
To relate the chemical potentials to experimentally accessible quantities, it is useful to translate them to hydrogen (H 2 ) partial pressures. In the following, we use a subscript H 2 to discriminate between quantities related to the gaseous hydrogen molecules and the monatomic hydrogen in solid solution. From the ideal gas equation and the Maxwell relation we obtain by integration the chemical potential of H 2 up to an additive term µ 0 H 2 , with an arbitrary reference pressure p ref . From the equilibrium between H 2 and dissociated hydrogen follows This leads to We therefore need an expression for µ 0 H 2 (T) to link the chemical potential to a hydrogen partial pressure.
In the dilute limit c 1, the chemical potential is (see also Figure 5) In combination with Equations (12) and (13), we therefore get Sievert's law This is in agreement with experimental findings [26] for the solubility of H 2 in Ni expressed in terms of the quantities used here, with c 0 = 5.8694 · 10 −7 and p 0 = 133.3224 Pa (the prefactor 2 stems from the dissociation of a H 2 molecule in two H atoms). This allows to identify an expression for the offset potential µ 0 H 2 (T) where we have chosen p ref = p 0 /(4c 2 0 ). With these identifications, Equation (14) is the desired relation, which is shown in Figure 10. Apparently, the H 2 partial pressures needed for the higher concentrations are extremely high and not reachable in experiments under normal conditions, making computer simulations particularly useful there. In addition, we note that for such high pressures the used ideal gas description of H 2 breaks down. On the other hand, the H 2 partial pressures needed for formation of a hydride are much lower for free volume situations, which will be discussed below (see Figure 11). We can therefore conclude that volumetric constraints, or-more generally-compressive stresses can suppress the formation of hydrides, which can be responsible for the embrittlement of the metal. Nonlinear elastic MC Free volume Figure 11. Comparison of the chemical potential for fixed volume and in a stress-free situation, P = 0, T = 300 K. The volume constraint leads to a positive slope of the chemical potential in the two-phase region, thus stabilizing phase separation of given chemical potential. In contrast, two phase states are unstable for given µ H , as shown here using the continuum prediction.

Volume Relaxation
Instead of a fixed volume constraint, one may also consider a situation with free expansion, P = 0. This does not mean that elastic effects disappear, because in the two-phase region internal coherency stresses still arise. Thus, the same elastic barrier appears for the first nucleation of the precipitate, but thereafter the elastic contribution to the chemical potential decreases with concentration, see Figure 11. Since the stress effects are much lower, elastic nonlinearities are negligible here, and the curve is calculated in linear isotropic approximation. The qualitatively different behavior shows the important role of external boundary conditions. In contrast to scenarios with free expansion, the volume constraint stabilizes two-phase equilibria. This behavior can intuitively be understood as follows: When a precipitate is inserted into the matrix, coherency stresses arise, which increase the elastic energy. If the precipitate grows and finally fills the whole sample, the freely expanding material is homogeneous and therefore stress free. This implies that the elastic energy decays with the hydride volume fraction, and therefore the elastic contribution to the chemical potential appears with negative sign. Therefore, we get a hydrogen chemical potential with a negative slope in the two-phase region, which corresponds to an unstable situation. We have verified this prediction using Monte Carlo simulations for fixed and free volume at T = 300 K, see Figure 12.
We note that volume constraints can therefore stabilize precipitates. In a simulation, this is useful for probing situations with phase coexistence, which would otherwise not appear, as either the system is homogeneously in the dilute state or fully saturated with H for given chemical potential. This is expressed here by a negative slope of the chemical potential of H. The present physically motivated approach is therefore an alternative to bias potentials [12] to probe configurations with phase coexistence.

Phase Diagram
Based on the temperature dependence of F c , we can predict the entire bulk phase diagram without elastic effects via the Maxwell construction, see Figure 13. We note that this phase diagram is derived solely from T = 0 K data for the H-H interaction. To verify it also for higher temperatures, we have compared the analytical expression for the chemical potential to the single-phase data from grand-canonical Monte Carlo simulations with free volume relaxation, where elastic and interfacial effects do not appear, and find an excellent agreement (see Figure 6). We note that the purpose of the present work is to match a continuum description of Ni-H to Monte Carlo simulations based on an empirical interatomic potential. The agreement between the two theoretical approaches is excellent. Compared to experimental values, however, the description is less suitable for high temperatures, and therefore it is not accurate in this regime. The reason is-apart from potential limitations of the interatomic potential-the use of the (lattice) Monte Carlo molecular statics description, which neglects in particular phonon excitations.  Figure 11) stable two phase states with arbitrary average concentrations can be obtained in thermodynamic equilibrium. (Bottom) For free volume, the miscibility gap is inaccessible for the given chemical potential in the grand canonical simulations, and instead the equilibrium state is almost free of hydrogen or almost fully saturated, as the solubilities are low at low temperatures (see phase diagram Figure 13). The reason for this behavior (which is accompanied by a hysteresis) is the negative slope of the chemical potential in the two-phase region (see Figure 11), which implies that these configurations are unstable. All simulations are started as a hydrogen free system, and formation of the hydride therefore demands to overcome a nucleation barrier. An application that we can readily extract from the above results is the phase diagram taking into account elastic effects. In particular, we consider the situation of fixed volume and determine the two-phase region by minimization of the total free energy. For the sake of simplicity, we assume a slab geometry with the phase normal in [100] direction (see Figure 14), which is the most relevant pattern in the medium concentration regime, as discussed above in Figure 10. In contrast to other geometries, homogeneity implies the constancy of concentrations in each phase. We neglect the elastic nonlinearities, which allows for a straightforward solution of the underlying one-dimensional elastic problem. As the elastic constants were determined mainly for c = 0 and c = 1, we assume additionally that they interpolate linearly for arbitrary concentrations 0 < c < 1. This is mainly relevant for higher temperatures, where the appearing phases have larger solubilities for hydrogen or vacancies, thus they appear with intermediate concentrations. As will be discussed below, it turns out that the concentration dependence of the elastic constants is not important, and therefore this approximation is legitimate.

x=[100]
L X 0 Phase 1 Phase 2 Figure 14. Sketch of the hydride phase (light grey) forming as a slab in the nickel matrix with fixed volume. Due to translational invariance, we can assume the hydride phase to be located in the left part of the cube, starting at x = 0.
We inspect the case of a constant total volume which is stress free for pure nickel, i.e., the lattice constant for a homogeneous system is a Ni . Consequently, for the two-phase system with cubic symmetry, the only non-vanishing displacement component is u x , which is a linear function of x in each phase, u (i) x = a i x + b (the phases are enumerated by i). Hence, the only nontrivial strain component is (i) xx = a i ; stresses follow from Hooke's law with cubic symmetry, Because the concentration is homogeneous in each phase for this effectively one-dimensional situation, the linear bulk elastic equations are automatically fulfilled by the above ansatz. The coefficients a 1 , b 1 , a 2 , b 2 are determined by the boundary and matching conditions, u (2) x (x = 0) = 0, u (1) x (x 0 ) = u (2) x (x 0 ) and σ (1) xx (x 0 ) = σ (2) xx (x 0 ). The elastic energy density is in each phase From that, we can calculate the elastic energy per unit area, The total free energy is then F = µ 0 N H + F c + F H-H + F el . Additionally, the concentrations c 1 and c 2 in the two phases are related to the (fixed) average concentration c via the lever rule, We can therefore investigate F/L as function of c 1 , c 2 for given values of T, c, and the phase fractions are determined by the lever rule. The concentration domain can be restricted to c 1 < c < c 2 , and the free energy is minimized with respect to c 1 and c 2 in equilibrium. Whenever the minimum is located inside this domain, the system is in a two-phase state, whereas the minimum is located on the border c 1 = c or c 2 = c for a single phase state. By performing this minimization in the entire c, T plane one obtains the elastic phase diagram.
Already at this level we see the tremendous influence of the elastic effects in Figure 13, with an enlargement of the two-phase region towards higher temperatures and the hydrogen rich side. This unusual behavior-one would intuitively expect the suppression of phase separation since coherency stresses are energetically unfavorable-is due to the fact that the eigenstrain is a concave function of the concentration, ε 0 (c) < 0 (see Figure 7), thus volumetric deviations from Vegard's law play the dominant role here.
To investigate this peculiarity in more detail, we look at the elastic energy depending on the phase concentrations. We start with a situation which differs from the true Ni-H problem. Instead, we assume a linear dependence of the eigenstrain on concentration (see Figure 7) (Vegard's law). Here, we find that single-phase states have the lowest elastic energy in a wide region of the phase diagram and are therefore energetically favorable. Thus, phase separation is suppressed by elastic effects in this fictitious situation, leading to an increased H solubility limit, and this is shown in Figure 15. In contrast, for the true concentration dependence of the lattice constant on the H concentration in Figure 7, we obtain a reduced solubility limit of hydrogen in the nickel matrix. As mentioned above, this is a counterintuitive result, as the appearance of phase separation should usually increase the elastic energy and therefore make this state thermodynamically unfavorable, in contrast to our observation. Here, the situation is opposite: Single-phase states are more expensive from point of view of elastic energy, and the hydrogen solubility limit is therefore reduced. The explanation is that the single-phase case has a higher average eigenstrain, and therefore the mechanical volume work is higher than for the phase separated case for the anticipated fixed volume ensemble. As a result, we find that the 0 (c) dependence has a tremendous influence on the phase diagram. This effect is much more pronounced that the concentration dependence of the elastic constant. If we assume, instead of the linear interpolation between the values for Ni and Ni-H, just concentration independent constants using the values of Ni, we find only a small shift of the binodal (see Figure 15). We can therefore conclude that the concentration dependence of the lattice constant is the most relevant quantity for the influence of elastic effects on hydrogen solubility limits.

Discussion and Conclusions
Understanding the thermodynamics of phase separation across the scales is a key for deriving macroscopic material properties from fundamental microscopic descriptions. Often, such a transfer is not straightforward, but, with careful descriptions, it is possible to obtain not only macroscopic parameters but also describe energy functionals in agreement with the underlying behavior on atomic scales. We demonstrated such an approach for the Ni-H system with a special focus not only on classical thermodynamics but also strong mechanical and interfacial effects.
In detail, we have seamlessly matched the results of hybrid molecular static and Monte-Carlo simulations in a binary system with full consideration of elastic effects to a continuum description for large scale simulations. The accurate identification of configurational, interaction, nonlinear elastic and surface effects allows for a thorough understanding, which will be equally useful and applicable to other mechanisms or material systems. Presently, the description is most useful for the low temperature regime below the Debye temperature. The consideration of vibrational degrees of freedom and defect formation, as well as the consideration of more complex phase diagrams will be the subject of future investigations. It should be pointed out that the accuracy of the predictions is limited by the atomic scale description, to which the continuum energy functional is tailored. The inclusion of the effects mentioned before will presumably require the consideration of additional energy terms and variables such as defect concentration.
The determination of phase diagrams demonstrates the opportunities arising from the transfer of atomistic data to the mesoscale: The quantitative determination of the free energy allows for example to construct the entire phase diagram in Figure 13 computationally much more efficiently than using atomistic simulations which require large system sizes and long relaxation times, especially in the room temperature regime. Moreover, the description serves as an accurate basis for mesoscale simulation methods such as phase field [27].
The equality of grand potentials becomes Notice that the elastic terms do not appear in the common tangent construction in the framework of the above approximations, since the concentration is assumed to be fixed in this approximation, which holds here for low temperatures.
Next, we calculate the chemical potential of the hydrogen in the two-phase region. It is given by  Figure A1. Energy density as function of the external hydrostatic strain . With the relation in Equation (A7), we get the bulk modulus K in quadratic approximation. The atomistic data are taken from simulations at T = 0 K. For high strains deviations from the linear elastic response can be seen.

Appendix C.2. Monoclinic Distortion
The monoclinic distortion describes shear effects to get the elastic constant C 44 both for Ni and the hydride. One gets the translational movement of an atom from the original location r = r(x, y, z) to the new position r = r (x , y , z ), with a deformation parameter α. With the displacement u = r − r and the linear strain tensor ij = (∂ i u j + ∂ j u i )/2, we obtain whereas all other strains are zero and the quadratic term in α can be neglected in a linear approximation. Thus, the elastic energy is (A10) Figure A2 shows the energy density of both phases depending on the strain xy for the monoclinic distortion, and the elastic constant C 44 is obtained from the curvature at the origin. We obtain a value of C 44 = 134.3 GPa for Nickel and C 44 = 33.5 GPa for the hydride. Notice that for pure shear the quadratic energy-deformation dependence holds up to several percent of strain.

Appendix C.4. Nonlinear Elastic Coefficients
To demonstrate the extraction of the elastic nonlinearities, we inspect the mechanical response of pure Ni under bulk compression or expansion, as shown in Figure A4.
For larger strains, we see deviations from the linear response, i.e., a quadratic energy-strain relationship. For tensile stresses the material appears effectively softer, whereas it is stiffer under compression. In the strain regime shown in the plot, the elastic response is well described by in the spirit of Equation (8). The coefficients d K i are given in Table 2. For pure shear during a monoclinic transformation, the energy is well described by a quadratic strain dependence, as shown before in Figure A2, and the same holds for the volume preserving orthorhombic transformation (see Figure A3). This demonstrates that the nonlinearities depend only on the trace of the strain tensor. We have performed monoclinic and orthorhombic transformations at different lattice units (i.e., bulk compression or expansion), to extract the coefficients in the spirit of Equation (8). This is exemplarily shown in Figure A5 for the monoclinic deformation of pure Ni. This allows to extract the elastic constants, here e.g., C 44 as function of the trace of the strain tensor. This dependence is shown in Figure A6. The resulting data for both phases are then described by Equation (8) with the coefficients in Table 2. The corresponding plot for the orthorhombic deformation is shown in Figure A7.  Figure A5. Monoclinic deformation of pure nickel for different bulk lattice constants. The fits (curves) are parabolic, and the data (points) are well described without further nonlinear corrections, in agreement with the statement that corrections depend only on tr . The latter is contained in the lattice constant dependence of the curvatures, leading to strain-dependent elastic constants C 44 (tr ). The energy of the bulk compression is subtracted in the plots. These are altogether three conditions to determine the constants a d , b d and a h . The elastic energy densities (per unit volume) are The total free energy is which is a lengthy expression. Using the conversions with the lattice unit a Ni , we can calculate the elastic contribution to the chemical potential In the general case, it becomes In particular for equal elastic constants (λ d = λ h = λ, G h = G d = G) µ el = a 3 Ni (∆ε 0 ) 2 (3λ + 2G)(2G + c(3λ + 2G)) (λ + 2G)N 0 . (A30)

Appendix E. Interfacial Effects
The shape of the precipitate forming in the two-phase region is mainly determined by interfacial energy. For isotropic surface energy, a spherical inclusion minimizes the interfacial energy for small radii R H L, with L being the edge length of the enclosing box with periodic boundary conditions. However, for R H ∼ L, the formation of other interface contours is energetically favorable, since we use periodic boundary conditions in the atomistic and FEM simulations.
To find the different ranges and the limiting concentrations at which the forms change, we have to compute the minimum surface energy. The surface energy is defined as where S defines the size of the interfaces and γ the surface energy coefficient. First, we look for the spherical inclusion of the hydride with the radius R H . In this case, the surface energy is The concentration (volume fraction) is given as so the surface energy is F s (c) = γL 2 (36c 2 π) 1/3 . (A34) The surface energy for a hydride tube with radius r and length L is with the respective concentration we obtain F s (c) = 2γL 2 (cπ) 1/2 .
In the case of slabs, the surface energy is independent of the concentration c: The cases with c > 1/2 are analogous, with the role of the matrix and the hydride being exchanged. Figure A8 shows the surface energy depending on the concentration for the respective case. Up to c = 4π/81 ≈ 0.16, it should be a spherical inclusion of the hydride in the nickel matrix. For 4π/81 ≤ c ≤ 1/π, the inclusion of the hydride should be a tube followed by plates up to c = 1 − 1/π. The analytical results predict a symmetric behavior for concentrations c > 0.5, so for concentrations 1 − 1/π ≤ c ≤ 1 − 4π/81 there should be a tube and for c > 1 − 4π/81 we expect a nickel sphere. With the calculated surface energy F s , we compute now the chemical potential with µ s = dF s /dN H . Then, we have to rewrite the surface energy in an expression which depends on the number of hydrogen atoms N H . For the first case, the hydride sphere inclusion, the radius can be expressed as