Extended Regression Analysis for Debye–Einstein Models Describing Low Temperature Heat Capacity Data of Solids

Heat capacity data of many crystalline solids can be described in a physically sound manner by Debye–Einstein integrals in the temperature range from 0K to 300K. The parameters of the Debye–Einstein approach are either obtained by a Markov chain Monte Carlo (MCMC) global optimization method or by a Levenberg–Marquardt (LM) local optimization routine. In the case of the MCMC approach the model parameters and the coefficients of a function describing the residuals of the measurement points are simultaneously optimized. Thereby, the Bayesian credible interval for the heat capacity function is obtained. Although both regression tools (LM and MCMC) are completely different approaches, not only the values of the Debye–Einstein parameters, but also their standard errors appear to be similar. The calculated model parameters and their associated standard errors are then used to derive the enthalpy, entropy and Gibbs energy as functions of temperature. By direct insertion of the MCMC parameters of all 4·105 computer runs the distributions of the integral quantities enthalpy, entropy and Gibbs energy are determined.


Introduction
Some sets of model parameters used to fit low temperature heat capacity measurements are purely empirical (see, e.g., [1][2][3][4]), other parameters are at least partly motivated by theory, e.g., following the Debye-Einstein approach [5][6][7][8][9][10][11][12].A Debye-Einstein integral to describe heat capacities of several compounds is introduced by Kelley and King [5].They proposed a theory-based heat capacity function that contains two parts, where a p-atomic isotropic crystal consists of one Debye term and (p − 1) Einstein terms.The value p equals the number of atoms in a molecule or else the number of atoms in the simplest chemical formula that may be written to represent the composition.A modification of this heat capacity function can be found in Wu et al. [10].It is demonstrated in [12] that the temperature-dependent standard molar heat capacities (C o p,m (T)) data for many crystalline solids can be described by means of these Debye-Einstein integrals over the low temperature range (0-300 K).
The Debye-Einstein regression analysis only requires a small number of thermodynamically motivated fitting parameters.Unlike using fitting polynomials or splines, it is possible to extrapolate the Debye-Einstein model to zero Kelvin, even if experimental C o p,m data are lacking for ultra-low temperatures, e.g., below 50 K (details can be found in [9]).Compared to various models (e.g., polynomial fits) that require a large number of fit parameters over the range of 0-300 K, the Debye-Einstein approach offers the advantage of easy tabulation of a few (four to six) fit parameters in a systematic manner for different crystalline solids.Moreover, models with a high number of fit parameters often encounter the issue of overfitting, which becomes evident when the uncertainties associated with these fit parameters are comparable to or even exceed the magnitude of the parameters themselves.In general, C o p,m values increase strongly, ranging from very small values close to 0 K to increasingly high values at T = 298.15K.In addition, the C o p,m values increase with the number of atoms in the formula unit.This possible wide range of variation indicates that the standard errors associated with the residuals also change with temperature.
Thermodynamic data are often fitted by means of various simply empirical or thermodynamically motivated models, where the uncertainties of the fit parameters are not provided (see, e.g., [13]).Frequently, the correlation of the fit parameters is missing too.However, uncertainty estimation is now the focus of recently published papers, e.g., in [14], where experimental data and atomistic experiments are included in the model development.In general, thermodynamic modeling relies on the significance of the experimental data, their availability and on the adequacy of the model functions.In this context, it is pointed out in Honarmandi et al. [15] that uncertainty quantification of phase diagrams is of paramount importance for decision making in materials design.Paulson et al. [16] state that uncertainty quantification in combination with CALPHAD is not yet widely adopted, and present uncertainty quantification of the properties from CALPHAD modeling and make their program codes available.Uncertainty quantification in thermodynamic modeling follows either from classical (frequentist) statistics (e.g., [17]) or when Bayesian statistics (e.g., [15,16,18]) are applied.One advantage of Bayesian interference is that prior knowledge can be introduced into the calculations updated by the likelihood function, which is influenced by new measurements in order to obtain the posterior probability.A direct calculation of the posterior probability in the n-dimensional problem set (n parameters are sought) is almost impossible.However, by means of modern data sampling techniques using the Markov chain Monte Carlo (MCMC) method the posterior probability distribution is approximated (see, e.g., Vrugt and Ter Braak [19]).The MCMC method leads to a comparatively fast convergence of high-dimensional problems and is, therefore, becoming increasingly popular.
In this work, the estimated distribution of the parameters of the Debye-Einstein fit are obtained from the local Levenberg-Marquardt least squares minimization and from the global MCMC regression tool.In the case of MCMC regression the standard errors of the experimental heat capacities are described as a function with parameters derived from the residuals (i.e., the deviations of the measured values to the calculated values at the measured heat capacities).In this way, unknown systematic errors are indirectly introduced into the error calculations [20].Unknown systematic errors due to a specific model occur since all models are only a simplification of reality [21].The parameters of the function describing the standard errors are simultaneously optimized with the parameters of the Debye-Einstein heat capacity function.Thereby, the Bayesian credible interval is calculated for the temperature-dependent heat capacities C o p,m (T).
The standard errors of the parameters are calculated by means of the classical least squares method and the best estimates of the fit parameters are obtained by using the global MCMC method.It is shown that the estimated distributions of the parameters calculated from the MCMC regression fit well to the standard errors which follow from the classical least squares method.

Literature Data and Theory
Calorimetric measurements of heat capacity data for various crystalline solids have been extensively reported in the literature (see e.g., [2,10,[22][23][24][25][26]).These measurements cover a temperature range from ultra-low temperatures, such as 2 K, up to 300 K or slightly higher temperatures.The heat capacity data are commonly obtained using the relaxation method, e.g., provided by Quantum Design Physical Property Measurement Systems, San Diego, CA, USA [27].The uncertainties associated with the measurements, including instrumental errors and statistical fluctuations, were rigorously evaluated and reported in [28] and usually come along with a relative uncertainty σ C p /C P = 5 • 10 −3 for T < 50 K and σ C p /C P = 3 • 10 −3 for T > 50 K.
In the following, the Debye-Einstein integral fit of heat capacity data of SrMoO 4 is presented, which is described by six parameters only in the low temperature range (0-320 K).The underlying experimental data are published in Morishita and Houshiyama [26].

Debye-Einstein Integral
In our study, we employ the Debye-Einstein integral fit (Equation (1)) proposed by Wu et al. [10] and extend its application to the entire measurement range (see also Gamsjäger and Wiessner [12] and Ogris and Gamsjäger [29]): where the Debye integral D(T D /T) is given by: and the two Einstein terms E 1 (T E1 /T) and E 2 (T E2 /T) are given by: These equations involve the following fit parameters: the Debye temperature T D , the Einstein temperatures T Ei with i = 1 or i = 2, and the prefactors m, n 1 and n 2 .According to theory, the sum of the prefactors are equal to the number of atoms in the formula unit as can be found in [5].Thus, the low temperature C o p,m measurements are described by the Debye-Einstein model.Starting from the temperature-dependent standard molar heat capacity C o p,m (T), the molar enthalpy H o m (T) and the molar entropy S o m (T) are obtained by integration as follows: and It is worth noting that the Debye-Einstein integral approach allows for extrapolation to absolute zero in the case that ultra-low temperature data are missing; e.g., data are compiled for T > 50 K only (see, e.g., [9,29]).The prefactors m, n 1 and n 2 describe the weight of the Debye integral and the Einstein terms, respectively, where the sum is not a priori fixed in the fitting algorithm, but appears to be close to the number of atoms in the formula unit of the investigated compound as it should be from a theoretical point of view (see also [12]).This is an indication that the fit parameters used in the Debye-Einstein integral are not only of an empirical nature, but are relevant with respect to the theory behind.In addition to the Levenberg-Marquardt least squares calculations for finding the optimized parameters of the Debye-Einstein integral, the optimal values of these fit parameters within their distributions are estimated by means of the MCMC method within a Bayesian framework.

Bayesian Framework and MCMC Regression
Unlike traditional optimization methods that rely solely on minimizing the sum of least squares of the residuals, the Bayesian approach offers a probabilistic framework for incorporating prior knowledge, estimating model parameters and quantifying uncertainties.
Bayes' theorem is derived from the product rule of conditional probabilities (see, e.g., [30]).The conditional probability P(B|A) represents the probability of event B given that event A is true.Applying Bayes' theorem for our case, the posterior probability P ⃗ ξ D, H follows from the following Equation ( 6): The prior probability distribution P ⃗ ξ H incorporates any existing prior knowledge about the vector of the parameters ⃗ ξ within the hypothesis space H.
The prior distribution is updated with new data D by using the likelihood function In our investigation, the vector ⃗ ξ consists of the six parameters of the Debye-Einstein integral, i.e., the prefactors m, n 1 , n 2 , the Debye temperature T D and the Einstein temperatures T E1 and T E2 , and the parameters s 0 and s 1 of the function of the standard errors.It is assumed that the hypothesis space H remains constant, which implies that the distribution P(D|H), commonly referred to as evidence, also remains constant.In our case, boundaries are imposed on the hypothesis space H, since the model parameters have to be positive; The Bayesian equation is often transformed into its logarithmic form due to numeric advantages: For the likelihood function in Equation (6) or Equation (7), it is commonly used to employ the Gaussian distribution Gauss(y i , y c,i , σ i ) to calculate the probability for each measurement point.
Furthermore, in a "naive" manner, it is commonly assumed that the residuals, i.e., the measured heat capacities minus the calculated heat capacity values (y i − y c,i ) for all i data points, are independent of each other.Therefore, the likelihood function is obtained by multiplying the individual probabilities of each Gaussian distribution.In this Bayesian framework, the standard errors of the residuals are estimated.To reduce the number of fit parameters to a manageable level, we propose to apply a simple function with the parameters ⃗ ξ 2 that describes these standard errors.The standard errors are influenced by the uncertainties in the measurements and unknown errors brought in by the model.The logarithmic posterior distribution that accounts for the function of the standard errors is written as Assuming a flat prior (no prior knowledge), the posterior distribution simplifies to Equation (10): The Markov chain Monte Carlo (MCMC) sampling method can be used effectively within the Bayesian framework.MCMC allows us to explore the parameter space with the aim to eventually converge to the joint posterior probability.Thereby, robust estimates of the parameters are possible and their associated standard errors are obtained from the posterior probability density distributions.
In our study, we employ an advanced version of the Metropolis-Hastings algorithm, known as the Differential Evolution Adaptive Metropolis (DREAM) algorithm, which was initially developed by Braak [31] and further enhanced by Vrugt and Ter Braak [19].The DREAM algorithm is specifically designed for Bayesian optimization and incorporates multiple chains with differential evolution and adaptive Metropolis-Hastings steps.This MCMC approach substantially enhances the exploration of the parameter space by dynamically adjusting the step sizes, leading to improved convergence and efficiency in the optimization process.The relative frequencies of parameter occurrences within the parameter range directly correspond to their probability density distribution.The standard errors of the parameters equal the standard deviations of these parameters and are calculated by considering all values from all Markov chains.These MCMC-based standard errors can be compared with the standard errors obtained from the error propagation rule from classical statistics.
In the following, the correlations between the parameters are estimated.The covariance between two parameters can be calculated using the following formula: Here, N represents the total number of samples or observations of all Markov chains, A i and B i are the values of parameters A and B for the ith element of the Markov chains, and Ā and B denote the mean values of the parameters A and B of all entries in the Markov chains, respectively.
The correlation coefficient, denoted as r, is defined as: In this equation, σ A and σ B represent the standard errors of parameters A and B, respectively.

Results and Discussion
The molar heat capacities C o p,m of SrMoO 4 have been measured over a temperature range from 2 K to 320 K by Morishita and Houshiyama [26] using a relaxation method instrument.As an example, these data, i.e., 81 C o p,m (T) data pairs, are evaluated in this work by means of the Debye-Einstein approach using both methods, least squares minimization and Bayesian statistics, with the help of Monte Carlo Markov chains (MCMC).Regression by the latter method is based on the DREAM algorithm.For the analyses, we used 10 chains, each consisting of 5 • 10 4 iterations, with the initial 1 • 10 4 iterations per chain discarded as burn-in.This means that a total of 4 • 10 5 parameter sets were available for analysis.
The probability density distributions of the simulated heat capacities follow from the MCMC approach by fitting the experimental heat capacities.These probability density distributions are presented at selected temperatures.It is worth noting that these distributions, that follow from the six-parameter Debye-Einstein integral, can be extrapolated to lower temperatures in the case of lacking experimental data.The probability density distribution of C o p,m at T = 15.0K is presented in Figure 1a, the probability density distribution of C o p,m at T = 98.1 K is shown in Figure 1b, and the probability density distributions at T = 248.6K and T = 298.15K can be seen in Figure 1c,d, respectively.The probability density distribution of C o p,m at T = 248.6K (Figure 1c) exhibits two distinct maxima.It can be speculated that these two maxima occur due to correlation between the parameters induced by the non-linear behavior of the Debye-Einstein approach.In such a case, the Bayesian approach results in more realistic error estimations compared to classical error propagation analysis.In the case of many local minima, global regression analysis is recommended to reliably estimate the error of the regression, as is shown for an example from X-ray diffraction data analysis in [32].
The experimental data for C o p,m of SrMoO 4 , taken from [26], are plotted versus T in Figure 2. The solid line in Figure 2   The six parameters of the Debye-Einstein fit and their standard errors are calculated for both methods-the classical least squares (LSQ) method and the MCMC approach-and are listed in Table 1.In the case of the MCMC approach, the mean value of the probability density distribution and the highest probability is calculated, as well as the standard errors which follow from the standard deviations of all values calculated in the MCMC approach.The parameters obtained by both, completely different, regression approaches result in values for the parameters that are very close and even the estimated standard errors are similar.The probability density distributions of the model parameters are calculated, and presented in Figure 3.The probability density distribution for the Debye temperature T D is presented in Figure 3a, for the Einstein temperatures T E1 in Figure 3b and T E2 in Figure 3c, respectively.The probability density distributions of the prefactors m, n 1 and n 2 are shown in Figure 3d, Figure 3e and Figure 3f, respectively.
In addition, the parameters of the Debye-Einstein model function, as provided in Table 1, and the parameters of the function describing the errors of the heat capacities are simultaneously optimized.

Estimating the Uncertainties of Each Measurement Point
Our objective is to identify a function as simply as possible to approximate the standard errors of the data points investigated, where the experimental data are provided in [26].This function must obey the following two criteria:

•
The temperature dependency of the residuals should adequately describe the temperature dependency of the standard errors and vice versa, e.g., as the residuals increase with increasing heat capacities, the standard errors should also increase with increasing heat capacities.

•
The correlations between the parameters in the function for the standard errors should not be excessively high (e.g., above 90 percent), as this indicates the potential for using a simpler standard error function without significant data loss.
Since the distribution of the residuals is not known beforehand, the evaluation is carried out iteratively, and in case of failure, the entire analysis must be repeated by using another function for the standard errors.The following functions may be considered for describing the standard errors: It is worth noting that the function s must remain positive over the whole range of C o p,m .The simplest approach is to assign an equal, i.e., constant, standard error to all data points (Equation ( 13)).However, a better choice may consider the increase of the residuals with increasing heat capacities.A high correlation between s 0 and s 1 is observed when using a linearly increasing function (Equation ( 14)).When Equation ( 15) is used to describe the standard error function, a more realistic distribution of residuals is obtained.Figure 4 displays the residuals versus T. The residuals are calculated from the Markov chain containing the parameters with the highest probability.The function of the standard error s together with −s, i.e., the credible interval, is calculated from Equation ( 15) and plotted versus T in Figure 4.The parameters s 0 and s 1 of the function describing the standard error of C o p,m versus T are listed in Table 2.

Determining the Correlation between Parameters
In this section, the correlation between all parameters is determined using Equation (12).The resulting correlation matrix, which is symmetric, is presented in Table 3.The values specified in the correlation matrix can be visualized by scatter plots (Figure 5) showing the correlation of two selected parameters.The points in these scatter plots are color-coded.The color of the points changes with the frequency of hits (axis on the right) in the range of the parameter space represented by the point.
As a representative example, the correlation between the Debye temperature T D and the prefactor m is illustrated in Figure 5a, which is very high and at 0.98 close to 1.These two parameters are almost linearly related.However, neither of these two parameters can be omitted, since the Debye integral has to have a certain weight m not known before the regression analysis.The parameters T D and n 2 are slightly anti-correlated, with a value of −0.31, as can be seen in Figure 5b.Whereas the Einstein temperatures are strongly correlated at 0.83, as shown in Figure 5c, the prefactors n 1 and n 2 are practically not correlated and the value of 0.09 results in a scatter plot which is symmetric to the abscissa (Figure 5d).
In this example, the condition number of the correlation matrix r is calculated to be 5700, which indicates a high value.This high condition number suggests that the equation system is poorly conditioned.Therefore, from this perspective, the use of a more complex model (e.g., incorporating additional Einstein terms) is not recommended.
Based on the residuals analysis discussed in the section "Estimating the Uncertainties of Each Measurement Point", it can be concluded that underfitting is not observed in the examined dataset.Moreover, the inclusion of additional Einstein terms would not lead to a substantial reduction in the residuals.The question may arise if a four-parameter (4p) Debye-Einstein integral with a Debye temperature T D and an Einstein temperature T E and their prefactors m and n suffices to describe C o p,m (T) of SrMoO 4 from 0 K to 300 K.
Thus, the residuals following from a simpler four-parameter Debye-Einstein approach are calculated by the MCMC approach and presented in Figure 6.Compared to the residuals of the 6p-Debye-Einstein fit, shown in Figure 4, the 4p-fit results in almost five times larger residuals (Figure 6), which are not randomly distributed for this compound.This means that the 6p-Debye-Einstein fit seems to be the better option for describing the heat capacities C o p,m of SrMoO 4 than the simpler 4p-fit.In addition, it is shown in [12] that the 6p-Debye-Einstein approach leads to a heat capacity description with small standard errors of the fit parameters for many compounds.

Thermodynamic Functions
The molar entropy S o m (T) and molar enthalpy H o m (T) can be determined by using Equations ( 16) and (17), i.e., integrating the simulated molar heat capacities C o p,m (T) numerically.
The derived function −G o m /T, with G o m being the molar Gibbs energy, is obtained from: The values of the thermodynamic functions S o m , H o m and S o m (T) − H o m (T)/T of SrMoO 4 at T = 298.15K are presented in Table 4.The values are obtained from Levenberg-Marquardt least squares analysis (see also [12]).These values are compared to those of the highest probability calculated with the MCMC approach and to values from [26].This approach allows for the evaluation of entropy as a function of temperature for each set of fit parameters obtained from the Monte Carlo Markov chains.The individual entropy profiles serve as the basis for generating histograms at specific temperatures.
The histogram (probability density distribution) of the molar entropy S o m at T = 298.15K is shown in Figure 7 providing insights into the distribution of entropy values.The mean entropy value is determined to be 136.5196J mol −1 K −1 .Additionally, the entropy value with the highest probability corresponds to 136.5195 J mol −1 K −1 , representing the most likely entropy state of the system at the given temperature.
Furthermore, the standard deviation, a measure of the uncertainty of entropy values around the mean, can be calculated as ∆S o m = 0.033 J mol −1 K −1 .The probability density distribution of the enthalpy H o m at T = 298.15K is shown in Figure 8.
The probability density distribution of the function S o m − H o m /T at T = 298.15K is presented in Figure 9.It is worth mentioning that in some cases the highest probability lies not exactly at the position where the distribution evaluated by the "naked eye" expects the highest probability to be.This point can be explained as follows: The probability density distribution is obtained from all 4 • 10 5 Markov chain entries.A certain entry (six parameters of the Debye-Einstein fit and the two parameters describing the function of the standard errors of the heat capacities) has the highest value of (ln P), Equation (10).This maximum probability corresponds to the minimum obtained by the Levenberg-Marquardt approach, assuming that the same function is used for the standard errors.
In summary, determining thermodynamic functions within the Bayesian framework does not pose any difficulties.The Bayesian approach allows for the calculation of thermodynamic properties as a function of temperature based on the obtained sets of fit parameters.The resulting histograms provide insights into the distribution of values of the thermodynamic functions at specific temperatures.The mean value, along with the value corresponding to the highest probability, can be determined from the histograms.Additionally, the standard errors can be estimated from the probability density distributions of the thermodynamic functions at a specific temperature.

Conclusions
For the example of SrMoO 4 , it is again demonstrated in this work that the sixparameter Debye-Einstein fit for molar heat capacities C o p,m (T) works very well in the range 0-300 K, where SrMoO 4 could be replaced by many crystalline solids.Two different regression methods are applied for this task; the first is based on frequentist statistics using classical least squares, the second is an application of Bayes' theorem, numerically treated by the MCMC method.It is demonstrated that both completely different approaches not only lead to comparable results for the values of the parameters, but also to similar uncertainties of these parameters.
In addition, this investigation showcases the efficacy of the Bayesian framework to determine thermodynamic functions and their uncertainties.Based on the residuals, the parameters for the temperature-dependent function of the standard errors are optimized together with the model parameters and the Bayesian credible interval is obtained as the result.
From the correlation matrix of this example, it can be deduced that no more fitting parameters should be used in this temperature range as the correlation between these physically based parameters is partially very high.It can be seen as an advantage of the MCMC approach that the probability density distributions of the model parameters and of the derived quantities, such as the entropy S, enthalpy H and other thermodynamic functions, are revealed.Based on the results of this extented regression analysis of the molar heat capacities of SrMoO 4 , it can be recommended to use the 6p-Debye-Einstein integral approach as a standard method to fit heat capacities of many crystalline solids in the range between 0 K and 300 K.
H .The likelihood function evaluates how close the fit function containing the parameters ⃗ ξ approaches the experimental data.The product of the likelihood function P D ⃗ ξ, H and the prior probability P ⃗ ξ H is then normalized by the evidence P(D|H), resulting in the posterior probability P ⃗ ξ D, H .

Figure 1 .
corresponds to the least squares Levenberg-Marquardt fit of the six-parameter Debye-Einstein integral, computed by means of Origin2022b [33].The mean values of the C o p,m probability distribution densities are also plotted in Figure 2. Both, the classical least squares six-parameter (6p)-Debye-Einstein fit and the MCMC calculation mimic the experimental data almost perfectly well.Probability densities of the heat capacities C o p,m for SrMoO 4 at selected temperatures T. (a) Probability density of C o p,m at T = 15.0K. (b) Probability density of C o p,m at T = 98.1 K. (c) Probability density of C o p,m at T = 199.9K.(d) Probability density of C o p,m at T = 199.9K.

Figure 2 .
Figure 2. Experimental C o p,m values from [26] versus T approximated by the 6p-Debye-Einstein integral fit; the mean values of the C o p,m values obtained by the MCMC-method are also shown.

Figure 3 .
Probability densities of the 6 parameters of the Debye-Einstein integral.(a) Probability density of the Debye temperature T D .(b) Probability density of the Einstein temperature T E1 .(c) Probability density of the Einstein temperature T E2 .(d) Probability density of the factor m. (e) Probability density of the factor n 1 .(f) Probability density of the factor n 2 .

Figure 4 .
Figure 4. Residuals (difference between simulated and measured heat capacities, C p,sim − C p,meas ) as a function of temperature T. The function s of the standard errors and this function mirrored at the abscissa, i.e., −s, are also plotted versus T.

Figure 5 .
Color -coded scatter plots visualizing the correlation between certain parameters.(a) Correlation between T D and m.(b) Correlation between T D and n 2 .(c) Correlation between T E1 and T E2 .(d) Correlation between n 1 and n 2 .

Figure 6 .
Figure 6.Residuals (difference between simulated and measured heat capacities, C p,sim − C p,meas ) in the case of a simpler 4p-Debye-Einstein fit as a function of temperature T. The function s of the standard errors and this function mirrored at the abscissa, i.e., −s are also plotted versus T.

Figure 7 .
Figure 7. Probability density distribution of S o m at T = 298.15K.

Figure 8 .
Figure 8. Probability density distribution of H o m at T = 298.15K.

Figure 9 .
Figure 9. Probability density distribution of S o m − H o m /T at T = 298.15K.

Table 1 .
Six-parameter Debye-Einstein fit by classical LSQ and by MCMC.

Table 2 .
The function for the standard error of the heat capacity.

Table 4 .
Thermodynamic functions of SrMoO 4 at T = 298.15K derived from molar heat capacity.