New Methods for Calculating the Free Energy of Charged Defects in Solid Electrolytes

A methodology for calculating the contribution of charged defects to the configurational free energy of an ionic crystal is introduced. The temperature-independent Wang-Landau Monte Carlo technique is applied to a simple model of a solid electrolyte, consisting of charged positive and negative defects on a lattice. The electrostatic energy is computed on lattices with periodic boundary conditions, and used to calculate the density of states and statistical-thermodynamic potentials of this system. The free energy as a function of defect concentration and temperature is accurately described by a regular solution model up to concentrations of 10% of defects, well beyond the range described by the ideal solution theory. The approach, supplemented by short-ranged terms in the energy, is proposed as an alternative to free-energy methods that require a number of simulations to be carried out over a range of temperatures.


I. INTRODUCTION
The presence of charged defects gives rise to useful electrolytic behaviour in many materials of practical interest. In solid electrolytes, such as yttria-stabilized zirconia (YSZ), vacancies are present in concentrations well above 1%. (This is in contrast to simple oxides where the vacancy concentrations are usually very low.) This high concentration of (locally charged) vacancies is essential in making these materials useful in applications; lattice diffusion of highly concentrated charged vacancies is primarily responsible for the useful electrolytic behaviour in such materials [1,2].
More generally, the behaviour of solid electrolytes of this type is determined by gradients in electrochemical potential, which provide the driving forces for point-defect diffusion. The chemical potentials in turn depend on the free energy as a function of defect concentrations.
Accordingly, for modelling such systems, e.g., using phase-field methods, or other implementations of irreversible thermodynamics, a detailed understanding of the thermodynamic behaviour of these defects is highly desirable.
A major challenge to the modelling of the thermodynamics of solid-electrolyte systems is the long-ranged nature of the Coulombic defect-defect interactions. For example, the longranged character of the Coulombic interactions significantly affects defect mobilities. [3,4] However, It is unclear how best to incorporate these interactions in an empirical description of the free energy as a function of concentrations, such as that used in the Calphad methodology [5,6]. We propose an approach that overcomes this challenge through the use of an efficient atomistic-simulation method to characterise the thermodynamic functions.
We use an implementation of Wang-Landau (WL) [7,8] sampling, a biased Monte Carlo method, in a lattice-based simulation to obtain the Helmholtz free energy of the system at any temperature as a function of concentration. This free energy can then be used to parameterise macroscopic thermodynamic approaches, such as the subregular solution (see, e.g., [6]) with Redlich-Kister [9] polynomials.
We choose the WL method as it is a sampling method that produces a temperatureindependent density of states for the system within a single simulation. This allows the partition function, and hence all the statistical-thermodynamic properties of the system, to be evaluated as functions of temperature using simple quadratures.
This paper is organised in as follows. Our model system based on YSZ is described in section II. We outline our methodology, including our implementation of the WL method in section III. Section IV reports the results of applying this methodology to the model system described in section, most importantly we describe and analyse our calculations of the Helmholtz free energy of the system as a function of temperature and defect concentration. This free energy is compared to that obtained using the ideal-solution theory (see e.g., [10]), and to that from Debye-Hückel theory [11]. As a proof of principle for future applications, we show that the excess free energy obtained from simulation is accurately represented as a function of defect concentration by low order polynomials with temperature dependent coefficients, in the manner familiar to practitioners of phase diagram assessment in CALPHAD [5,6]. Finally, our conclusions are presented in section V.

II. THE MODEL
The high-temperature fluorite structure of YSZ is chosen as the basis for our simple model as a prototypical example of a material in which charged defects not only exist in high concentrations, but also play a principal role in its use as an ionic conductor. When zirconia (ZO 2 ) is doped with yttria (Y 2 O 3 ) the yttrium ions (Y 3+ ) replace zirconium ions (Zr 4+ ) on a face-centred cubic (FCC) lattice. As there is a difference in the net charge of these ions, charge neutrality in the system is preserved by the creation of locally-charged oxygen vacancies on a simple-cubic (SC) sub-lattice nested within the FCC lattice. Substitutional yttrium ions and oxygen vacancies are assumed to carry local excess charges of -1 and +2 respectively with respect to the perfect lattice, so to preserve overall charge neutrality the ration of oxygen vacancies to yttrium ions should be 1:2. The high concentration and mobility of these +2 charged oxygen vacancies gives YSZ its useful electrolytic properties.
We assume that the yttrium substitutional ions and oxygen vacancies behave as −1 and +2 point charges and Coulomb's law is treated with periodic boundary conditions by Ewald summation (see, e.g., [12]). The volume, the number of lattice sites and the numbers of charges in each simulation of the electrolyte remains constant. The only length scale in this model system is the lattice parameter. For convenience we use the physical value of 0.51 nm, appropriate to YSZ [13]. In the present case its role, along with a dielectric constant of 30 (relative to vacuum) [14], is to scale the energies; this could equally well be done after the simulation, however, it is convenient to generate energies and temperatures directly in the familiar units. Moreover, when a more-realistic potential is used to augment the Coulomb potential, in order to take into account the short-ranged interactions, the correct choice of lattice parameter is essential. Unless otherwise stated, all simulations are carried out within supercells of 4 × 4 × 4 fluorite unit cells with periodic boundary conditions one unit cell of which is depicted in Figure 1.

III. METHODOLOGY
In classical statistical mechanics we are interested in evaluating configurational integrals or partition functions Z of the form These may be written more concisely as where β = (1/k B T ), k B is Boltzmann's constant, T is the absolute temperature, E(x) is the potential energy of a microstate at x, which represents a point in the 3N-dimensional space of possible atomic positions. In a complete partition function the kinetic energy is added to the potential, with integration extending also over the 3N momental; this is done analytically to give a familiar temperature-dependent de Broglie prefactor in front of the above integral. In our model, however, the kinetic energy is zero, and will not concern us here.
All the thermodynamic functions can be expressed either directly as Boltzmann averages, or in terms of derivatives of Z with respect to β. For example the Helmholtz free energy (F ) of the system is and the internal energy (U) is Evaluating the configurational integral (1) or its derivatives by a brute-force method is usually highly inefficient, because the variation in energy of a small number of the possible microstates at low energy may have a disproportionate impact on the results. Moresophisticated methods are therefore highly desirable. Traditionally, Metropolis Monte Carlo (MMC) [15] has been the method of choice for evaluating thermodynamic functions of a given system. However, MMC is carried out in the canonical ensemble and, as a consequence, a simulation is required for every temperature one wishes to examine; consequently, the approach can be computationally expensive. The formulation of the method also means that evaluating properties that are not characterized as Boltzmann distributions, notably F or the entropy S, is more challenging, and it is normally done by integrating a function that can be calculated as a configurational average, which requires a whole series of long simulations to obtain the integrand at each point on the path.
The problems associated with MMC have meant that in recent years there has been a renewed interest in developing Monte Carlo methods that allow the full thermodynamics of a system to be calculated from a single simulation [16][17][18]. One of the more-popular of these methods is an approach developed by Wang and Landau [7,8], which enables the evaluation of the configurational integral through the determination of the density of states of the system. This method is widely used, has been generalized for off-lattice systems [19], and has now been applied to systems ranging from metals [20] to proteins [21]. It has also been incorporated into molecular-dynamics simulations [22].

A. Wang-Landau Sampling
Within the WL formalism one proceeds from a reformulation of the configurational integral, equation (1), in terms of the normalised density of states of the system g(E) at energy E: where is the normalised density of states, Ω is the total number of microstates of the system, and δ is the Dirac delta function. It is apparent from Equation (6) that given g(E) (and knowledge of Ω) one could calculate the configurational integrals and, thus, the thermodynamics of the system at any temperature by simple quadrature. Unfortunately, in most cases of interest, one has no a priori knowledge of g(E). The aim of the WL method is to sample the system in such a way that g(E) can be accumulated as a single simulation progresses. As in conventional Monte-Carlo methods, a 'simulation' refers to the process of evolving the system from microstate to microstate in a series of random steps that are accepted or rejected according to specific criteria. Thus it does not imply any attempt to simulate a physical trajectories in phase space. The term 'walker' is used to refer to a member of an ensemble of microstates that evolve simultaneously by steps of the walkers during the simulation.
Following [23], one first seeks to obtain a Metropolis-like acceptance criterion from the detailed-balance condition requiring that, at equilibrium, the net flux of random walkers between the set of microstates with energy E i and those with energy E i+1 must be zero.
This condition can be expressed in terms of the density of walkers, n (E i ), the probability, , of attempting to make a move from an initial microstate among those with energy E i to a final microstate among those with energy E i+1 , and the probability, P A E i →E i+1 , that this move will be accepted: We observe that if the system is sampled evenly in energy space, i.e., if one obtains Conversely, if we can define an acceptance criterion that will ensure (10) is satisfied, then from (8) the sample population n(E i ) must be uniformly distributed in energy. This is rather straightforward to achieve, as follows.
If the microstate to which a random walker will attempt to move is chosen in an unbiased manner, the probability of attempting to move to any given microstate with energy E is equal to the normalised density of states at that energy (i.e., By enforcing a Metropolis-style sampling criterion, the probability of accepting a move to a given microstate will be given by where w E j is a weight associated with each energy. (In conventional MMC these weights are set equal to the Boltzmann factor). Substituting (12) into Equation (11) one obtains a condition that the weights must satisfy: This holds for all energies only if This result implies that Monte Carlo moves based on the acceptance criterion, will satisfy detailed balance in the form required, implying that the walkers are distributed evenly in energy. The strategy of WL is to ensure that this acceptance criterion is approached as the system evolves, so that the distribution of walkers will converge to a constant value of n(E i ) for all i.
In the WL method one proceeds by systematically improving an initial estimate of the density of states (g 0 (E)) (usually this is chosen such that g 0 = 1). A random walk is taken according to (14) and, after every attempted move, the estimate of the density of states for the current energy of the system is increased by some multiplicative factor f . In practice one maintains a record of ln g (E) and increases it according to where j indicates a point along the random walk. In theory this process should be repeated until the walkers have visited all regions of energy space "equally", i.e., Equation (9) is satisfied. This, in turn, means that Equation (10) holds and, thereby, satisfying Equation (13) implies that the estimate of the density of states has converged to the true value.
However, the probability of every region of energy space being visited an equal number of times is extremely small so, in practice, simulations would continue indefinitely if this convergence criterion were applied. Instead, in the original implementation of WL, energy space is divided into windows, and a histogram is updated every time a window is visited by a walker. The simulation is taken to be converged when this histogram satisfies some (user-defined) flatness criterion. A proof of convergence of this procedure has been given in [24].
It was also suggested in [24] that, as the statistical fluctuations in the histogram are proportional to 1 √ ln f , a more rigourous convergence criterion than this arbitrary, used-defined flatness criterion should be applied. They proposed that convergence should be deemed adequate when each histogram bin of energy space has been visited at least 1 √ ln f times. This is the criterion for convergence we have adopted in the present work.
Knowledge of the statistical fluctuations in the histogram, which is what we want to approach the density of states, also motivates the choice of value for the update factor f .
In principle one would like to use a value of f which is as small as possible to minimise this error. However, such a choice renders the simulation very computationally expensive, since the walkers take more cycles to explore phase space. On the other hand, since the value of f implicitly limits the accuracy of the convergence, a large value may not allow the density of states to converge with sufficient accuracy. In order to address this issue one repeats the WL procedure a number of times, using the density of states from the previous iteration obtained using a larger value of f , as the starting point for the next. By systematically reducing the update factor for each WL iteration one successively gains a better approximation of the density of states of the system whilst minimising the computational cost.
In practice one starts with f = e (i.e., ln f = 1) and, once the convergence criterion has been met, this value of f is reduced by an order of magnitude and the process is repeated, reducing the update factor in the same manner each time the convergence criterion is met. This is carried out until the update factor has a minimal impact on the previous estimate of the density of states. It is at this point that one assumes an accurate estimate of the density of states has been obtained.
As mentioned in [7], rather than calculating the exact density of states of the system, one in fact calculates a relative density of states; in practice one calculates the density of states of the system multiplied by a simulation-dependent constant. However, one can readily obtain a normalised density of states from a WL simulation. Thus the only remaining unknown in Equation (6) is the total number of microstates of the system (Ω). As all the simulations to be carried out will take place on a lattice, Ω can be directly enumerated through combinatorics (see e.g., [10]). Thus, having calculated g (E) and Ω, one can evaluate Equation (6) and hence all required thermodynamic functions through simple quadrature.

IV. RESULTS FOR LATTICE MODELS OF YSZ
We begin by illustrating, in Figure 2, the evolution of the estimate of g (E) after each successive WL iteration for a 2: Once calculated, g (E) can be analysed by examining representative microstates over the full energy range. Unsurprisingly, the very lowest-energy microstates of our system correspond to trimers of one oxygen vacancy and two nearest-neighbour yttrium substitutions. This is illustrated in Figure 3 for the simple system containing just three defects. By choosing a system with very few defects, the individual features in g (E) are preserved and easier to discern, whereas they are "smoothed out" when examining a system containing many defects. Knowledge of g (E) allows one to calculate the full thermodynamics of the system by evaluating the configurational integral, equation (6), thereby obtaining the configurational internal energy, equation (5) and the Helmholtz free energy, equation (3) over the entire temperature range. In Figure 4 for the internal energy of the three-defect system we include data points corresponding to internal energies evaluated using traditional MMC simulations in addition to the WL results; good agreement between the two methods is apparent. Also illustrated in the figure are the changes that can be attributed to the three different regimes the system passes through as the temperature is increased. As one might expect from the analysis of representative microstates at different regions of g (E), with increasing energy the system evolves from a fully associated to a fully dissociated regime. For all sizes of system, structures composed of charge-neutral trimers are stable at lower temperatures.
These dissociate successively into dimers and monomers as the temperature increases and entropic effects begin to dominate.
One can now calculate the free energy of the model YSZ system as a function of temperature and defect density. In Figure 5 the density dependence is illustrated for isotherms at  Figure 6 the temperature dependence of the Helmholtz free energy for the 1% defect isochore is displayed. The slope of this curve is the negative entropy, and its change with temperature is indicative of the dissociation of the charge-neutral trimers, discussed in relation to Figure 4.
It is interesting to compare the relative efficiency of the WL and MMC routines in obtaining convergence of the internal energy. In Table I we provide the number of energy evaluations (i.e., calculations of the total energy of the system) required, using each technique, to obtain convergence of the internal energy to an accuracy of 0.001eV per defect. It is worth reiterating here that each simulation using the MMC technique relates to a single, fixed temperature, whereas the WL simulation technique provides a description over the full temperature range. One can observe that the relative efficiency of the two techniques is concentration dependent as one would expect. For low concentrations of defects the WL method requires fewer energy evaluations than MMC, so that WL obtains the thermodynamics of the system at all temperatures at less computational cost than a single MMC simulation at fixed T . At higher concentrations the advantage is less spectacular; although with 10% yttria a single WL simulation makes about ten times more energy evaluations than a single MMC simulation, one may require many more than ten MMC simulations to obtain the desired information. Furthermore, a number of strategies for improving the performance of the WL method exist [24][25][26]; a discussion of the merits of these strategies is beyond the scope of our current work, but a more rigorous analysis of the scaling behaviour of the methodology will be outlined in a forthcoming publication.

A. Comparison with an ideal solution
One can readily observe from Figure 7 that, as the temperature increases, the behaviour of the system tends towards that of an ideal solution [10], which one can calculate using simple combinatorics for this lattice system, in which the only contribution to the free energy is the configurational entropy. This model thus takes the form where n i is the total number of a given charged species i, N i are the total number of lattice sites available to charged species i and N c is the total number of charged species in the system.  Figure 7 is what one would expect, given the behaviour already illustrated in Figures 4 and 6 wherein, with increasing T , the system undergoes a transition from a regime in which the energetics (Coulombic interactions) dominate to one that is predominantly characterised by its entropy. It is, however, important to point out that this behaviour is present only in the low-density regime. This also is to be expected as the ideal solution model is, by construction, exact in the limit of infinite dilution.

The behaviour illustrated in
Since at the lower temperatures dimers and trimers are increasingly prominent, we can expect large departures from ideal behaviour. To understand this in more detail it is helpful to study the behavour of the entropy, evaluated as (U − F )/T , which is plotted versus concentration at four different temperatures in Figure 8. For comparison on this figure we have also plotted the entropy for two extremes of ideal solution, in which the charges are 100% associated as neutral trimers, which is the low temperature limit, and 100% dissociated,

B. Comparison with Debye-Hückel theory
For systems characterised by a low density of charged defects an alternative approximation to the free energy can be obtained with Debye-Hückel (DH) theory [11,27]. The DH contribution to the free energy is given by (17) and the internal energy by where V is the volume of the system. The state-dependent parameter κ is the reciprocal of the Debye length, and is defined as as before N c is the number of charged species in the system, n i is the total number of charged species i, of charge z i ; ǫ is the dielectric constant. Other approaches that should be more accurate than the ideal solution mode could be considered, such as the mean spherical approximation (MSA) [28][29][30], but we do not consider these in the present work.
In Figure 9 we compare the DH free energy with that computed for our system with the WL approach. There is reasonable agreement at higher temperatures, where both the DH description and the WL results tend to the ideal limit, however at low temperatures the overall agreement between DH theory and the WL simulations, which we consider to be the 'exact' results for this model system, becomes worse. On further analysis it becomes evident that this disagreement is due to the poor description of the internal energy of the system with the DH approach; this is illustrated in Figure 10. The poor quality of the DH description can be understood in terms of the relatively high density of defects in our system which have a far from random distribution as they tend to lower the internal energy by forming dimers and trimers; this undermines the assumptions made in the formulation of DH theory, a mean-field theory that is strictly valid only for low charge densities.
In contrast with the internal energy, one can observe from Figure 11 that the entropic contribution obtained with the DH theory is in rather better agreement with the WL values, correcting some of the error made by the ideal approximation. Thus given an MMC calculation of an accurate internal energy at a single temperature one could make a rapid and reasonably good estimate of the free energy by using the DH value of entropy, as illustrated in Figure (12).

C. Regular solution Models
It is useful to explore how an analytic description of the free energy can be developed beyond the ideal solution model, which is only accurate at very low densities or high temperatures. In thermodynamic modelling of solutions, as in the CALPHAD approach [5,6] where x i and x j are the concentrations of two components and the L ij are adjustable parameters; in our implementation these are treated as temperature dependent. One can observe from Figure 13 that a simple one parameter RS model of this type can be employed to provide a very good description of the thermodynamics of the system.
Although for our simple system, a good representation of the thermodynamic properties can be achieved with the RS model, for more-complicated systems it may be beneficial to go a step further and use empirical extensions such as the Redlich-Kister (RK) [9] model, in which the parameters of the regular-solution model are taken to depend on the concentration of the species, i.e., where A ij; ν are adjustable parameters. This would enable accurate simulations to be fitted analytically with arbitrary precision.

V. DISCUSSION AND CONCLUSIONS
We have proposed a methodology for the description of the thermodynamics of charged defects in crystals, which enables the free energy of a model of interacting point defects to be calculated accurately even at defect concentrations up to about 10%, as exhibited by complex oxides of practical interest, such as yttria stabilized zirconia. This is far above the range of validity of the ideal solution model that is commonly used in discussions of point defect equilibrium in ceramic materials. In our approach, the temperature-independent Wang-Landau (WL) Monte Carlo method is used to calculate the free energy of the system at the atomistic level of an electrolyte lattice. As a test bed for this approach, we apply it to a very simple model of yttria stabilized zirconia, in which only the Coulomb interactions, dominant at all but nearest neighbour sites, are considered. The approach is much more efficient than conventional Metropolis Monte Carlo, since following a single, temperature-independent simulation one can obtain the Helmholtz free energy and entropy at all temperatures of interest by simple quadrature over the density of states.
We have shown that the free energy tends to ideal solution behaviour in the hightemperature, low-density limit, as expected. A comparison of our WL calculations with the free energy obtained analytically with Debye-Hückel theory shows that the latter provides a reasonable description of the non-ideal entropy, but a poor description of the internal energy.
Our resulting free energies as a function of temperature and defect concentration can be accurately fitted with a regular solution model using one temperature-dependent parameter. This suggests a promising route to providing input from physical models into the data needed for the calculation of phase diagrams. A more sophisticated model of the crystal and its defects could be parameterised from, e.g., DFT calculations, incorporating correctly the long-ranged Coulomb interaction, together with short ranged potentials or cluster models, and even the elastic interactions between defects. A free energy calculation using the methodology introduced here would then correctly include all the correlation between defect positions, and the energy and entropy therefrom that is neglected in the dilute or ideal solution model.