Combining thermodynamics with tensor completion techniques to enable multicomponent microstructure prediction

Abstract

Multicomponent alloys show intricate microstructure evolution, providing materials engineers with a nearly inexhaustible variety of solutions to enhance material properties. Multicomponent microstructure evolution simulations are indispensable to exploit these opportunities. These simulations, however, require the handling of high-dimensional and prohibitively large data sets of thermodynamic quantities, of which the size grows exponentially with the number of elements in the alloy, making it virtually impossible to handle the effects of four or more elements. In this paper, we introduce the use of tensor completion for high-dimensional data sets in materials science as a general and elegant solution to this problem. We show that we can obtain an accurate representation of the composition dependence of high-dimensional thermodynamic quantities, and that the decomposed tensor representation can be evaluated very efficiently in microstructure simulations. This realization enables true multicomponent thermodynamic and microstructure modeling for alloy design.

Introduction

A number of recent discoveries has largely increased the interest in multicomponent alloy design. High entropy or multi-principle element alloys, for example, can give access to a wide range of properties and combinations of properties that cannot be obtained in alloys based on a single major element.¹ Moreover, the large range over which the diffusion properties of the different elements can vary gives rise to complex precipitation, interdiffusion, and coarsening paths,^2,3,4 resulting in new ways to stabilize metastable phases or structures for extended times.⁵ Adding too many alloying elements may, however, also unnecessarily complicate material recovery and recycling.⁶ Hence, in order to exploit in an optimal and responsible way the opportunities multicomponent alloy design can bring, a profound understanding of the effects of compositional variations and heat treatment on microstructure evolution in these alloys is required. Microstructure simulations are indispensable to obtain these insights given the immense parameter space, including not only the current temperature and the amount of each element, but also many parameters describing the history of the material, such as cooling rate and the duration of a heat treatment or mechanical loading.

Phase-field models (PFM) are generally considered as most suitable to study multicomponent microstructure evolution as diffusion and phase morphology often play a crucial role.^7,8,9 They describe complex morphological changes such as precipitate shape evolution during growth and coarsening and dendrite growth. The effects of various competing thermodynamic driving forces and transport processes can be considered. For multicomponent alloys, the chemical bulk driving force is an important contribution. CALPHAD thermodynamic models are usually adopted to describe the composition and temperature-dependent Gibbs free energy expressions of the different phases. For many technologically relevant alloy systems, CALPHAD models have been assessed to reproduce experimental, theoretical and ab initio data on thermodynamic quantities and phase diagrams.^10,11,12

Different approaches have been proposed to include the CALPHAD thermodynamic descriptions in PFM.^{13,14,15,16,17,18,19,20,21,22,23,24,25,26} A classification and discussion of the different approaches is given further in the paper. However, to date, only a few results for quaternary or higher-order systems are reported in the literature. The reason is that the proposed methods were either for a specific case and cannot be generalised towards systems containing other types of phases, or become prohibitively complex or computation and data-intensive when the number of elements in the alloy is increased. There is thus no robust simulation approach available facilitating in silico microstructure design for multicomponent alloys.

Here, we introduce the use of tensor completion^27,28,29,30 in materials science to enable the representation and use of high-dimensional data sets efficiently. More specifically, we used a canonical polyadic decomposition (CPD) with the factor vectors constrained to polynomial expressions³¹ to include high-dimensional thermodynamic data sets obtained from a CALPHAD model in phase-field simulations. Using CPD, the multicomponent composition dependence can be described using a limited number of coefficients. We show that the accuracy obtained with the proposed approach is comparable or better than that obtained with previous approaches, whereas the applicability is general and the computational cost to evaluate data is low. Moreover, the addition of more elements does not lead to a substantial increase of the complexity, and the number of coefficients and the computational cost to evaluate a CPD depend only linearly on the number of components in the system, in contrast to an exponential increase of the amount of thermodynamic data represented. The observed advantages will thus become increasingly prominent when adding more elements. The capabilities of the approach are illustrated for phase-field simulations of spinodal decomposition and subsequent coarsening, an important phenomenon in multicomponent alloys,³² for Ag-Cu-Ni-Sn liquid alloys.

Results and discussion

Tensor completion of thermodynamic data

In material science, we are familiar with the tensorial character of direction-dependent material properties, like stresses and strains. However, in mathematical engineering, the term tensor has been used much more widely. In general, tensors are higher-order generalisations of vectors (first-order tensors) and matrices (second-order tensors).^33,34 They can be seen as multiway arrays of any order holding numerical values. The term thermodynamic data tensor (TDT) is introduced here to denote a tensor that is obtained by calculating and storing the data set of any thermodynamic quantity of a phase, evaluated from a CALPHAD model. The TDTs typically needed in PFM implementations are the molar Gibbs free energy, ${\mathcal{G}}$, the diffusion potential of each component $n$, ${{\mathcal{P}}}^{(n)}$, and the derivative of the diffusion potential of the component $n$ with respect to the molar fraction of components $m$, ${{\mathcal{D}}}^{(n,m)}$.

Typically, the molar fractions ${x}_{n}$ of the components are used to indicate an alloy composition. For an alloy with $C$ components, the composition dependence of any thermodynamic quantity can be expressed as a function of $C-1$ molar fractions ${x}_{n}$, with $n=1,\ldots ,C-1$ referring to $C-1$ of the $C$ components. The order $N$ of the TDT, this is the number of dimensions of the tensor, then equals $C-1$. A graphical illustration of a third-order TDT, representing composition dependence of a thermodynamic quantity is shown in Fig. 1b. As by definition ${x}_{C}+{\sum }_{n=1}^{N}{x}_{n}=1$ and $0\,<\,{x}_{n}\,<\,(1-{x}_{C})$, tensor entries (these are the values in the tensor) for a combination of molar fraction values disobeying this constraint are not physical and cannot be calculated. A TDT is therefore necessarily incomplete. Besides composition dependence, other dependencies of the thermodynamic function, such as temperature, can be considered, increasing the order of the tensor with one for each extra dependency.

Fig. 1: Characteristics of a thermodynamic data tensor.

a The number of entries in a TDT is plotted as a function of the step size ${\delta }_{x}$ (used for data collection) and the number of components in the system. b Visualisation of two third-order incomplete TDTs sampled with different step size. b The dashed box represents what would be the shape of a full third-order tensor; however, owing to the constraint ${\sum }_{i=1}^{C}{x}_{i}=1$ on the molar fractions, the TDT is incomplete and only the entries within the indicated tetrahedron, i.e., the entries satisfying ${\sum }_{i=1}^{N}{x}_{i}\le 1$ can be evaluated. The step size ${\delta }_{x}$ determines the density and hence the number of entries in the TDT. On the other hand, the number of TDT entries depends exponentially on the number of components, or, more generally, on the order $N$ of the tensor, making storage, and computation quickly intractable, as can be seen in a.

If a thermodynamic model is available, TDTs can be precalculated and used in phase-field simulations to get the required thermodynamic quantities. However, the great challenge of dealing with TDTs is that their number of entries increases exponentially with the order of the tensor, as illustrated in Fig. 1a. Therefore, when more components are considered, it quickly becomes impossible to generate, store, or handle the huge data sets. The computational difficulties caused by this exponential dependence are known as the curse of dimensionality.²⁸

Moreover, as illustrated in Fig. 1a, b, the step size, ${\delta }_{{x}_{n}},$ with which the molar fraction ${x}_{n}$ is sampled, obviously also affects the number of entries in the TDT, as the size of the tensor along each dimension ${I}_{n}$ is inversely proportional to the step size. The step size can be chosen differently along each dimension, but in this study we take a same step size for all molar fractions, i.e., ${\delta }_{x}\mathop{=}\limits^{{\rm{def}}}{\delta }_{{x}_{1}}=\ldots ={\delta }_{{x}_{N}}$. If TDTs are used to provide thermodynamic quantities in phase-field simulations, a small step size ${\delta }_{x}$ is desirable to avoid sudden molar fraction jumps and the need for multi-directional interpolation between entries.

In practice, the use of TDTs in microstructure evolution simulations is thus limited to lower-order systems. The hypothesis motivating our work is that, if a higher-order TDT is too large to be constructed or used, the information present in the TDT may still be accessible if a CPD of the TDT is obtained. See Methods section Tensor Decomposition for a general introduction to this methodology. When applied to a quaternary thermodynamic system ($N=3$), such as the liquid Ag-Cu-Ni-Sn alloys under study, our technique decomposes the TDT into a sum of $R$ rank-1 terms, of which each term is the outer product of $N$ factor vectors, and with $R$ called the rank of the decomposed tensor model. This is graphically illustrated in Fig. 2a, b. An alternative representation of the same CPD is obtained by collecting all factor vectors of the CPD related to each dimension $n$ into a factor matrix ${{\bf{A}}}^{(n)}$. In this case, the dimensions $n$ refer to the molar fractions ${x}_{n}$, with $n=1,2,3$, and factor matrices ${{\bf{A}}}^{(1)}$, ${{\bf{A}}}^{(2)}$, and ${{\bf{A}}}^{(3)}$ are obtained, as shown in Fig. 2a–c.

Fig. 2: Canonical polyadic decomposition of a thermodynamic data tensor (TDT).

a–c Visualisation of a canonical polyadic decomposition (CPD) when applied to a third-order TDT with dimension ${I}_{1}\times {I}_{2}\times {I}_{3}$. The third-order TDT a is written as a sum of $R$ rank-1 terms b, each of which is represented by the outer product of three nonzero factor vectors ${{\bf{a}}}_{r}^{(1)}$, ${{\bf{a}}}_{r}^{(2)}$, and ${{\bf{a}}}_{r}^{(3)}$. The collection of $R$ factor vectors into factor matrices c with dimensions ${I}_{1}\times R,{I}_{2}\times R,$ and ${I}_{3}\times R$, provides a convenient representation of the CPD. Each factor matrix contains the coefficients necessary to describe the contribution of the molar fraction of one of the components to the TDT. The entire tensor in a corresponds to the dashed box in Fig. 1b; the CPD in b will be determined from only the tetrahedron part in Fig. 1b. d The number of entries in the TDT and the number of coefficients in the CPDs with rank values R = 3, 6, 10 are plotted as a function of the number of components, for step size ${\delta }_{x}=0.0001$. The exponential dependence of the number of entries in the TDT on the number of components is broken when a CPD representation is used instead. Indeed, the CPD needs only $({I}_{1}+{I}_{2}+{I}_{3})R$ coefficients, which depends linearly on the number of components. e The same information as in d is plotted, but excluding the TDT size for a better visualisation of the number of coefficients in the TTM as a function of the number of system components.

A major benefit of this approach is demonstrated in Fig. 2d, where the number of entries in a TDT is shown to grow exponentially with the number of components in the system, whereas the number of parameters in a CPD, representing the same TDT only grows linearly. For alloys with four or more components, the reduction in coefficients when using a CPD representation of a full tensor is gigantic, enabling multicomponent microstructure evolution simulations.

While determining the rank value, $R$ is needed for an exact decomposition of a TDT is difficult, previous applications in other fields have shown that tensors representing physical data can often be well approximated using a low-rank CPD.^35,36,37,38 In this work, we applied tensor decomposition to high-dimensional thermodynamic data sets and we show that a low-rank approximation can be used to obtain an accurate decomposed tensor model describing the dependence of Gibbs energies and derived thermodynamic functions on the molar fractions, representing the alloy composition.

Moreover, when a priori knowledge on solution thermodynamics is taken into account, a CPD can be made even more compact and a continuous description of the composition dependence can be obtained. In this study, for example, we exploited the polynomial dependence of thermodynamic quantities on the molar fractions when the ideal mixing term is subtracted. This approach is reasonable as the contribution from ideal mixing depends only on known information (namely, the molar fractions of all components and the temperature) and can be calculated separately and added to the values evaluated from the CPD in the phase-field simulation (see Methods sections CALPHAD thermodynamic model and Thermodynamic tensor model). A polynomial constraint could therefore be imposed on the factor vectors,^31,39,40 and the CPDs representing the molar Gibbs free energy, ${}^{\ast}\!{\mathcal{G}}$ (the superscript ${}^{* }$ denotes that the ideal mixing contribution is subtracted from the thermodynamic quantity), the diffusion potential of component $n$, ${}^{\ast}\!{\mathcal{P}}^{(n)}$, (i.e., derivative of the molar Gibbs free energy) and the derivative of the diffusion potential of component $n$ with respect to the molar fractions of the component $m$, ${}^{\ast}\!{\mathcal{D}}^{(n,m)}$, (i.e., the second derivatives of the molar Gibbs free energies) were coupled. A single set of factor matrices modeling the composition dependence of all the required thermodynamic quantities at once was then obtained. Such a set containing all factor matrices necessary to completely describe the thermodynamic quantities of the system in a phase-field simulation is further referred as a thermodynamic tensor model of rank $R$. For a quaternary system, it is represented as

$${\rm{TTM}}\,R=\left\{{{\bf{A}}}^{(1)},{{\bf{A}}}^{(2)},{{\bf{A}}}^{(3)},{\dot{{\bf{A}}}}^{(1)},{\dot{{\bf{A}}}}^{(2)},{\dot{{\bf{A}}}}^{(3)},{\ddot{{\bf{A}}}}^{(1)},{\ddot{{\bf{A}}}}^{(2)},{\ddot{{\bf{A}}}}^{(3)}\right\}.$$

(1)

The factor matrices ${{\bf{A}}}^{(n)}$ give the contribution related to ${x}_{n}$ to the Gibbs free energy TDT. The ${\dot{{\bf{A}}}}^{(n)}$ model, the contribution related to taking the first partial derivative of the Gibbs energy with respect to ${x}_{n}$, is needed to represent the TDT of the diffusion potential of component $n$. The ${\ddot{{\bf{A}}}}^{(n)}$ model contributions related to taking the second partial derivative of the Gibbs energy with respect to ${x}_{n}$. A cross partial derivative with respect to ${x}_{n}$ and ${x}_{m}$ is obtained using ${\dot{{\bf{A}}}}^{(n)}$ and ${\dot{{\bf{A}}}}^{(m)}$.

From these factor matrices, the thermodynamic tensor models representing the required TDTs can be obtained, for instance, for a quaternary system

$$\begin{array}{ll}{}^{* }{\mathcal{G}}\approx [\kern-2pt[ {{\bf{A}}}^{(1)},{{\bf{A}}}^{(2)},{{\bf{A}}}^{(3)}]\kern-2pt] ,\\ {}^{* }{{\mathcal{P}}}^{(1)}\approx [\kern-2pt[ {\dot{{\bf{A}}}}^{(1)},{{\bf{A}}}^{(2)},{{\bf{A}}}^{(3)}]\kern-2pt] ,\\ {}^{* }{{\mathcal{D}}}^{(1,2)}\approx [\kern-2pt[ {\dot{{\bf{A}}}}^{(1)},{\dot{{\bf{A}}}}^{(2)},{{\bf{A}}}^{(3)}]\kern-2pt] ,\\ {}^{* }{{\mathcal{D}}}^{(3,3)}\approx [\kern-2pt[ {{\bf{A}}}^{(1)},{{\bf{A}}}^{(2)},{\ddot{{\bf{A}}}}^{(3)}]\kern-2pt] ,\end{array}$$

(2)

giving an extremely compact representation of an immense amount of data. Taking, for example, $N=3$, ${\delta }_{x}=0.0001$, and $R=7$, a TDT of about ${(1/{\delta }_{x})}^{N}/N!\,\approx \,1.6\,\times \,1{0}^{11}$ entries (the factor $1/N!$ takes into account that only the entries within the tetrahedron in the TDT in Fig. 1a, b are calculated) is represented using factor matrices that contain only $R\,\times \,N\,\times \,1/{\delta }_{x}\,\approx \,2.1\,\times \,1{0}^{5}$ coefficients. The CPD model is thus about a million times more compact than the TDT. Furthermore, the coefficients of the degree $d$ polynomial constraint alone already suffice for the representation of the TTM. The Gibbs free energy and derived TDTs can thus be represented with $R\times N\times (d+1)=105$ coefficients only. Such impressive compression ratios have recently revolutionised tensor-based scientific computing;^35,36,37,38 in materials science applications, the approach is new.

Moreover, the coefficients in the factor vectors can be obtained using only a limited set of entries from the original TDT, called the training set (see Methods section Thermodynamic tensor model training for more details). There is thus no need to compute the full TDT at any point in the derivation.

Another major advantage of the approach is that values approximating individual entries of the TDT can be calculated very efficiently and easily from the TTM. For example, the diffusion potential of the first component ${\tilde{\mu }}_{1}$ at a composition given by the molar fraction values ${x}_{1}$, ${x}_{2}$, and ${x}_{3}$, for which the entry in the tensor is indexed by ${i}_{1}$, ${i}_{2}$, and ${i}_{3}$, is calculated as

$${\tilde{\mu }}_{1}({x}_{1},{x}_{2},{x}_{3})\approx \sum _{r=1}^{R}{\dot{a}}_{{i}_{1},r}^{(1)}\ {a}_{{i}_{2},r}^{(2)}\ {a}_{{i}_{3},r}^{(3)}+\bar{R}{T}({\log}({x}_{1})-{\log}(1-{x}_{1}-{x}_{2}-{x}_{3})).$$

(3)

Therefore, again we avoid the need to construct explicitly the full TDT ${}^{\ast}\!{\mathcal{P}}^{(1)}$ before data can be extracted. This provides us with an efficient framework to supply PFM with CALPHAD thermodynamic quantities.

First TTMs with rank values $R=3,\ldots ,10$ are optimised using Tensorlab,⁴⁰ modeling the composition dependence of the molar Gibbs energy, the diffusion potentials of the elements Ag, Cu, and Ni, and the partial derivatives of these diffusion potentials with respect to the molar fractions Ag, Cu, and Ni. The data used for training were obtained using Thermo-Calc 2017b with the COST 531 database.⁴¹ Further details are given in Methods section Thermodynamic tensor model training. This particular system was chosen for validation purposes, as we have access to all coefficients of the CALPHAD model for this quaternary system. The rank $R$ of the TTM is an important parameter in this study. A rank as low as possible is preferred, as the size of the factor matrices scales linearly with $R$ (see Fig. 2b, c), and consequently the number of coefficients that has to be optimised and the number of terms that has to be evaluated when extracting data from the tensor model (equation (3)) in the phase-field simulations, increase linearly with $R$ as well.

The dependence of the accuracy of the TTMs on the rank value $R$ is evaluated by comparing thermodynamic quantities approximated with TTMs $R=3,\ldots,10$ with the corresponding entries in the validation TDTs (see Methods section Validation). The empirical cumulative density function (ECDF) of the relative errors (${\rm{RE}}$) on the thermodynamic quantities calculated from the TTMs is plotted for the different rank values in Fig. 3a. The interquartile range (IQR) and $Q=0.95$ and $Q=0.98$ quantiles of the relative error are given for each rank value $R$ in the table in Fig. 3b. The plot and table show that up to $R=7$, the accuracy improves for increasing $R$, whereas for higher values of the rank further improvement is limited. Using the TTM with rank value $R=7$, $98 \%$ of the data points are represented with more than four digits of accuracy, i.e., ${\rm{RE}}\,<\,1{0}^{-4}$. For a small fraction ($<0.1 \%$) of points, the relative error remains large. However, further inspection of these larger RE showed that they are only found when the considered thermodynamic quantity itself has a value close to zero. It is verified, for instance, that TTMs with $R=7$ or higher all have a maximum absolute error smaller than 32 J/mol, which is very small compared with the range over which the Gibbs free energy data vary, namely between −1.02 × 10⁵ J/mol and 3.82 × 10⁶ J/mol.

Fig. 3: Thermodynamics tensor model accuracy.

a Empirical cumulative distribution function (ECDF) of the relative error of TTMs for rank values $R=3,\ldots ,10$. b Table with the interquartile range (IQR) and quantiles $Q=0.95$ and $Q=0.98$ of the relative error of TTMs for rank values $R=3,\ldots ,10$. In both cases a, b the relative error is calculated for the Gibbs free energy, diffusion potentials, and derivatives taking data collected using TC-toolbox with COST 531 database as the reference. There is little improvement of the accuracy of the TTM for rank values $R\ge 7$.

Phase-field simulations

Next, the TTMs for different rank values are used to evaluate the diffusion potentials as a function of the local composition in a multicomponent phase-field model simulating spinodal decomposition and further coarsening in Ag-Cu-Ni-Sn liquids, using the model described in Methods section Phase-field microstructure evolution model. The accuracy and advantages of the approach using TTMs are compared with those of the existing approaches for using CALPHAD thermodynamic models in phase-field simulations. Fig. 4 gives an overview of the most common coupling approaches.^{13,14,15,16,18,19,20,22,23,24,25,26} The methods using direct implementation of the CALPHAD Gibbs free energy expressions (b) or an interface with a thermodynamic software (d), give an exact and continuous representation of the composition dependence of the CALPHAD Gibbs free energies. They are considered as a reference against which the accuracy of the simulations using TTMs is evaluated. In comparison with the method using lookup tables (e–f), the curse of dimensionality is broken and an arbitrarily fine resolution, and even continuous representation, of the composition dependence can be obtained.

Fig. 4: Schemes coupling CALPHAD and phase-field method.

Using the CALPHAD method, an expression for the bulk free energy of the existing phases a is coupled with the phase-field model i via different approaches. b The expression of the CALPHAD Gibbs free energy model and its derivatives with respect to the composition variables are inserted directly in the phase-field model. It is simple if all phases show the behaviour of a substitutional solution. However, most solid phases are described using a sublattice model for which it is complex to relate the phase-field and the CALPHAD variables. c A paraboloid expression fitted to data calculated with the CALPHAD method is used to approximate the composition dependence of the Gibbs free energy. However, this approximation can only describe the Gibbs free energy accurately over limited composition ranges. Moreover, for higher-order systems, it becomes hard to prevent molar fractions from taking nonphysical values below 0 or above 100%. d External software is used to evaluate thermodynamic quantities as required by the phase-field simulation. The time spent in the communication between software is the main disadvantage of this approach, making it less efficient than a. e, f The sampling of data as a function of composition using a thermodynamic software into lookup tables that are consulted in the phase-field simulation, which is heavily affected by the curse of dimensionality. g Proposed methodology in this work: a canonical polyadic decomposition is applied to describe efficiently the thermodynamic quantities, resulting in factor matrices h from which the thermodynamic quantities can be evaluated in PFM. The alloy compositions considered for the phase-field simulations are given in i.

First, 1D simulations were conducted using TTMs and validated against the approach using an interface with thermodynamic software. The conclusions are discussed in the Supplementary Information. Next, 2D simulations were conducted for the six alloy compositions given in Fig. 4i, using TTMs with rank values $R=3,\ldots ,10$, and validated against the approach where the CALPHAD substitutional model is directly implemented in the PFM. For each of the six alloy compositions an initial condition is created by adding small random noise at each position in the system to the values of the molar fractions as given in Fig. 4i. The initial condition and all phase-field simulation parameters are the same for all simulations. The only difference is the approach used to include the CALPHAD diffusion potentials in the PFM.

A selection of the 2D-simulated microstructures is shown in Fig. 5a, c, e. Results obtained using the direct implementation of the CALPHAD model expressions (CE) are shown in the first column and results obtained using TTMs $R=4$ and $R=6$ in the second and third column. Even for $R=4$ and $R=6$, the microstructures obtained using the TTMs look very similar to those obtained using the direct implementation of the CALPHAD model. More results can be found in the Supplementary Information.

Fig. 5: 2D microstructures and probability distribution function of the molar fraction profiles.

a, c, e show a selection of microstructures resulting from the 2D simulations of three different alloys for visual comparison. The first column shows the microstructures obtained with direct use of the CALPHAD full expressions in the phase-field model. The second and third column show the results obtained using the coupling through decomposed TTMs for rank values $R=4$ and $R=6$. Each row represents one alloy and the molar fraction of one component as described by the label. a, c, e show the PDF of the molar fraction values in the different grid cells of the numerical domain. The two means ${m}_{1}$ and ${m}_{2}$ of the bimodal distribution are indicated. The amplitude $A$ of the molar fraction variation in the microstructure is defined as the difference between ${m}_{1}$ and ${m}_{2}$.

For a more quantitative comparison, the corresponding probability density functions (PDF) of the molar fraction fields are visualised in Fig. 5b, d, f. The distributions are bimodal with two means, ${m}_{1}$ and ${m}_{2}$, corresponding to the different compositions of the co-existing liquid phases. According to the Langer-Bar-On-Miller (LBM) method,^42,43,44 the microstructure characteristics can be quantified based on the amplitude $A$, defined as $A=| {m}_{2}-{m}_{1}|$ (see Methods subsection Validation). For all six alloys, the bimodal PDFs obtained from the simulations using the TTM with rank value $R=6$, coincide with those obtained evaluating the full CALPHAD model expressions. When using the TTM $R=4$, a deviation of the PDF from that obtained for CE is observed in some cases, see Fig. 5d. These deviations are also reflected in the values of the amplitude $A$. Therefore, the relative error on the values obtained for the amplitude $A$ from those measured from the profiles obtained using the direct implementation of the CALPHAD model (CE) are considered as an appropriate measure to quantify the accuracy of the phase-field simulations using TTMs with a different rank value $R$.

In Fig. 6a, the CDFs of the RE on the amplitudes of the bimodal distributions of the molar fraction values are plotted for varying rank. Up to rank value $R=6$, the accuracy improves when the rank value of the TTM is increased. For $R\,\ge\, 6$, no further improvement of the accuracy is obtained. The analysis of the IQR given in Fig. 6b indicates a low dispersion in the relative error, even for $R\,<\,6$. The $Q=0.95$ and $Q=0.98$ quantile show that even when approaching the upper limit of the distribution, high accuracy is maintained.

Fig. 6: Rank dependence of the accuracy of the 2D simulations using thermodynamic tensor models.

a The ECDF of the relative errors on the amplitude of the molar fraction variation plotted in different line colours for TTM rank values $R=3,\ldots ,10$. For $R=6$ and above, markers are also used for better identification of overlapping lines. b Table showing the interquartile region (IQR) and quantiles $Q=0.95$ and $Q=0.98$ of the relative errors for rank values $R=3,\ldots ,10$. c, d Evolution of the volume fraction of the liquid phases as a function of time steps obtained from simulations using directly the CALPHAD model expressions (CE) and using TTMs with rank values $R=3,\ldots ,6$. No further improvement in the accuracy of the simulations using TTMs is obtained when using a rank value higher than $R=6$. It should be noticed, however, that $R=5$ already gives high accuracy. Moreover, for certain applications, lower rank values might already give sufficient accuracy.

As an additional comparison, the volume fraction of one of the phases is measured as a function of time from the simulations. In Fig. 6c, d, results from simulations using the full CALPHAD expression (CE) are plotted along results from simulations using TTMs with rank values $R=3,\ldots ,6$. The volume fractions obtained using a TTM with rank value $R=6$ are almost identical to those obtained with the direct implementation of the CALPHAD expression in the PFM. However, if the main interest is in the volume fraction measurements, a TTM with rank value $R=5$ or even $R=4$ may already give sufficient accuracy (see Fig. 5c, d). Depending on the application and the desired accuracy, a lower rank value $R$ can thus be chosen to limit the number of coefficients in the TTM and the computation time needed to extract data from the tensor by evaluation of equation (3).

In conclusion, exploiting the low-rank tensor character of thermodynamic properties for higher-order multicomponent systems, a highly accurate and compact representation of their composition dependence can be obtained, which can be evaluated easily and efficiently in microstructure evolution simulations. Although tensor decomposition and tensor completion can be applied as a purely data-driven method, we illustrated that there is the flexibility to also account for a priori knowledge on the models or physical phenomena behind the data by applying functional constraints to the factor vectors. It is possible to extend the proposed method to other types of thermodynamic solution models, to include more elements or the dependence on other variables than composition, or to represent different types of composition-dependent phase properties, such as diffusion mobilities, other kinetic coefficients, interface properties, and mechanical properties. The rank value may be slightly higher for more complex dependencies. It is also possible to adjust the degree of the polynomial constraint or even use different types of constraints to model thermodynamic functions with certain features. The statement that the number of coefficients in the CPD depends only linearly on the number of components in the system is an intrinsic property of the CPD and will thus remain valid. The important advantage is that the efficiency and lower complexity is not compromised when adding more elements or variables, which opens up the possibility for detailed higher-order multicomponent microstructure evolution studies, essential for future alloy development.

Methods

CALPHAD thermodynamic model

Ag-Cu-Ni-Sn liquid alloys are considered in this study. According to the CALPHAD method,^10,11,12,45 the required thermodynamic quantities are described assuming a random substitutional solution model. They are a function of $C=4$ molar fractions ${x}_{i}$, with $i=1,\ldots ,C$, and the temperature $T$. The molar Gibbs free energy consists of three parts

$${G}_{{\rm{m}}}={G}^{{\rm{o}}}+{G}_{{\rm{mix}}}^{{\rm{id}}}+{G}_{{\rm{mix}}}^{{\rm{xs}}},$$

(4)

in which ${G}^{{\rm{o}}}$ refers to the Gibbs free energy of a reference surface, ${G}_{{\rm{mix}}}^{{\rm{id}}}$ to the contribution from ideal random mixing (i.e., owing to the configurational entropy), and ${G}_{{\rm{mix}}}^{{\rm{xs}}}$ to the excess energy owing to deviations from ideal behavior. The term ${G}^{{\rm{o}}}$ is given by

$${G}^{{\rm{o}}}=\sum _{i=1}^{C}{x}_{i}{G}_{i}^{{\rm{o}}},$$

(5)

where ${G}_{i}^{{\rm{o}}}$ is the temperature-dependent molar Gibbs free energy of pure component $i$ with respect to the SGTE⁴⁶ reference state. ${G}_{{\rm{mix}}}^{{\rm{id}}}$ is given by

$${G}_{{\rm{mix}}}^{{\rm{id}}}=\bar{R}T\sum _{i=1}^{C}{x}_{i}{\log}{x}_{i},$$

(6)

in which $\bar{R}$ is the ideal gas constant. ${G}_{{\rm{mix}}}^{{\rm{xs}}}$ is expressed using a Redlich–Kister–Muggianu polynomial⁴⁷

$$\begin{array}{ll}{G}_{{\rm{mix}}}^{{\rm{xs}}}=&\sum _{i=1}^{C}\sum _{j> i}^{C}{x}_{i}{x}_{j}\sum _{v}{L}_{ij}^{(v)}{({x}_{i}-{x}_{j})}^{v}\\ &+\sum _{i=1}^{C}\sum _{j> i}^{C}\sum _{k> j}^{C}{x}_{i}{x}_{j}{x}_{k}({x}_{i}{L}_{ijk}^{(0)}+{x}_{j}{L}_{ijk}^{(1)}+{x}_{k}{L}_{ijk}^{(2)}),\end{array}$$

(7)

in which ${L}_{ij}^{(v)}$ are binary interaction parameters and ${L}_{ijk}^{(0,1,2)}$ ternary interaction parameters.

As the molar fractions should always sum to one, i.e., ${\sum }_{i=1}^{C}{x}_{i}=1$, one molar fraction, e.g., the $C$th, is dependent on the other fractions, i.e., ${x}_{C}=1-{\sum }_{n=1}^{N}{x}_{n}$. Therefore, the molar Gibbs free energy can be expressed as a function of $N=C-1$ molar fractions ${x}_{n}$, with $n=1,\ldots ,N$, and the temperature $T$.

The diffusion potential of the $n$th component, needed in multicomponent phase-field simulations, is given by the first partial derivative of the molar Gibbs free energy with respect to ${x}_{n}$,

$$\begin{array}{l}{\tilde{\mu }}_{n}=\frac{\partial {G}_{{\rm{m}}}}{\partial {x}_{n}},\end{array}$$

(8)

which is related to the chemical potentials ${\mu }_{n}$ of component $n$ and ${\mu }_{C}$ of component $C$, as defined in standard solution thermodynamics, as¹⁴

$$\begin{array}{l}{\tilde{\mu }}_{n}={\mu }_{n}-{\mu }_{C}.\end{array}$$

(9)

The partial derivatives of the diffusion potential of the $n$th component with respect to ${x}_{m}$,

$$\displaystyle\frac{\partial {\tilde{\mu }}_{n}}{\partial {x}_{m}}=\frac{\partial {\mu }_{n}}{\partial {x}_{m}}-\frac{\partial {\mu }_{C}}{\partial {x}_{m}}=\frac{{\partial }^{2}{G}_{{\rm{m}}}}{\partial {x}_{n}\partial {x}_{m}},$$

(10)

are often necessary to numerically solve phase-field equations.

Inspection of the different terms in equation (4) shows that only the contributions ${G}^{{\rm{o}}}$ and ${G}_{{\rm{mix}}}^{{\rm{xs}}}$ contain model parameters from the thermodynamic database, namely ${G}_{i}^{{\rm{o}}}$, ${L}_{ij}^{(v)}$, and ${L}_{ijk}^{(0,1,2)}$. As these parameters may not be accessible by the user or the PFM, the terms ${G}^{{\rm{o}}}$ and ${G}_{{\rm{mix}}}^{{\rm{xs}}}$ must be modelled. The contribution ${G}_{{\rm{mix}}}^{{\rm{id}}}$ given by equation (6), however, is a function of known variables only, namely the independent molar fractions ${x}_{n}$ and the temperature $T$. Therefore, the contribution of ${G}_{{\rm{mix}}}^{{\rm{id}}}$ to equation (4) can be calculated at any moment and for any entry of the TDT without consulting the thermodynamic database. As the logarithmic dependence on the molar fractions of this term complicates the formulation of an accurate TTM, we define a starred version of the Gibbs free energy, diffusion potentials, and their partial derivatives with respect to the molar fractions in which the contribution of ${G}_{{\rm{mix}}}^{{\rm{id}}}$ is not included:

$${}^{* }{G}_{{\rm{m}}}={G}_{{\rm{m}}}-\bar{R}{T}\left[\left(\sum _{n=1}^{N}{x}_{n}{\log}{x}_{n}\right)+{x}_{C}{\log}{x}_{C}\right],$$

(11)

$${}^{* }{\tilde{\mu }}_{n}={\tilde{\mu }}_{n}-\bar{R}{T}({\log}{x}_{n}-{\log}{x}_{C}),$$

(12)

$$\frac{{\partial\ }^{* }{\tilde{\mu }}_{n}}{\partial {x}_{m}}=\left\{\begin{array}{ll}\displaystyle\frac{\partial {\tilde{\mu }}_{n}}{\partial {x}_{m}}-\bar{R}T\left(\frac{1}{{x}_{n}}+\frac{1}{{x}_{C}}\right),&n\;=\;m\\ \displaystyle\frac{\partial {\tilde{\mu }}_{n}}{\partial {x}_{m}}-\bar{R}T\left(\frac{1}{{x}_{C}}\right),&n\ne m.\end{array}\right.$$

(13)

Equations (5) and (7) show that ${}^{* }{G}_{{\rm{m}}}$, ${}^{* }{\tilde{\mu }}_{n}$, and $\frac{{\partial\ }^{* }{\tilde{\mu }}_{n}}{\partial {x}_{m}}$ depend polynomially on the molar fractions of the different components. As discussed in Methods section Thermodynamic tensor model, only these polynomial contributions to the TTMs are modelled, knowing that equation (7) for ${G}_{{\rm{mix}}}^{{\rm{id}}}$ and its derivatives can be evaluated and added during the phase-field simulation.

For the simulations, temperature-dependent expressions for all binary and ternary interaction parameters ${L}_{ij}^{(v)}$, ${L}_{ijk}^{(0,1,2)}$, and the Gibbs free energies ${G}_{i}^{{\rm{o}}}$ of the pure elements for the liquid phase are obtained from the COST 531 thermodynamic database.⁴¹

Tensor decomposition

A CPD writes a $N$th-order TDT ${\mathcal{T}}$ as a sum of $R$ rank-1 terms,^28,33,34 each of which is the outer product, denoted by ${\scriptsize \otimes}$, of $N$ nonzero factor vectors ${{\bf{a}}}_{r}^{(n)}$, with $r=1,\ldots ,R$ and $n=1,\ldots ,N$:

$${\mathcal{T}}={\sum _{r=1}^{R}}{\bf{a}}_{r}^{(1)}{\scriptsize\otimes} \cdots {\scriptsize\otimes}\ {\bf{a}}_{r}^{(N)}.$$

(14)

Elementwise, this means that each entry indexed with ${i}_{n}$, for $n=1,\ldots ,N$, of ${\mathcal{T}}$ is given by

$${\mathcal{T}}({i}_{1},\ldots ,{i}_{N})=\sum _{r=1}^{R}{a}_{{i}_{1},r}^{(1)}\ldots {a}_{{i}_{N},r}^{(N)}.$$

(15)

Collecting all factor vectors ${{\bf{a}}}_{r}^{(n)}$ into factor matrices ${{\bf{A}}}^{(n)}=\left[{{\bf{a}}}_{1}^{(n)},\ldots ,{{\bf{a}}}_{R}^{(n)}\right]$, the CPD can be compactly written as

$${\mathcal{T}}=[\kern-2pt[{{\bf{A}}}^{(1)},\ldots ,{{\bf{A}}}^{(N)}]\kern-2pt] .$$

(16)

There exist various techniques to determine the coefficients in the factor matrices based on a limited set of data entries of the original tensor.²⁸ Moreover, when available, a priori knowledge on the model or physical process behind the data can be taken into account, resulting in an even more compact and continuous tensor model of the data. In this work, for example, we exploited the polynomial dependence of ${}^{* }{G}_{{\rm{m}}}$ on the molar fractions, and the fact that the ${}^{* } {\tilde{\mu }}_{n}$ and $\frac{{\partial\ }^{* } \!{\tilde{\mu }}_{n}}{\partial {x}_{m}}$ are the first and second partial derivatives of ${}^{* }{G}_{{\rm{m}}}$. A polynomial constraint could therefore be applied on the factor vectors and the CPDs for ${}^{* }{G}_{{\rm{m}}}$ and its first and second partial derivatives ${}^{* }{\tilde{\mu }}_{n}$ and $\frac{{\partial\ }^{* }{\tilde{\mu }}_{n}}{\partial {x}_{m}}$ could be coupled within a single tensor model, as discussed in Methods Thermodynamic tensor model.

Thermodynamic tensor model

We first show how a TDT can be written as a CPD in which each factor vector ${{\bf{a}}}_{r}^{(n)}$ is a discretized smooth function in a single variable ${x}_{n}$, with $n=1,\ldots ,N$. Moreover, these functions are linear combinations of a set of basis functions common to each factor vector. This way, the problem of modeling all TDTs is recast into finding $N$ coefficient matrices. The method is built upon two mild assumptions: separability and a linear combination of basis functions.

The first assumption is that the underlying function $f({x}_{1},\ldots ,{x}_{N})$ is a sum of separable functions, i.e., the terms can be written as the product of $N$ functions in a single variable:

$$f({x}_{1},\ldots ,{x}_{N})=\sum _{r=1}^{R}{f}_{r}^{(1)}({x}_{1})\cdots {f}_{r}^{(N)}({x}_{N}),$$

(17)

in which $R$ is preferably low. This is exactly the continuous formulation of equation (15). For example, ${G}_{{\rm{m}}}({x}_{1},\ldots ,{x}_{N})$ is not separable with low $R$ because of the term ${G}_{{\rm{mix}}}^{{\rm{id}}}$ (equation (6)), which contains ${x}_{C}\mathrm{log}{x}_{C}$ with ${x}_{C}=1-{\sum }_{n=1}^{N}{x}_{n}$, but ${}^{* }{G}_{{\rm{m}}}({x}_{1},\ldots ,{x}_{N})$ is, as it is the sum of products of polynomials in ${x}_{n}$. As ${G}_{{\rm{mix}}}^{{\rm{id}}}$ (equation (6)) contains only known constants, we can subtract this term and work with ${}^{* }{G}_{{\rm{m}}}({x}_{1},\ldots ,{x}_{N})$ (equation (11)) without loss of generality. If the function $f$, and therefore also each ${f}_{r}^{(n)}({x}_{n})$, is evaluated on (a subset of points of) an $N$-dimensional grid, an (incomplete) tensor is obtained. For example, in this paper, an equidistant grid is defined by ${{\bf{x}}}_{1}=\cdots ={{\bf{x}}}_{N}={\left[\begin{array}{cccc}{\delta }_{x}&2{\delta }_{x}&\ldots &(1-N{\delta }_{x})\end{array}\right]}^{\top }$ for some ${\delta }_{x}$ and the TDT ${}^{\ast}\!{\mathcal{G}}$ is ${}^{* }{G}_{{\rm{m}}}({x}_{1},\ldots ,{x}_{N})$ evaluated on this grid. If we collect the functions ${f}_{r}^{(n)}({x}_{n})$ evaluated in ${{\bf{x}}}_{n}$ as columns of ${{\bf{A}}}^{(n)}$, $n=1,\ldots ,N$, we obtain the CPD:

$${}^{* }{\mathcal{G}}=[\kern-2pt[ {{\bf{A}}}^{(1)},\ldots ,{{\bf{A}}}^{(N)}]\kern-2pt],$$

(18)

where ${}^{* }{\mathcal{G}}$ is the tensor containing the values ${}^{* }{G}_{{\rm{m}}}({x}_{1},\ldots ,{x}_{N})$. Similarly, the TDTs for ${}^{* }{\tilde{\mu }}_{n}$ and $\frac{{\partial }^{* }{\tilde{\mu }}_{n}}{\partial {x}_{m}}$ can be written as a CPD. Note that the grid points do not need to be chosen in the same way for different components $n$, nor do they have to be equidistant.

Second, we assume that each ${f}_{r}^{(n)}$ can be expressed as a linear combination of $K$ basis functions ${b}_{k}^{(n)}({x}_{n})$, $k=1,\ldots ,K$:

$$\begin{array}{ll}{f}_{r}^{(n)}({x}_{n})&={u}_{1,r}^{(n)}{b}_{1}^{(n)}({x}_{n})+\cdots +{u}_{K,r}^{(n)}{b}_{K}^{(n)}({x}_{n})\\ &=\left[\begin{array}{ccc}{b}_{1}^{(n)}({x}_{n})&\ldots &{b}_{K}^{(n)}({x}_{n})\end{array}\right]{\left[\begin{array}{ccc}{u}_{1,r}^{(n)}& \ldots & {u}_{K,r}^{(n)}\end{array}\right]}^{\top }.\end{array}$$

(19)

By evaluating ${f}_{r}^{(n)}({x}_{n})$ on the grid, we obtain

$${{\bf{A}}}^{(n)}=\left[\begin{array}{ccc}{{\bf{b}}}_{1}^{(n)}&\ldots &{{\bf{b}}}_{K}^{(n)}\end{array}\right]{{\bf{U}}}^{(n)}={{\bf{B}}}^{(n)}{{\bf{U}}}^{(n)},$$

(20)

in which ${{\bf{B}}}^{(n)}$ collects the basis functions evaluated in ${{\bf{x}}}_{n}$ and ${{\bf{U}}}^{(n)}$ the coefficients. As the basis functions are fixed, only the coefficients ${{\bf{U}}}^{(n)}$ need to be determined. In this paper, we choose a basis of monomial functions evaluated in the same points in all dimensions, i.e., ${b}_{k}^{(n)}({x}_{n})={b}_{k}(x)={x}^{(k-1)}$ for $k=1,\ldots ,d+1$, and $n=1,\ldots ,N$, which is reasonable given the polynomial nature of the theoretical models; see equations (11), (12), and (13). This way, each TDT can be written as a CPD in which each factor matrix ${{\bf{A}}}^{(n)}$ is subject to a linear constraint:

$${}^{\ast}\!{\mathcal{G}}= [\kern-2pt[{\bf{B}}{{\bf{U}}}^{(1)},\ldots ,{\bf{B}}{{\bf{U}}}^{(N)}]\kern-2pt] \quad {\rm{with}}\quad {\bf{B}}=\left[\begin{array}{ccccc}{\bf{1}}&{\bf{x}}&{{\bf{x}}}^{2}&\ldots & {\bf{x}}^{d}\end{array}\right]$$

(21)

in which the powers are elementwise. Note that the coefficient matrices ${{\bf{U}}}^{(n)}$ do not necessarily contain the coefficients ${G}_{i}^{o}$, ${L}_{ij}^{(v)}$, and ${L}_{ijk}^{(0,1,2)}$ as their entries.

A key property following from the definition of the TDTs and the CPD formulation with linear constraints is that computing derivatives only requires the computation of derivatives of the known basis functions ${b}_{k}({x}_{n})$:

$$\frac{{\rm{d}}{f}_{r}^{(n)}({x}_{n})}{{\rm{d}}{x}_{n}}={u}_{1,r}^{(n)}\frac{{\rm{d}}{b}_{1}({x}_{n})}{{\rm{d}}{x}_{n}}+\cdots +{u}_{K,r}^{(n)}\frac{{\rm{d}}{b}_{K}({x}_{n})}{{\rm{d}}{x}_{n}}.$$

(22)

For example, as ${}^{* }{\tilde{\mu }}_{1}=\frac{{\partial }^{* }{G}_{{\rm{m}}}}{\partial {x}_{1}}$, the corresponding TDT ${}^{* }{{\mathcal{P}}}^{(1)}$ can be represented as

$${}^{\ast}\!{\mathcal{P}}^{(1)}=[\kern-2pt[ {\dot{{\bf{A}}}}^{(1)},{\bf{B}}{{\bf{U}}}^{(2)},\ldots ,{\bf{B}}{{\bf{U}}}^{(N)}]\kern-2pt] \quad {\rm{with}}\quad {\dot{{\bf{A}}}}^{(1)}=\dot{{\bf{B}}}{{\bf{U}}}^{(1)},$$

(23)

as only ${{\bf{A}}}^{(1)}$ depends on ${x}_{1}$. $\dot{{\bf{B}}}$ contains the evaluated derived functions ${\dot{b}}_{k}(x)=(k-1){x}^{k-1}$:

$$\dot{{\bf{B}}}=\left[\begin{array}{ccccc}\bf{0} & \bf{1} & {2}\bf{x} & \ldots & {d}{\bf{x}}^{d-1}\end{array}\right].$$

(24)

Similarly, the derivatives of the potentials depend on $\dot{{\bf{B}}}$ and $\ddot{{\bf{B}}}$, in which the latter matrix contains the second-order derivatives:

$$\ddot{{\bf{B}}}=\left[\begin{array}{cccccc}\bf{0}&\bf{0}&{2}\times\bf{1}&6\bf{x}&\ldots& d(d-1){\bf{x}}^{d-2}\end{array}\right].$$

(25)

For example, the CPDs of the TDTs corresponding to $\frac{{\partial\ }^{* }{\tilde{\mu }}_{1}}{\partial {x}_{1}}$ and $\frac{{\partial\ }^{* }{\tilde{\mu }}_{1}}{\partial {x}_{2}}$ are given by

$${}^{* }{{\mathcal{D}}}^{(1,1)}= \left[\kern-2pt\left[\ddot{{\bf{B}}}{{\bf{U}}}^{(1)},{\bf{B}}{{\bf{U}}}^{(2)},\ldots ,{\bf{B}}{{\bf{U}}}^{(N)}\right]\kern-2pt\right],$$

(26)

$${}^{* }{{\mathcal{D}}}^{(1,2)}= \left[\kern-2pt\left[\dot{{\bf{B}}}{{\bf{U}}}^{(1)},\dot{{\bf{B}}}{{\bf{U}}}^{(2)},\ldots ,{\bf{B}}{{\bf{U}}}^{(N)}\right]\kern-2pt\right],$$

(27)

respectively. Hence, once the coefficient matrices ${{\bf{U}}}^{(n)}$, $n=1,\ldots ,N$, are determined, all TDTs are modeled. Moreover, as the coefficients are independent of the grid, different grids can be used in the training phase, i.e., to find ${{\bf{U}}}^{(n)}$ and when actually using the model, i.e., during the phase-field simulations.

Thermodynamic tensor model training

In the optimisation routine for modeling the TDTs of a quaternary system ($N=3$), we want to find the coefficients in the matrices ${{\bf{U}}}^{(1)}$, ${{\bf{U}}}^{(2)}$, ${{\bf{U}}}^{(3)}$ that minimise the least squares cost function

$$\begin{array}{ll}{{\rm{min}}}\atop{{{\bf{U}}}^{(1)},{{\bf{U}}}^{(2)},{{\bf{U}}}^{(3)}}\,&\frac{1}{2}{\left\Vert {\mathcal{S}}* \left({\,}^{*}{{\mathcal{P}}}_{{\rm{T}}}^{(1)}- \left[\kern-2pt\left[\dot{{\bf{B}}}{{\bf{U}}}^{(1)},{\bf{B}}{{\bf{U}}}^{(2)},{\bf{B}}{{\bf{U}}}^{(3)}\right]\kern-2pt\right] \right)\right\Vert }^{2}\\ &+\frac{1}{2}{\left\Vert {\mathcal{S}}* \left({\,}^{*}{{\mathcal{P}}}_{{\rm{T}}}^{(2)}-\left[\kern-2pt\left[{\bf{B}}{{\bf{U}}}^{(1)},\dot{{\bf{B}}}{{\bf{U}}}^{(2)},{\bf{B}}{{\bf{U}}}^{(3)}\right]\kern-2pt\right] \right)\right\Vert }^{2}\\ &+\frac{1}{2}{\left\Vert {\mathcal{S}}* \left({\,}^{*}{{\mathcal{P}}}_{{\rm{T}}}^{(3)}- \left[\kern-2pt\left[{\bf{B}}{{\bf{U}}}^{(1)},{\bf{B}}{{\bf{U}}}^{(2)},\dot{{\bf{B}}}{{\bf{U}}}^{(3)}\right]\kern-2pt\right]\right)\right\Vert }^{2}.\end{array}$$

(28)

${}^{* }{{\mathcal{P}}}_{{\rm{T}}}^{(1)}$, ${}^{* }{{\mathcal{P}}}_{{\rm{T}}}^{(2)}$, and ${}^{* }{{\mathcal{P}}}_{{\rm{T}}}^{(3)}$ are training TDTs, which are computed with a step size of ${\delta }_{x}=0.025$ using the TC-Toolbox (Thermo-Calc 2017b) with the COST 531 thermodynamic database⁴¹ at 1400 K, each with 11,482 entries. The residual between the tensor and its low-rank, polynomial-constrained CPD model, e.g., ${}^{* }{{\mathcal{P}}}_{{\rm{T}}}^{(1)}-[\kern-2pt[\dot{{\bf{B}}}{{\bf{U}}}^{(1)},{\bf{B}}{{\bf{U}}}^{(2)},{\bf{B}}{{\bf{U}}}^{(3)}]\kern-2pt]$, is multiplied elementwise (denoted by $*$) with a binary tensor ${\mathcal{S}}$ which is one for each computed training entry and zero otherwise, such that only the limited set of training entries contribute to the objective function. ${\mathcal{S}}$ can be different for each term of the objective problem during the optimisation, but here we choose to keep it constant. Only the diffusion potentials are used in the optimisation routine since other approaches, e.g., using ${}^{\ast}\!{\mathcal{G}}$ or a combination of all TDTs, have shown to be less accurate. Therefore, as only derivatives are used, the Gibbs free energy is determined up to a constant $\alpha$, which can be estimated as the mean of the residual and must be added to the computed CPD for the Gibbs free energy.

In this work, a conservative, i.e., small enough, value for ${\delta }_{x}$ is chosen when selecting the training data to guarantee a highly accurate model. No further investigations are conducted on the influence of the step size used for the training data on the accuracy of the tensor model or phase-field simulations. Various other aspects can affect the accuracy of the tensor model, such as the degree of the polynomial constraint, the rank value of the CPD and the distribution of the training data over the system. However, the grid with which data are sampled is not very relevant, as we are interested in the optimisation of the coefficient in the ${{\bf{U}}}^{(n)}$ matrices and not the factors in the ${{\bf{A}}}^{(n)}$ matrices. The factors in the latter are only computed after the optimisation routine with equation (20) and the grid spacing can be chosen according to the necessity. This procedure also assures that the thermodynamic quantities are accurately modelled over the entire composition domain (including the dilute solution regions, where at least one of the molar fractions takes a very low value).

The least squares cost function is optimised in Tensorlab⁴⁰ using the data independent algorithm for the CPD of incomplete tensors subject to linear constraints (CPDLI DI).³¹ A few decompositions are computed starting from random initial guesses for ${{\bf{U}}}^{(1)}$, ${{\bf{U}}}^{(2)}$, ${{\bf{U}}}^{(3)}$, given $R$ and $d$, and the best result is retained. CPDs for rank values $R=3,\ldots ,10$ are computed. (Models with maximal degrees per variable $d=3,4,5$ have been computed as well to validate that the correct degree $d=4$ can indeed be recovered. This degree $d=4$ can be obtained from equation (4) and the parameters obtained from the thermodynamic database COST 531.⁴¹)

Once the matrices with the unknown coefficients ${{\bf{U}}}^{(1)}$, ${{\bf{U}}}^{(2)}$, ${{\bf{U}}}^{(3)}$ are found, the factor matrices ${{\bf{A}}}^{(n)}$, ${\dot{{\bf{A}}}}^{(n)}$ and ${\ddot{{\bf{A}}}}^{(n)}$, for $n=1,2,3$, are constructed using equation (20) and basis matrices ${\bf{B}}$, $\dot{{\bf{B}}}$, $\ddot{{\bf{B}}}$, which are evaluated on a grid with ${\delta }_{x}=0.0001$.

Validation

The TC-toolbox (Thermo-Calc 2017b) and MATLAB are used to generate validation TDTs (indicated with subscript ${\rm{V}}$) of the Gibbs free energy ${{\mathcal{G}}}_{{\rm{V}}}$, diffusion potentials ${{\mathcal{P}}}_{{\rm{V}}}^{(1)}$, ${{\mathcal{P}}}_{{\rm{V}}}^{(2)}$, and ${{\mathcal{P}}}_{{\rm{V}}}^{(3)}$ and of the derivatives of the diffusion potentials ${{\mathcal{D}}}_{{\rm{V}}}^{(1,1)}$, ${{\mathcal{D}}}_{{\rm{V}}}^{(2,2)}$, ${{\mathcal{D}}}_{{\rm{V}}}^{(3,3)}$, ${{\mathcal{D}}}_{{\rm{V}}}^{(1,2)}$, ${{\mathcal{D}}}_{{\rm{V}}}^{(2,3)}$, and ${{\mathcal{D}}}_{{\rm{V}}}^{(1,3)}$ of the liquid phase at 1400 K, with a step size of ${\delta }_{x}=0.0001$ and using the COST 531 thermodynamic database.⁴¹

The accuracy of the TTMs with rank value $R=3,\ldots ,10$ is verified based on the RE on the entries of the validation TDTs. For example, for the Gibbs free energy tensor, the relative error ${E}_{R}^{({\mathcal{G}})}$ given by

$${E}_{R}^{({\mathcal{G}})}=\left| \frac{{{\mathcal{G}}}_{{\rm{V}}}({i}_{1},{i}_{2},{i}_{3})-\left[\sum _{r=1}^{R}{a}_{{i}_{1},r}^{(1)}\ {a}_{{i}_{2},r}^{(2)}\ {a}_{{i}_{3},r}^{(3)}+\bar{R}T\left(\sum _{i=1}^{N}{x}_{i}{\log}{x}_{i}+{x}_{C}{\log}{x}_{C}\right)+\alpha \right]}{{{\mathcal{G}}}_{{\rm{V}}}({i}_{1},{i}_{2},{i}_{3})}\right| .$$

(29)

The contribution from the ideal mixing term ${G}_{{\rm{mix}}}^{{\rm{id}}}$ is calculated separately and added to the values calculated from the TTMs, which only include the ${}^{* }{G}_{{\rm{m}}}$ in equation (11) (see Methods section CALPHAD thermodynamic model). As discussed in Methods Thermodynamic tensor model training, a constant value $\alpha$ has to be added to the Gibbs energy values computed from the TTMs. Analogously, the relative error ${E}_{R}^{({\mathcal{P}})}$ of the diffusion potentials and ${E}_{R}^{({\mathcal{D}})}$ of the partial derivatives of the diffusion potentials with respect to the molar fractions are calculated including their ideal mixing contributions: see equations (12) and (13). As $\alpha$ is a constant, it does not appear in the expressions for the diffusion potentials and their partial derivatives with respect to the molar fractions.

An ECDF is then determined for each rank value $R$ considering all the values obtained for ${E}_{R}^{({\mathcal{G}})}$, ${E}_{R}^{({\mathcal{P}})}$, and ${E}_{R}^{({\mathcal{D}})}$ together. The MATLAB fuction ecdf is used for this. The ECDFs are plotted in Fig. 3a; the IQR, and quantiles $Q=0.95$ and $Q=0.98$ are presented in Fig. 3b. Additional quantile values can be found in the Supplementary Information.

The validation of the 2D simulations using TTMs is conducted based on a comparison with simulations performed using the direct implementation of the CALPHAD model (approach b in Fig. 4). The reason for this choice is the fact that most of the alternatives, presented in Fig. 4, cannot be applied to a quaternary system owing to certain limitations. For instance, in the approach that uses polynomial fitting (Fig. 4c), only a narrow region of the thermodynamic function can be accurately described and since diffusion paths can become complex for multicomponent systems, it is difficult to identify which composition region will be most relevant for a simulation. Moreover, polynomial fits cannot prevent molar fractions from taking an unphysical value below zero or above one. Although this can be avoided in binary system, by designing the simulation in such a way that the molar fraction does not take a negative value, when considering four or more elements, there is always a time or location in the simulation where one of the molar fraction fields goes below zero, making the remainder of the simulation unstable. The approach that uses lookup tables (Fig. 4f) is also unfeasible in this case owing to the size of the sampled data sets given the required high resolution, as seen in (Fig. 1a). These approaches, however, are suitable and straightforward to be applied to binary or ternary systems, whereas the proposed use of tensor completion requires at least three independent variables in the thermodynamic descriptions.

The approach that uses the direct implementation (Fig. 4b) of the CALPHAD model expression and the one that uses an interface with a thermodynamic software (Fig. 4d) can be applied for quaternary and higher-order systems. With our implementation of the method that uses an interface with a thermodynamic software we could only run 1D simulations in a reasonable amount of time (these are presented in the Supplementary Information). It should be mentioned that optimisations (mainly based on the idea of interpolation between previously calculated data) of phase-field codes using an interface with a thermodynamic software have been reported in the literature;^18,19 however, to the best of our knowledge, only for binary and ternary systems.

For the method using the direct implementation of the CALPHAD model expression (CE), we found that (for our implementation) the efficiency is similar to that using the TTM for the quaternary system. However, such a comparison is often not relevant. The direct implementation of the CALPHAD model expression requires access to an open thermodynamic database (while today, the most advanced thermodynamic databases are commercial). Moreover, if the phase under study is modelled using sublattices, it is mostly impossible (except for some special cases) to relate the phase-field molar fractions with the site fractions in the CALPHAD expression. Hence, in these cases, the direct implementation of the CALPHAD model expression is not even possible, whereas the tensor model can still be formulated. We also observed that the current implementation of the TTM in the phase-field model is memory bound, which can be addressed with further optimisation studies. Theoretically, the TTM method should be cheaper to evaluate than CE.

For the validation of the 2D simulation, the LBM method is used to estimate the amplitudes of the molar fraction distributions of the simulated microstructures during spinodal decomposition and consequent coarsening.^42,43,44 For each simulation, the microstructures obtained at every thousandth time step are considered. For each microstructure, a Gaussian mixture model is used to fit the distributions of the molar fraction values for each of the three independent components with the MATLAB function gmdistribution. A bimodal distribution is obtained with two means ${m}_{1}$ and ${m}_{2}$, corresponding to the compositions of the two liquid phases, as shown in Fig. 5b–d. The amplitude $A$ is then calculated as the difference between these two means, i.e., $A=| {m}_{1}-{m}_{2}|$. The amplitude ${A}^{R}$ from the simulation using the TTM with rank value $R$ is compared with the amplitude $A$ obtained from the simulation using the direct implementation of the CALPHAD model (approach b in Fig. 4). The relative error ${E}_{R}^{{\rm{2D}}}$ is defined by

$${E}_{R}^{{\rm{2D}}}= \left| \frac{A-{A}^{R}}{A}\right| ,$$

(30)

and is determined for each molar fraction field and each alloy composition and for rank values $R=3,\ldots ,10$. Then, the ECDF is calculated for each rank value $R$, considering all the RE determined for the different alloy compositions and the three independent molar fraction fields, obtained for the given rank value $R$, together. The ECDFs are plotted in Fig. 6a and the information of the IQR, $Q=0.95$ and $Q=0.98$ are presented in Fig. 6b.

In addition, the volume fractions of the two liquid phases are measured at every thousandth time step as follows. First, the mean of one of the molar fraction distributions is computed. Then, the grid cells of the simulation domain are divided into two groups, one containing the grid cells with a molar fraction value greater than the mean and the other with the grid cells with a molar fraction value smaller than the mean. The volume fractions of the two phases are then obtained by dividing the number of grid cells in each group by the total number of cells in the simulation domain. This information is presented in Fig. 6c, d.

Phase-field microstructure evolution model

Assuming a volume-fixed frame of reference ⁴⁸, the temporal and spatial evolution of the non-equilibrium microstructure of the quaternary Ag-Cu-Ni-Sn system ($C=4$), is described with $N=C-1=3$ conserved field variables, namely the molar fractions ${x}_{1}({\bf{r}},t)$, ${x}_{2}({\bf{r}},t)$, and ${x}_{3}({\bf{r}},t)$ of silver, copper, and nickel, respectively. The molar fraction of tin is then obtained as ${x}_{C}=1-{x}_{1}-{x}_{2}-{x}_{3}$.

Following the Cahn-Hilliard⁴⁹ formalism, the non-uniform free energy functional is defined as

$$F({x}_{1},{x}_{2},{x}_{3})={\int }_{V}\left[{f}_{0}({x}_{1},{x}_{2},{x}_{3})+\sum _{n}^{N}\frac{{k}_{n}}{2}{(\nabla {x}_{n})}^{2}\right]{\rm{d}}{\bf{r}},$$

(31)

where ${f}_{0}({x}_{1},{x}_{2},{x}_{3})$ is the homogeneous free energy density and ${k}_{n}$ are positive gradient energy coefficients for $n=1,2,3$. The variational derivative of $F$ with respect to ${x}_{n}$ is given by

$$\frac{\delta F}{\delta {x}_{n}}=\frac{\partial {f}_{0}({x}_{1},{x}_{2},{x}_{3})}{\partial {x}_{n}}-{k}_{n}{\nabla }^{2}{x}_{n}.$$

(32)

The partial derivative of the homogeneous free energy is obtained given the thermodynamic relation¹⁴ in equation (9) as

$$\frac{\partial {f}_{0}}{\partial {x}_{n}}=\frac{1}{{V}_{{\rm{m}}}}({\mu }_{n}-{\mu }_{C})=\frac{{\tilde{\mu }}_{n}}{{V}_{m}},$$

(33)

with ${\tilde{\mu }}_{n}$ the diffusion potential and ${\mu }_{n}$ the chemical potential of component $n$ ($n=1,2,3$) and ${\mu }_{C}$ the chemical potential of Sn.

The diffusion equation describes the evolution of the conserved field variables for the multicomponent system⁷ and is given by

$$\left(\frac{1}{{V}_{{\rm{m}}}}\right)\frac{\partial {x}_{n}({\bf{r}},t)}{\partial t}=-\nabla \times {{\bf{J}}}_{n},$$

(34)

in which ${V}_{{\rm{m}}}$ is a constant coefficient approximating the molar volume of the alloy. ${{\bf{J}}}_{n}$ is the diffusion flux of the $n$th component, which is defined with the Onsager reciprocal relations⁵⁰ as

$${{\bf{J}}}_{n}=\sum _{m}^{N}-{M}_{nm}\nabla \frac{\delta F}{\delta {x}_{m}},$$

(35)

in which ${M}_{nm}$ is the mobility of component $n$ taken with component $m$ as reference. As information on the diffusion coefficients for quaternary systems are scarce, equal atomic mobilities $\beta$ for all components are adopted, giving⁴⁸

$${M}_{nm}=\left\{\begin{array}{ll}\displaystyle\frac{1}{{V}_{{\rm{m}}}}\beta {x}_{n}(1-{x}_{m}),&n=m\\ -\frac{1}{{V}_{{\rm{m}}}}\beta {x}_{n}{x}_{m},&n\;\ne\;m.\end{array}\right.$$

(36)

Equation (34) is implemented in MATLAB for 1D and 2D simulations using finite difference discretization with explicit Euler time stepping and a central three-point (1D) or five-point (2D) scheme for the spatial coordinates. The parameters β = 10⁻¹² m⁵/(Js), Vm = 10⁻⁵ m³/mol, k_n = 5 × 10⁻¹¹ J/m are used. A domain size of 5 × 10⁻⁸ m for the 1D and 25 × 10⁻¹⁶ m² for the 2D simulations is considered with a grid spacing of 2 × 10⁻¹⁰ m. A time step of 10⁻¹³ s is used. For each simulation, an evolution time of $1{0}^{5}$ simulation time steps, equivalent to10⁻⁸ s, is considered.