A comparison of RedlichKister polynomial and cubic spline representations of the chemical potential in phase field computations
Abstract
Free energies play a central role in many descriptions of equilibrium and nonequilibrium properties of solids. Continuum partial differential equations (PDEs) of atomic transport, phase transformations and mechanics often rely on first and second derivatives of a free energy function. The stability, accuracy and robustness of numerical methods to solve these PDEs are sensitive to the particular functional representations of the free energy. In this communication we investigate the influence of different representations of thermodynamic data on phase field computations of diffusion and twophase reactions in the solid state. Firstprinciples statistical mechanics methods were used to generate realistic free energy data for HCP titanium with interstitially dissolved oxygen. While RedlichKister polynomials have formed the mainstay of thermodynamic descriptions of multicomponent solids, they require high order terms to fit oscillations in chemical potentials around phase transitions. Here, we demonstrate that high fidelity fits to rapidly fluctuating free energy functions are obtained with spline functions. Spline functions that are many degrees lower than RedlichKister polynomials provide equal or superior fits to chemical potential data and, when used in phase field computations, result in solution times approaching an order of magnitude speed up relative to the use of RedlichKister polynomials.
keywords:
Free energy; Spinodal decomposition; Phase transformation1 Introduction
A free energy potential obtained by applying appropriate Legendre transforms to the internal energy encapsulates all the thermodynamic information that can be known about a bulk solid. Its first derivatives with respect to experimentally controlled state variables yield the equilibrium values of conjugate state variables that are not controlled experimentally. Its second derivatives are related to response functions such as heat capacities, compressibilities, elastic moduli etc. A free energy potential also codifies information about phase stability, and its rate of decrease points the direction of nonequilibrium processes in kinetic theories of phase transformations and microstructure evolution. The CALPHAD approach from its inception recognized the central role that free energies should play to ensure that experimental and calculated data is collected and organized in a thermodynamically selfconsistent way (Hillert, 2007; Kaufman and Agren, 2014).
The Gibbs and Helmholtz free energies are especially important in the study of multicomponent solids. The Gibbs free energy is the potential used to calculate temperature composition phase diagrams at ambient pressure when assuming incoherent multiphase coexistence (Hillert, 2007). The Helmholtz free energy is more appropriate when analyzing coherent multiphase equilibrium, where the pressure is no longer uniform throughout the solid (Voorhees and Johnson, 2004).
A direct experimental measurement of a free energy is often not possible. Instead, free energies must be determined indirectly by integrating over measurable state variables. Chemical potentials, for example, can be determined electrochemically or by measurements of partial pressures. The integration of this data with respect to composition then yields a free energy (DeHoff, 1993). Free energies can also be calculated using firstprinciples statistical mechanics approaches that rely on effective Hamiltonians to extrapolate computationally expensive electronic structure calculations within Monte Carlo simulations (Sanchez et al., 1984; Fontaine, 1994; Laks et al., 1992; Ceder, 1993; Van der Ven et al., 1998; van de Walle and Ceder, 2002; van de Walle and Asta, 2002; Arroyave et al., 2005; Van der Ven et al., 2010a; Thomas and Ven, 2013; Puchala and Van der Ven, 2013; Chen et al., 2015; Natarajan et al., 2016). While such approaches are rarely quantitatively accurate, they yield free energy descriptions that are physics based and that rigorously account for vibrational and configurational entropy.
Independent of how a free energy description has been obtained there is a need to represent it mathematically. The CALPHAD approach maps the composition dependence of free energies on a polynomial expansion. A commonly used polynomial expansion is that introduced by Redlich and Kister (Redlich and Kister, 1948). While other polynomial expansions can also be used, “mathematically, the choice of basis (for a finitedimensional space) makes no difference” according to Dahlquist and Björck (Dahlquist and Björck, 2008). It is known that sampling data at Chebyshev points significantly improves a polynomial interpolation, where the function passes through every data point. However, fitting to measured or calculated free energies generally involves many data points, making a least squares fit more appropriate than a polynomial interpolation. Chebyshev points are, therefore, not necessary (Dahlquist and Björck, 2008). Because of the mathematical equivalence of the various polynomial expansions (RedlichKister expansion, Legendre polynomial series, simple power series, etc.) (Tomiska, 1984; Ouerfelli et al., 2010), the effectiveness of a polynomial expansion can be assessed by considering one basis set. Even so, Dahlquist and Björck (2008) point out that some functions are “not at all suited for approximation by one polynomial over the entire interval. One would get a much better result using approximation with piecewise polynomials.”
In this study, we compare the use of the RedlichKister polynomials with cubic splines (piecewise cubic polynomials with global continuity; i.e. continuous second derivatives) in fitting free energy data. This is done in the context of phase field modeling using the CahnHilliard equation (Cahn and Hilliard, 1958). One physical phenomenon that this model captures is spinodal decomposition, where a material separates into two distinct phases. Spinodal decomposition arises when the free energy is concave with respect to composition resulting in a negative thermodynamic factor that causes uphill diffusion. Capturing this physics in a phase field model requires an accurate representation of the free energy and its higher order derivatives. We find that there are cases where even loworder splines are much more effective at representing the physics of the problem than are global polynomials, especially within the spinodal regions. We also see that the high polynomial degree sometimes required by the RedlichKister expansions significantly increases computation time.
The thermodynamics of our model, the titaniumoxygen system is discussed in Section 2. A detailed treatment of the different functional representations of the chemical potential and the free energy follows in Section 3, supported by fits to a data set generated from first principles. Phase field computations is the subject of Section 4, where we focus on the consequences of the different chemical potential and free energy representations in terms of the physics and numerical performance. We place our results in perspective and make concluding remarks in Section 5.
2 Model system: oxygen dissolved in HCP Ti
We use the TiO binary as a model system to explore different representations of the composition dependence of free energies. In contrast to many common metals, HCP titanium is capable of dissolving a high concentration of oxygen before other Ti suboxides having very different crystal structures form (Cancarevic et al., 2007; Okamoto, 2011). The dissolved oxygen fills octahedral interstitial sites of the HCP Ti crystal structure to form TiO where measures the fraction of filled octahedral sites (Burton and van de Walle, 2012). The oxygen solubility limit in HCP Ti is as high as =1/2 (Cancarevic et al., 2007; Okamoto, 2011).
Since oxygen dissolves interestitially in HCP Ti, its chemical potential is related to the slope of the free energy according to
(1) 
provided that the total free energy of the solid, , is normalized by the number of Ti atoms, , with and . A measurement of the oxygen chemical potential as a function of oxygen concentration , enables a determination of the free energy through integration of Eq. (1).
An important source of entropy in HCP TiO arises from configurational disorder among oxygen and vacant interstitial sites. A firstprinciples statistical mechanics approach based on the cluster expansion is well suited to predict the composition dependence of its free energy(Burton and van de Walle, 2012). We used the CASM software package (Van der Ven et al., 2010a; Thomas and Ven, 2013; Puchala and Van der Ven, 2013) to calculate a firstprinciples phase diagram with a cluster expansion (Sanchez et al., 1984; Fontaine, 1994) parameterized with density functional theory (DFT) energies and Monte Carlo simulations. The VASP planewave DFT code (Kresse and Furthmüller, 1996; Kresse and Joubert, 1999) was used to calculate the energies of many different oxygenvacancy orderings within HCP Ti. Details of these calculations and the parameterization of the cluster expansion will be published elsewhere.
The calculated phase diagram of Figure 1 shows that ordered phases are stable at low temperature at stoichiometric compositions of =1/6, =1/3 and =1/2. The ordered phases are shown in Figure 2. These ordered phases disorder upon heating to form a solid solution at high temperature. The octahedral interstitial sites of HCP form twodimensional triangular lattices that are stacked directly on top of each other along the HCP axis. The ordered phases at =1/6 and =1/3 have a supercell on the triangular lattices of octahedral sites within basal planes of HCP as illustrated in Figure 2. Both ordered phases are staged in the sense that every other layer of interstitial octahedral sites is empty. Oxygen in the =1/2 ordered phase arrange in a zigzag pattern separated by a zigzag pattern of vacant sites. The ordered phase at =1/6 disorders around 1000K when at its stoichiometric composition while the ordered phase at =1/2 disorders around 1100K. The ordered phases have a wide stability range along the composition axis due to their ability to tolerate antisite defects. The ordered phase at =1/6, for example, is predicted to remain thermodynamically stable to concentrations as high as =0.25 at 600 K. The excess oxygen in the =1/6 ordered phase is accommodated in a disordered fashion on the vacant sites within the filled layers. A second order phase transition separates the stability domain of the =1/6 and =1/3 orderings. The calculated phase diagram of Figure 1 is qualitatively very similar to that predicted by Burton and van de Walle (Burton and van de Walle, 2012) using a similar approach, although lower orderdisorder transition temperatures are predicted in Figure 1. The phase diagram in Figure 1 also shows the stability range of the BCC based TiO phase observed experimentally (Cancarevic et al., 2007; Okamoto, 2011).
Figure 3 shows the calculated oxygen chemical potentials and free energies for HCP TiO at 1800 K and 800 K. The free energies were calculated by integrating the chemical potential versus the oxygen concentration (Dalton et al., 2012). At 1800 K, the disordered solid solution is stable between =0 and =1/2 as reflected by the sloping oxygen chemical potential and the convex nature of the free energy at this temperature. At 800 K, both the =1/6 and =1/2 ordered phases are stable with twophase regions separating the ordered phases from the solid solution. The free energy within the twophase regions is continuous due to the group/subgroup symmetry relation between the ordered phases and the disordered solid solution. We calculated the chemical potential within the twophase region using variance constrained Monte Carlo simulations (Sadigh and Erhart, 2012), which permits access to compositions where the free energy is concave and the solid is susceptible to decomposition.
The ordered phases at =1/6 and =1/2 have a lower symmetry than the disordered solid solution TiO. Order parameters, , therefore exist that can distinguish the ordered phases from the disordered solid solution. Usually order parameters are defined such that they are zero in the disordered state () and have finite values in the ordered phase (i.e. for ordered phase ). The free energy of the solid can therefore be expressed not only as a function of concentration, but also as a function of order parameters . When a solid solution is thermodynamically stable, this free energy will exibit a minimum at . Likewise, at concentrations where a particular ordered phase is stable, the free energy will have a minimum in the vicinity of . The free energies depicted in Figure 3 can be viewed as the minimum of with respect to order parameters at each concentration , i.e.
(2) 
3 Different free energy representations
We compare the ability of a truncated RedlichKister series and a spline fit to faithfully represent the free energies of TiO. We first consider the high temperature free energy curve corresponding to a solid solution. The low temperature free energy poses more challenges due to the presence of several twophase regions where the free energy is concave.
3.1 Methodology
Following standard CALPHAD practice, we express the free energy as
(3) 
where the first term corresponds to an ideal solution entropy and is an excess free energy representing a deviation from thermodynamic ideality. As usual is Boltzmann’s constant and natural logarithms are used. The inclusion of an ideal solution entropy term is especially useful in describing the free energy in the dilute limits (i.e. and ) where alloys behave as ideal solutions and where the excess free energy goes to zero. With the above free energy expression, the oxygen chemical potential in TiO then becomes
(4) 
where is the derivative with respect to of . For a substitutional binary AB alloy, the chemical potential in Eq (17) corresponds to the difference in chemical potentials between A and B (i.e. if measures the concentration of B). The ideal solution term captures the logarthmic divergences in the chemical potential in the dilute limits ( and ).
Below we compare the ability of a RedlichKister polynomial series and splines in fitting calculated values for the excess free energy . Since the free energy is often found by first measuring or computing the chemical potential and then integrating, we fit to the difference between the chemical potential data and the logarithmic term .
3.1.1 RedlichKister polynomial series
We consider a RK polynomial expansion for the function of the following form:
(5) 
Then the function used in the curve fitting of the chemical potential has the form:
(6) 
The coefficients can be found by a least squares method.
We used the chemical potential data for the titanium oxide system for compositions from to . To avoid the noise where the data is steepest near and , we only fit to the data between 0.001 and 0.499.
3.1.2 Splines
A spline is a piecewise polynomial with a specified order of continuity at the subdomain junctions. The endpoints of each subdomain are called knots or breaks. Quoting verbatim from Dahlquist and Björck (2008), a more formal definition is as follows:
Definition 1
A spline function of order (degree ), on a grid
of distinct knots is a real function with the following properties:

For ,, is a polynomial of degree .

For , is a piecewise constant function. For , and its first derivatives are continuous on , i.e.
Note that the order of the spline is equal to the degree of the polynomial plus one. For example, cubic splines (order = 4) are piecewise cubic polynomials that are continuous across all knots (see Figure 4). Splines are able to capture local features due to the local nature of piecewise polynomials. The continuity conditions at the knots maintain a specified global continuity across the entire function. These local and global characteristics of splines make them an important and often ideal tool in function interpolation and fitting.
There are multiple ways of representing splines. The simplest form is as a standard piecewise polynomial, where a standard polynomial is given for each subdomain. There is nothing in this structure that guarantees that the function is, in fact, a spline; the polynomial coefficients must be chosen to ensure the necessary continuity across knots. Splines may also be written as a linear combination of Bsplines (short for basis splines). We will denote a Bspline of order at knot by , and we let be the location of knot . Additional knots may be added at the endpoints, called exterior knots. The original set of knots will be called interior knots. The Bspline is strictly positive on the interval and zero outside the interval. Given interior knots, Bsplines have the summation property , for (Dahlquist and Björck, 2008). Bsplines can be represented using the following Coxde Boor recurrence relation (de Boor, 1972):
(7)  
(8) 
We can represent a spline function of order with interior knots as a linear combination of Bsplines:
(9) 
where are the appropriate Bspline coefficients. The piecewise polynomial representation is
(10)  
where are the coefficients for the polynomial in subdomain .
In the case of spline interpolation, each given data point becomes a knot in the spline function. For given data points, there are knots and subdomains. For a spline of order , there are polynomial coefficients within each subdomain, for a total of unknowns for the entire spline. The function value at each subdomain endpoint is specified, giving constraints. This gives at least continuity across the function. Specifying continuity of order at each subdomain junction adds total constraints, for a total of constraints. Since the number of constraints must be less than or equal to the number of unknowns, we find that function interpolation using ^{th} order splines allows . The interpolation is performed by solving for the unknown with these constraints. This interpolation process can be done using splines written as a linear combination of Bsplines or using the simple piecewise polynomial representation.
Where many data points are given, as in this study, it is useful to perform a spline fit instead of a spline interpolation. This allows fewer knots in the spline than data points, which can smooth out noise in the data and simplify the evaluation of the function. The curve fit can be done using a least square spline approximation. The Bspline representation of a spline allows us to express the least squares problem as follows:
(11) 
where is the vector of coefficients, is the Bspline order, is the number of interior knots, and is the ^{th} data point (Dahlquist and Björck, 2008).
We use cubic splines to fit in the chemical potential function. We perform a spline fit, as opposed to a spline interpolation, to smooth out the noise from the data and simplify the representation and evaluation of the function. The Octave splinefit function takes as input the data points and allows us to specify the location of the knots in the spline fit. It performs a least squares fit to the given data using Bsplines. However, this function returns the spline’s coefficients for the simple piecewise polynomial form shown in Eq. (10). This set of coefficients can then be evaluated for a given value of in Octave using the ppval function or manually as standard piecewise polynomials.
To accurately capture local features of the data while minimizing the number of knots, we first identify regions where the data are changing rapidly. In the final spline fit, we use a higher density of knots in those regions. To do this, we first fit a cubic spline to the data with knots at relatively small intervals of 0.005 to approximate . We use this initial spline fit of to approximate , and we identify any region where as a region with rapidly changing data. For the final fit, we specify knots at intervals of 0.004 within the regions of rapidly changing data and at intervals of 0.03 outside those regions. This allows us to capture all the physics while reducing the number of knots. The function splinefit is used with this new set of knots to fit the data and update our spline representation of .
3.2 Results
3.2.1 Data with no spinodal regions
At 1800 K, the titanium oxide system has a simple free energy landscape with no spinodals, and a degree three RK polynomial expansion gives a sufficient fit of the chemical potential, as does the cubic spline. The chemical potential data and curve fit are shown in the left plot of Figure 5. The center plot shows the derivative of the chemical potential fit and the numerically differentiated data. This is also the second derivative of the free energy with respect to composition, which shows where the free energy curve is concave or convex. When the second derivative of the free energy is negative, the free energy is concave and represents a twophase region. The numerical derivative was taken by first smoothing the data with a running average to remove some of the noise, then performing a central difference derivative. The plot on the right shows the free energy found by integrating the curve fit and data. It is plotted with respect to the end members at and . We also present a cubic spline fit of the data. Note that the cubic spline, while it also fits the data well, does not give a significantly better representation of this set of data.
3.2.2 Data with three spinodal regions
At 800 K, the HCP form of TiO exhibits three spinodal regions. To accurately capture this physics, the function representing the free energy must have three corresponding concave regions. Fewer data points were calculated with the variance constrained Monte Carlo within these spinodal regions. They were weighted to help capture the spinodals with the polynomial function. The data are fitted with a degree three RK polynomial as before. The data and fit are plotted in Figure 6 along with their corresponding derivatives and free energy. Note that the chemical potential data was smoothed before computing the numerical data to reduce the effects of noise. Also, as in all plots in this work, the free energies are plotted with respect to the end members at and .
The data diverge due to an ordering at , but the equation used here does not capture this divergence. We can represent the divergent behavior of the system by scaling the logarithmic term. This gives us the modified form
(12) 
to model the chemical potential curve over the domain . This corresponds to a free energy density of
(13) 
For consistency, we also rescale the RedlichKister polynomial to the same domain of .
We now fit to the difference between the chemical potential data and . The resulting fit is shown in Figure 7. Note that the divergence of the data at is now represented in the curve fit. It is clear, however, that a degree three RK polynomial cannot represent the spinodal regions. Polynomials of higher degree are used to fit the chemical potential and are plotted in Figure 8. The curve fit represents a spinodal region wherever the derivative of the chemical potential is negative, meaning the free energy is concave (see the center plots of Figure 8). Using the curve fit methods described here, the degree 21 polynomial capture all three spinodals. The degree 15 polynomial represents only two (the first and third), and it nearly produces a spurious spinodal at composition 0.05. The polynomial of degree nine captures one (the first spinodal), and degree five polynomial does not capture any spinodal regions. Clearly, to accurately represent all three spinodals with a RK polynomial, a polynomial of high degree must be used.
A cubic spline was also used to fit . Compare these results, shown in Figure 9, with the previous RK polynomial fits. The spline clearly captures all three twophase regions. This is done with minimal spurious oscillations in the derivative of the chemical potential and no spurious spinodal regions (see the center plot of Figure 9).
It should be noted that the failure of low order RedlichKister polynomials to capture the various spinodals of the firstprinciples free energy curve of TiO at 800K signifies that the initial model based on Eq. (13) and Eq. (12) is too simple to capture the complex entropic contributions when ordered and disordered phases can coexist. The logarithmic composition dependence of Eq. (3) and Eq. (17) represents configurational entropic contributions of a thermodynamically ideal solid solution. Strong deviations from thermodynamic ideality occur at stoichiometric compositions where ordering occurs. This is evident in the oxygen chemical potential of TiO at 800K shown in Figure 3 where steps appear at x=1/6 and 1/2. The nonideal behavior at compositions where ordered phases are stable are more accurately captured with compound free energy models Hillert (2001). However to date, these models have been developed on a case by case basis and can become quite cumbersome for complex ordered phases such as the substoichiometric TiO oxides at x=1/6, 1/3 and 1/2. They rely on sublattice concentration variables and approximate the entropy assuming random mixing within each sublattice. As such models are more accurate in describing entropy due to dilute mixing on the different sublattices of an ordered phase, they require less terms in Redlich Kister expansions that describe the excess free energy. While such an approach is desirable for a rigorous CALPHAD assessment of phase stabilityFischer (1997); Waldner and Eriksson (1999), it is less practical when developing free energy descriptions for phase field simulations due to the presence of additional sublattice composition variables Zhu et al. (2002); Kitashima et al. (2008); Zhang et al. (2015). In that respect, splines, as will become evident in the next section are more versatile in empirically representing experimental or firstprinciples free energy curves needed as input for kinetic modeling.
4 Phase field computations
In this section we explore the extent to which phasefield simulations are affected by different representations of the free energy. To this end we perform phase field simulations of oxygen indiffusion in metallic Ti using the various RedlichKister polynomial parameterizations and the spline fits of the free energy of HCP TiO described in the previous section.
Oxygen indiffusion in HCP Ti is an important process in a variety of applications, the most obvious one being during the oxidation of Ti (Unnam et al., 1986; Pouilleau et al., 1997). The exposure of Ti to an oxygen rich environment quickly results in the formation of oxides such as rutile (Ting et al., 2000; D. S. R. Krishna, 2005; Jamesh et al., 2013) or anatase TiO and rock salt TiO (Kao et al., 2011; Chung et al., 2011, 2012; Okazumi et al., 2011). A layered microstructure typically forms with the most oxygen rich phases appearing at the surface and the oxides having a lower oxygen to titanium ratio forming below the surface (Unnam et al., 1986). The high oxygen solubility within HCP Ti should also result in a sizable region below these distinct oxide phases that is characterized by a large oxygen concentration gradient that drives further oxygen diffusion towards the interior of the metal (Unnam et al., 1986; Dong et al., 1997; Dong and Li, 2000). While Tioxides with high oxygen content form at ambient oxygen partial pressures, they can be suppressed at sufficiently low oxygen partial pressures where they are no longer thermodynamically stable (Dong and Li, 2000). Under those conditions, oxygen will simply dissolve in HCP Ti and form an HCPbased TiO solid solution at the surface along with the ordered HCP TiO phases having =1/6, 1/3 and 1/2.
Similar oxygen diffusion processes in HCP Ti are relevant within the recently proposed Tiair battery, consisting of a Ti anode separated from an oxygen rich cathode by a solidoxide electrolyte phase (Van der Ven et al., 2013a). The externally imposed voltage (EMF) enables control over the oxygen chemical potential of the anode and high voltages and capacities are theoretically achievable when cycling the Ti anode between pure Ti and oxygen rich TiO solid solutions (). By limiting the voltage to a narrow window, it is possible to suppress the formation of the oxygen rich Tioxides.
Here we use phase field simulations to model oxygen diffusion in HCP Ti assuming zero flux boundary conditions on the surface. This approximates the process of oxygen diffusion through the metallic Ti substrate of an oxidizing sample and through the Ti anode during electrochemical cycling of a Tiair battery after the source of oxygen at the surface has been cut off.
The TiO solid solution as well as the ordered TiO, TiO and TiO suboxides all share the same HCP Ti sublattice. The ordered suboxides therefore have a group/subgroup symmetry relationship with the disordered TiO solid solution. This property makes it possible to describe oxygen indiffusion and orderdisorder phase transformations with a combined CahnHilliard and AllenCahn description (Cahn and Hilliard, 1958; Hilliard, 1970; Allen and Cahn, 1979; Chen, 2002; Balluffi et al., 2005). The overall free energy of the solid can be written as an integral over a free energy density that is a function of the local concentration and local values of the order parameters as well as the gradients of the local composition and order parameters according to:
(14) 
where and correspond to the Cartesian , and directions and where is the volume of TiO per Ti atom. The tensor elements, and , represent gradient energy coefficients. In the above free energy expression, we have neglected gradient energy terms that couple a gradient in concentration and a gradient in the order parameters , which may be allowed depending on how the order parameters transform under symmetry operations of the crystal.
The temporal evolution of the solid when the initial total free energy is not minimal can be described with coupled CahnHilliard and AllenCahn evolution equations. The CahnHilliard equation in this example describes longrange oxygen diffusion with the oxygen flux, , driven by gradients in the oxygen chemical potential according to (de Groot and Mazur, 1984)
(15) 
where is the Onsager transport coefficient (Van der Ven et al., 2010b, 2013b) for oxygen diffusion over the interstitial sites of HCP Ti. Combining this flux expression with the continuity equation (i.e. conservation of oxygen atoms)
(16) 
yields a CahnHilliard differential equation that links the time derivative of the oxygen concentration to spatial derivatives of the oxygen chemical potential . The oxygen chemical potential in a solid having a nonuniform concentration profile is the variational derivative of the total free energy
(17) 
These equations are sufficient to describe oxygen redistribution in the solid solution phase above the orderdisorder transition temperatures of TiO, TiO and TiO. However, at 800 K, a concentration gradient will result in local concentrations where the ordered TiO and TiO suboxides are stable. An AllenCahn description is then necessary to complement the above CahnHilliard equations to describe oxygen ordering tendencies. For each order parameter , there is an AllenCahn evolution equation of the form (Chen, 2002)
(18) 
where M is a kinetic coefficient that determines the rate of change in the degree of local ordering in the presence of a free energy driving force. It is determined by local atomic hop events that are required to change the local degree of ordering.
Since ordering processes at fixed concentration require only shortrange diffusion involving a small number atomic hops, they will relax substantially more rapidly than processes that require longrange diffusion as described by the Onsager transport coefficients, . Hence, while both CahnHilliard and AllenCahn processes occur simultaneously, the AllenCahn evolution equations will reach a local minimum much more rapidly than the longrange diffusion processes described by CahnHilliard. This therefore motivates an effective CahnHilliard dynamics, in which AllenCahn ordering processes occur sufficiently rapidly that the local order parameters, , are equal to the values that minimize the free energy at the local concentration . This has as a consequence that the order parameters are functions of the local concentration, i.e. . The above total free energy, Eq. (14), can then be converted into a free energy that only explicitly depends on the concentration profile
(19) 
where the homogeneous free energy is defined according to Eq. (13) and the gradient energy coefficients are functions of , and the equilibrium relationships between order parameters and composition .
The CahnHilliard phase field simulations of oxygen diffusion in HCP Ti implicitly rely on the above assumptions. We also neglect contributions to the free energy from coherency strains due to the dependence of lattice parameters on composition. These effects can be accounted for by explicitly making the free energy expression a function of local strains and variationally minimizing with respect to a displacement field to yield mechanical equilibrium criteria (Chen, 2002; Voorhees and Johnson, 2004; Rudraraju et al., 2014, 2016). While the equilibrium volume of HCP based TiO does depend on , this variation is at most only several percent (Van der Ven et al., 2013a).
The resulting equation for the chemical potential, equal to the variational derivative of the free energy in Eq. (19), is given by
(20) 
where is the homogeneous chemical potential as defined in Eq. (1) and (12). Considering the case where simplifies this equation further:
(21) 
By inserting this equation for the chemical potential into Eq. (15) and (16), we have the following CahnHilliard equation:
(22) 
The parameter and the average homogeneous chemical potential give a length scale for the interface . Substituting the expression for the homogeneous chemical potential in Eq. (12) gives the following:
(23) 
The corresponding boundary conditions are
(24)  
where the first boundary condition results from assuming equilibrium at the boundary, and the second represents an influx at the boundary.
Because the cubic spline representation of the excess chemical potential is central to this study, we insert the spline function into Eq. (23). The final form for the local CahnHilliard equation, using the piecewise polynomial spline representation from equation (10) for , is given by
(25) 
where the piecewise cubic polynomials are given by
(26) 
The polynomial coefficients used in the computations presented here were obtained using the Octave splinefit function.
As described previously, the spline function can also be expressed as a linear combination of Bsplines, following Eq. (9). This alternative form gives the following CahnHilliard equation:
(27) 
In addition to the phase field computations using a cubic spline fit for the excess chemical potential, we also performed simulations using the RedlichKister polynomial representation. A comparison of the resulting composition profiles, total energy, and computation times is presented here.
We solve the weak form of the CahnHilliard equation with Isogeometric Analysis (IGA), which is a meshbased numerical method that uses NURBS (NonUniform Rational BSplines) as basis functions (Hughes et al., 2005; Cottrell et al., 2009). The code was implemented using the PetIGA library (Collier et al., 2013).
4.1 Composition profiles
To demonstrate the effect of the curve fits from the previous section, we performed 2D simulations of the diffusion of oxygen within metallic Ti at 800K. We fit the data for the excess chemical potential with the cubic spline fit shown in Figure 9 and the RK polynomial fits of degree 15 and 21 shown in Figure 8. The domain was with a element mesh. Initial conditions were a uniform composition gradient between at and at . Zero flux boundary conditions were applied. The average value of the homogeneous chemical potential is . We use , which gives a length scale of about 0.02 . This is well resolved by the element length of 0.002 in the direction.
Figure 10 shows the composition profile after 100 s for the three simulations, as well as the derivative of the homogeneous component of the chemical potential with respect to composition, . Recall that the twophase regions correspond to the concave regions of the free energy, which is identified by a negative value for . While the 800K data has three spinodal regions between and , recall from Figure 8 that the degree 15 RK polynomial does not capture the middle twophase region. The effect of this is seen in Figure 10. Note that the cubic spline and the degree 21 polynomial have a welldefined interface near , but this interface is smoothed out when using the degree 15 polynomial. Additionally, because we are using a constant mobility , every oscillation in the excess chemical potential function is reflected in the composition profile. Some of these oscillations are in the data itself; others result from the curve fit and can appear to be interfaces. However, even when it is difficult to identify the spinodal regions from the composition profile alone, the derivative of the homogeneous chemical potential clearly represents what is a true twophase region.
4.2 Total free energy as a measure of accuracy
We compare the accuracy of the free energy representations by evaluating the time evolution of their total free energies, F, given by Eq. (19). Twodimensional simulations with the same conditions described previously ran for 1000 time steps of 10 s each, using a constant mobility of . A thickness of was used in the calculation of total energy. These results are compared to the simulation that used the cubic spline fit, which is considered to most accurately represent the chemical potential of the physical system (see Figure 11). We would expect a good fit of the chemical potential to compute a total energy that matches the energy when using the spline fit.
The zero flux boundary conditions, together with the initial uniform composition gradient, result in a composition profile that gradually approaches a uniform composition of in each simulation. The degree five RK polynomial clearly overestimates the magnitude of the total free energy. To a lesser degree, the degree nine order polynomial underestimates the magnitude of the total free energy. This is consistent with the computed homogeneous free energy curves shown in Figure (8), where the degree 5 polynomial fit for the excess chemical potential results in a free energy curve significantly lower than the reference curve. The degree 9 RK polynomial is slightly above the reference. The other two RK polynomial fits and the spline fit match the reference free energy curve relatively well, which is reflected in the simulation results of Figure 11.
4.3 Assembly time
We also compare the effect on computation time of using the cubic spline in comparison with the degree 21 polynomial. Evaluation of the chemical potential and its derivatives takes place when the residual and Jacobian are evaluated and assembled in the phase field code. The additional terms in the degree 21 polynomial relative to the cubic spline incur a significant increase in the number of floating point operations in the computation. To quantify this effect, 3D simulations of oxygen diffusion in metallic Ti were performed on a cube of length with a element mesh using 160 processors. Zero flux boundary conditions were applied, with initial conditions of a uniform composition gradient between at and at .
Wall times for assembly of the residual and the Jacobian for the first time step are shown in Table 1. There is a significant difference in the assembly times, with the polynomial fit taking more that 15 times as long to evaluate and assemble the residual and about 5 times as long for the Jacobian. This difference in computation time continues long after the first time step. The wall times, averaged over each time step, are plotted in Figure 12 for the first 100 time steps. The wall times for both simulations remain essentially constant over that time, with the residual and Jacobian taking about five and 15 times longer, respectively, when using the degree 21 polynomial instead of the cubic spline.
Note that for a system where a polynomial of lower degree provides a sufficient fit, the difference in wall time would decrease. For example, if we simulated the oxgyen diffusion in titanium at 1800K, as in Figure 5, a cubic polynomial would sufficiently represent the excess chemical potential, and there would be no decrease in computation time from using a cubic spline instead. However, for a complex free energy landscape, such as present in the diffusion of oxygen in Ti at 800K, the decrease in computation time achieved by using a cubic spline as opposed to a polynomial of high degree is significant.
Newton  Residual time (s)  Jacobian time (s)  
iteration  Polynomial  Spline  Polynomial  Spline 
0  21.0  1.37  1230  249 
1  21.0  1.40  1240  249 
2  20.9  1.32  1230  249 
3  21.3  1.34  1230  249 
4  20.9  1.32  1240  250 
5 Conclusion
RedlichKister polynomials have commonly been used to represent free energy functions. The global nature of the these polynomials often prevents them from capturing local phenomena. Spline functions present an effective alternative for representing chemical potential and free energy data. Splines are piecewise polynomials with a specified order of continuity across the entire domain. Because of their piecewise structure, they are able to accurately represent local features of the data. The inherent constraints prescribing that the spline and a certain number of its derivatives be continuous across all subdomain junctions give splines an order of global continuity. These local and global properties make splines an important tool in curve fitting.
For a simple landscape, the RK polynomial sufficiently models the chemical potential with a low polynomial degree. This is demonstrated by the high temperature data for diffusion of oxygen in metallic Ti, where a cubic RK polynomial gave a sufficient representation of the excess chemical potential. This contrasts, however, with more complex cases involving multiple regions of phase separation, where much higher polynomial degrees are needed when using the RK polynomials. At the lower temperature of 800K, for example, the diffusion of oxygen in Ti produces three spinodal regions, and the RK expansion requires a polynomial degree of around 20 to represent all three. These polynomials of high degree add complexity and can introduce spurious oscillations, but polynomials of lower degree potentially miss significant physics, such as the twophase regions corresponding to a negative curvature of the free energy. Splines, however, are able to represent the fitted data and the associated physics with high accuracy while using a low degree.
To demonstrate the effect of the function type used to fit the excess chemical potential, we performed phase field computations of the diffusion of oxygen in titanium using a cubic spline fit and RK polynomial fits of multiple degrees. Because the AllenCahn equations reach a minimum more rapidly than the CahnHilliard equation, our phase field model is based on an effective CahnHilliard equation. Using this model, we compared the performance of a cubic spline fit to the RK polynomial fits in three ways. First, we observed that, for this example, RK polynomials of degree less than 21 fail to represent all three welldefined phase interfaces, while the cubic spline does. We also compared the evolution of the total free energy of the system as it relaxes toward a homogeneous composition, and we found that the degree five and, to a lesser extent, the degree nine RK polynomials misrepresent the magnitude of the total free energy. From these two studies, we concluded that a degree 21 RK polynomial gives a sufficient representation of the excess chemical potential for this particular data. However, we then compared the computation time for the assembly of the residual and Jacobian using the cubic spline and the degree 21 RK polynomial. The computation with the cubic spline evaluated and assembled the residual with a speedup of 15 and the Jacobian with a speedup of 5 relative to the computation using the degree 21 polynomial. The ability of the cubic spline to capture all of the local physics of the system while minimizing spurious oscillations and computation time makes it a valuable tool in representing chemical potential and free energy functions.
Acknowledgements
This work was carried out under an NSF DMREF grant: DMR1436154 “DMREF: Integrated Computational Framework for Designing Dynamically Controlled AlloyOxide Heterostructures”. The computations were performed using resources supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award #DESC0008637 that funds the PRedictive Integrated Structural Materials Science (PRISMS) Center at University of Michigan.
References
References
 Hillert [2007] M. Hillert. Phase Equilibria, Phase Diagrams and Phase Transformations. Cambridge University Press, Cambridge, 2nd edition, 2007.
 Kaufman and Agren [2014] L. Kaufman and J. Agren. CALPHAD, first and second generation â Birth of the materials genome. Scripta Materialia, 70:3–6, 2014.
 Voorhees and Johnson [2004] P. W. Voorhees and W. C. Johnson. The thermodynamics of elastically stressed crystals. Solid State Phys. Adv. Res. Appl., 59:1–201, 2004.
 DeHoff [1993] R. T. DeHoff. Thermodynamics in Materials Science. McGrawHill, 1993.
 Sanchez et al. [1984] J. M. Sanchez, F. Ducastelle, and D. Gratias. Generalized cluster description of multicomponent systems. Physica A, 128:334–350, 1984. doi: 10.1016/03784371(84)900967.
 Fontaine [1994] D. De Fontaine. Cluster approach to orderdisorder transformations in alloys. volume 47 of Solid State Physics, pages 33 – 176. Academic Press, 1994. doi: http://dx.doi.org/10.1016/S00811947(08)606396. URL http://www.sciencedirect.com/science/article/pii/S0081194708606396.
 Laks et al. [1992] D. B. Laks, L. G. Ferreira, S. Froyen, and A. Zunger. Efficient clusterexpansion for substitutional systems. Phys. Rev. B, 46, 1992. doi: 10.1103/PhysRevB.46.12587.
 Ceder [1993] G. Ceder. A derivation of the ising model for the computation of phase diagrams. Computational Materials Science, 1:144–150, 1993.
 Van der Ven et al. [1998] A. Van der Ven, M. K.. Aydinol, G. Ceder, G. Kresse, and J. Hafner. Firstprinciples investigation of phase stability in LiCoO. Phys. Rev. B, 58:2975–2987, 1998.
 van de Walle and Ceder [2002] A. van de Walle and G. Ceder. The effect of lattice vibrations on substitutional alloy thermodynamics. Reviews of Modern Physics, 74:11–45, 2002.
 van de Walle and Asta [2002] A. van de Walle and M. Asta. Selfdriven latticemodel Monte Carlo simulations of alloy thermodynamic properties and phase diagrams. Modelling and Simulation in Materials Science and Engineering, 10:521–538, 2002.
 Arroyave et al. [2005] R. Arroyave, D. Shin, and Z. K. Liu. Ab initio thermodynamic properties of stoichiometric phases in the NiAl system. Acta materialia, 53:1809–1819, 2005.
 Van der Ven et al. [2010a] A. Van der Ven, J. C. Thomas, Q. Xu, and J. Bhattacharya. Linking the electronic structure of solids to their thermodynamic and kinetic properties. Mathematics and Computers in Simulation, 80(7):1393–1410, 2010a.
 Thomas and Ven [2013] John C. Thomas and Anton Van der Ven. Finitetemperature properties of strongly anharmonic and mechanically unstable crystal phases from first principles. Phys. Rev. B, 88:214111, Dec 2013. doi: 10.1103/PhysRevB.88.214111. URL http://link.aps.org/doi/10.1103/PhysRevB.88.214111.
 Puchala and Van der Ven [2013] B. Puchala and A. Van der Ven. Thermodynamics of the ZrO system from firstprinciples calculations. Phys. Rev. B, 88:094108, Sep 2013. doi: 10.1103/PhysRevB.88.094108. URL http://link.aps.org/doi/10.1103/PhysRevB.88.094108.
 Chen et al. [2015] M.H. Chen, B. Puchala, and A. Van Der Ven. Hightemperature stability of ZrO. Calphad: Computer Coupling of Phase Diagrams and Thermochemistry, 51:292–298, 2015. doi: 10.1016/j.calphad.2015.10.010.
 Natarajan et al. [2016] A. R. Natarajan, E. L. S. Solomon, B. P. Puchala, E. Marquis, and A. Van der Ven. On the precipitation sequence in dilute MgNd alloys. Acta Materialia, 108:367–379, 2016.
 Redlich and Kister [1948] Otto Redlich and A. T. Kister. Algebraic representation of thermodynamic properties and the classification of solutions. Industrial & Engineering Chemistry, 40(2):345–348, 1948. doi: 10.1021/ie50458a036. URL http://dx.doi.org/10.1021/ie50458a036.
 Dahlquist and Björck [2008] Germund Dahlquist and Åke Björck. Numerical methods in scientific computing, volume 1. Philadelphia: Society for Industrial and Applied Mathematics, 2008.
 Tomiska [1984] Josef Tomiska. Mathematical conversions of the thermodynamic excess functions represented by the RedlichKister expansion, and by the Chebyshev polynomial series to power series representations and viceversa. Calphad, 8(4):283 – 294, 1984. ISSN 03645916. doi: http://dx.doi.org/10.1016/03645916(84)900324. URL http://www.sciencedirect.com/science/article/pii/0364591684900324.
 Ouerfelli et al. [2010] N. Ouerfelli, O. Iulian, and M. Bouaziz. Competition between RedlichKister and improved HerrÃ¡ez equations of correlation viscosities in 1,4dioxane+water binary mixtures at different temperatures. Physics and Chemistry of Liquids, 48(4):488–513, 2010. doi: 10.1080/00319100903131559. URL http://dx.doi.org/10.1080/00319100903131559.
 Cahn and Hilliard [1958] John W. Cahn and John E. Hilliard. Free energy of a nonuniform system. I. Interfacial free energy. The Journal of Chemical Physics, 28(2):258–267, 1958. doi: http://dx.doi.org/10.1063/1.1744102. URL http://scitation.aip.org/content/aip/journal/jcp/28/2/10.1063/1.1744102.
 Cancarevic et al. [2007] M. Cancarevic, M. Zinkevich, and F. Aldinger. Thermodynamic description of the TiâO system using the associate model for the liquid phase. Calphad, 31(3):330 – 342, 2007. ISSN 03645916. doi: http://dx.doi.org/10.1016/j.calphad.2007.01.009. URL http://www.sciencedirect.com/science/article/pii/S0364591607000247.
 Okamoto [2011] H. Okamoto. OTi (oxygentitanium). Journal of Phase Equilibria and Diffusion, 2011.
 Burton and van de Walle [2012] Benjamin Paul Burton and Axel van de Walle. First principles phase diagram calculations for the octahedralinterstitial system TiO, . Calphad, 39:97 – 103, 2012. ISSN 03645916. doi: http://dx.doi.org/10.1016/j.calphad.2012.09.004. URL http://www.sciencedirect.com/science/article/pii/S0364591612000788.
 Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio totalenergy calculations using a planewave basis set. Phys. Rev. B, 54:11169–11186, Oct 1996. doi: 10.1103/PhysRevB.54.11169. URL http://link.aps.org/doi/10.1103/PhysRevB.54.11169.
 Kresse and Joubert [1999] G. Kresse and D. Joubert. From ultrasoft pseudopotentials to the projector augmentedwave method. Phys. Rev. B, 59:1758–1775, Jan 1999. doi: 10.1103/PhysRevB.59.1758. URL http://link.aps.org/doi/10.1103/PhysRevB.59.1758.
 Dalton et al. [2012] A. S. Dalton, A. A. Belak, and A. Van der Ven. The thermodynamics of lithium in TiO(B) from first principles. Chemistry of Materials, 24:1568–1574, 2012.
 Sadigh and Erhart [2012] Babak Sadigh and Paul Erhart. Calculation of excess free energies of precipitates via direct thermodynamic integration across phase boundaries. Phys. Rev. B, 86:134204, Oct 2012. doi: 10.1103/PhysRevB.86.134204. URL http://link.aps.org/doi/10.1103/PhysRevB.86.134204.
 de Boor [1972] Carl de Boor. On calculating with Bsplines. Journal of Approximation Theory, 6(1):50–62, July 1972.
 Hillert [2001] M. Hillert. The compound energy formalism. Journal of Alloys and Compounds, 320:161–176, 2001.
 Fischer [1997] E. Fischer. Thermodynamic calculation of the OTi system. Journal of Phase Equilibria, 18:338, 1997.
 Waldner and Eriksson [1999] P. Waldner and G. Eriksson. Thermodynamic modeling of the system titaniumoxygen. Calphad, 23(189218), 1999.
 Zhu et al. [2002] J. Z. Zhu, Z. K. Liu, V. Vaithyanathan, and L. Q. Chen. Linking phasefield model to CALPHAD: application to precipitate shape evolution in Nibase alloys. Scripta Materialia, 46:401–406, 2002.
 Kitashima et al. [2008] T. Kitashima, J. Wang, and H. Harada. Phasefield simulation with the CALPHAD method for microstructure evolution of multicomponent Nibase superalloys. Intermetallics, 16:239–245, 2008.
 Zhang et al. [2015] L. Zhang, M. Stratmann, Y. Du, B. Sundman, and I. Steinbach. Incorporating the CALPHAD sublattice approach of ordering into the phasefield model with finite interface dissipation. Acta Materialia, 88:156–169, 2015.
 Unnam et al. [1986] J. Unnam, R. N. Shenoy, and R. K. Clark. Oxidation of commercial purity titanium. Oxidation of Metals, 26:231, 1986.
 Pouilleau et al. [1997] J. Pouilleau, D. Devilliers, F. Garrido, S. DurandVidal, and E. Mahe. Structure and composition of passive titanium oxide films. Materials Science and Engineering, B47:235, 1997.
 Ting et al. [2000] C. C. Ting, S. Y. Chen, and D. M. Liu. Structural evolution and optical properties of TiO thin films prepared by thermal oxidation of sputtered Ti films. ournal of Applied Physics, 88:4628, 2000.
 D. S. R. Krishna [2005] Y. Sun D. S. R. Krishna. Effect of thermal oxidation conditions on tribological behavior of titanium films on 316L stainless steel. Surface & Coating Technology, 198:447–453, 2005.
 Jamesh et al. [2013] M. Jamesh, T. S. N. S Narayanan, and P. K. Chu. Thermal oxidation of titanium: Evaluation of corrosion resistance as a function of cooling rate. Materials Chemistry and Physics, 138:565, 2013.
 Kao et al. [2011] C. H. Kao, S. W. Yeh, H. L. Huang, D. Gan, and P. Shen. Study of the TiO to anatase transformation by thermal oxidation of Ti film in air. Journal of Physical Chemistry C, 115:5648–5656, 2011.
 Chung et al. [2011] Y. L. Chung, D. S. Gan, K. L. Ou, and S. Y. Chiou. Formation of anatase and TiO from Ti thin film after anodic treatment and thermal annealing. Journal of the Electrochemical Society, 158:C319, 2011.
 Chung et al. [2012] Y. L. Chung, D. S. Gan, and K. L. Ou. The Ti to TiO and TiO to anatase transformations induced by anodizing and annealing treatments of Ti thin film. Journal of The Electrochemical Society, 159:C133, 2012.
 Okazumi et al. [2011] T. Okazumi, K. Ueda, K. Tajima, N. Umetsu, and T. Narushima. Anatase formation on titanium by twostep thermal oxidation. Journal of Materials Science, 46:2998, 2011.
 Dong et al. [1997] H. Dong, A. Bloyce, P. H. Morton, and T. Bell. Surface engineering to improve tribological performance of Ti6Al4V. Surface Engineering, 13(5):402–406, 1997.
 Dong and Li [2000] H. Dong and X. Y. Li. Oxygen boost diffusion for the deepcase hardening of titanium alloys. Materials Science and Engineering A, 280:303–310, 2000.
 Van der Ven et al. [2013a] A. Van der Ven, B. Puchala, and T. Nagase. Ti and Zr based metalair batteries. Journal of Power Sources, 242:400–404, 2013a.
 Hilliard [1970] J. E. Hilliard. Spinodal decomposition. In H. I. Aaronson, editor, Phase Transformations, pages 497–560. American Society for Metals, Metals Park, Ohio, 1970.
 Allen and Cahn [1979] S. M. Allen and J. W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase boundary coarsening. Acta Metallurgica, 27:1085â1091, 1979.
 Chen [2002] L. Q. Chen. Phasefield models for microstructure evolution. Annual Review of Materials Research, 32, 2002. doi: 10.1146/annurev.matsci.32.112001.132041.
 Balluffi et al. [2005] R. W. Balluffi, A. M. Allen, and W. C. Carter. Kinetics of Materials. John Wiley & Sons, 2005.
 de Groot and Mazur [1984] S. R. de Groot and P. Mazur. NonEquilibrium Thermodynamics. Dover, 1984.
 Van der Ven et al. [2010b] A. Van der Ven, H. C. Yu, G. Ceder, and K. Thornton. Vacancy mediated substitutional diffusion in binary crystalline solids. Progress in Materials Science, 55(2):61–105, 2010b.
 Van der Ven et al. [2013b] A. Van der Ven, J. Bhattacharya, and A. Belak. Understanding Li diffusion in Liintercalation compounds. Accounts of Chemical Research, 46:1216–1225, 2013b.
 Rudraraju et al. [2014] S. Rudraraju, A. Van der Ven, and K. Garikipati. Threedimensional isogeometric solutions to general boundary value problems of Toupinâs gradient elasticity theory at finite strains. Comput. Methods Appl. Mech. Eng., 278:705â728, 2014.
 Rudraraju et al. [2016] S. Rudraraju, A. Van der Ven, and K. Garikipati. Mechanochemical spinodal decomposition: a phenomenological theory of phase transformations in multicomponent, crystalline solids. npj Computational Materials, 2, 2016. doi: 10.1038/npjcompumats.2016.12.
 Hughes et al. [2005] T. J. R. Hughes, J. A. Cottrell, and Yuri Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194:4135–4195, October 2005.
 Cottrell et al. [2009] J. Austin Cottrell, Thomas J. R. Hughes, and Yuri Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA. John Wiley & Sons, Ltd, Chichester, 2009.
 Collier et al. [2013] Nathan O. Collier, Lisandro Dalcín, and Victor M. Calo. PetIGA: Highperformance isogeometric analysis. CoRR, abs/1305.4452, 2013. URL http://arxiv.org/abs/1305.4452.