Efficient Debye Function Interpolation Formulae: Sample Applications to Diamond

Roland Pässler ^*

Technische Universität Chemnitz, Institut für Physik, D-09107 Chemnitz, Germany

* Correspondence: Roland Pässler

Academic Editor: Hossein Hosseinkhani

Special Issue: Quantum Mechanics in Solid State Systems

Received: June 17, 2021 | Accepted: September 22, 2021 | Published: November 22, 2021

Recent Progress in Materials 2021, Volume 3, Issue 4, doi:10.21926/rpm.2104042

Recommended citation: Pässler R. Efficient Debye Function Interpolation Formulae: Sample Applications to Diamond. Recent Progress in Materials 2021; 3(4): 042; doi:10.21926/rpm.2104042.

© 2021 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

A basic model for approximate numerical simulations of measured heat capacities had been constructed more than 100 years ago by Debye [1]. This familiar model was based on the assumption that the relevant phonon-densities-of-states (PDOS) spectra of metals as well as non-metals, should be describable to a good approximation in the form of purely quadratic phonon energy dependences, ${g_P^{ }}\left( \varepsilon \right) \propto {\varepsilon ^2}$, from 0 up to a certain cut-off (Debye) energy, ${\varepsilon _D}$. This primary heat theory for solids [1] was focused, above all, on derivations of usable analytical formulae for numerical calculations of isochoric heat capacities of crystal lattices (${C_V}\left( T \right)$), which are generally defined as being given by the partial derivatives, ${C_V}\left( T \right) \equiv {(\partial U\left( {T,V} \right)/\partial T)_V}$, of internal energies, $U\left( {T,V} \right)$, with respect to the lattice temperature, $T$ (at a fixed volume $V$). According to this theory [1], it should have been possible to describe the temperature dependence of ${C_V}\left( T \right)$ in terms of certain Debye function integrals [1], the $T$-dependent magnitudes of which are governed by a unique model-specific parameter, the so-called Debye temperature, ${\Theta _D} = {\varepsilon _D}/{k_B}$ (where ${k_B}$ represents the Boltzmann constant).

However, since the middle of the past century, a wealth of experimental studies on the thermal properties of non-metals (semiconductors as well as isolators) has been performed. From the numerous results of these studies, one could conclude that, as a rule, it becomes possible to simulate the measured heat capacity data using the Debye function integrals only under the assumption that the Debye temperature depends (more or less strongly) on the lattice temperature, $T$. The corresponding ${\Theta _D}\left( T \right)$ dependencies have been found to be very pronounced for typical semiconductor materials, in particular, for Si and Ge [2,3] as well as for numerous III-V materials [4,5,6] and II-VI materials [7,8,9,10,11].

These material-specific ${\Theta _D}\left( T \right)$ dependences show “snaky” shapes, which exhibit the following:

(i) Rapid fall from their $T \to 0$ limiting magnitudes [12], ${\Theta _D}\left( 0 \right)$, to significantly lower minima, ${\Theta _{D\;min}}$, which are usually located at temperatures of the order ${\Theta _D}\left( 0 \right)/12$

(ii) Subsequent rise up to a certain maxima, ${\Theta _{D\;max}}$, which are frequently located at temperatures of the order ${\Theta _D}\left( 0 \right)/2$

(iii) Final drop to zero, ${\Theta _D}\left( {{T_o}} \right) = \;0$, at the characteristic points on the $T$-scale, ${T_o}$, where the measured ${C_p}\left( T \right)$ curves are crossing the classical (limiting) Dulong-Petit value for harmonic lattice heat capacities, ${C_p}\left( {{T_o}} \right) = {C_{Vh}}\left( \infty \right) = \;3{n_A}R$ (where $R$ represents the gas constant and ${n_A}$ is the number of atoms per molecule of the material in consideration).

Thus, in the case of these typical semiconductor materials, as well as various alkali halides [13,14,15,16], the maxima of these “snaky” ${\Theta _D}\left( T \right)$ curves are significantly higher than the associated minima, ${\Theta _{D\;max}} > {\Theta _{D\;min}}$. Consequently, sufficiently broad temperature intervals do not usually exist within the confines of which the Debye temperature would indeed be approximately constant. This excludes the possibility of describing, at least, some limited sections of the given isochoric or isobaric heat capacity $T$-dependences, ${C_{V/p}}\left( T \right)$, in terms of the Debye function integrals governed by fixed, material-specific Debye temperatures, ${\Theta _D}$.

On the other hand, among a larger variety of wide-band-gap materials (${E_g}$ > 3 eV), one can find at least a few cases where the differences between ${\Theta _{D\;max}}$ and ${\Theta _{D\;min}}$ are relatively small (being limited to a few %). This indicates that the ${\Theta _D}\left( T \right)$ values within the respective $T$-intervals, from approximately ${\Theta _D}\left( 0 \right)/12$ up to approximately ${\Theta _D}\left( 0 \right)/2$, are more or less constant. Such a limited quasi-plateau behavior, ${\Theta _D}\left( T \right) \approx {\Theta _D} = const.$, of a ${\Theta _D}\left( T \right)$ curve above 200 K, was well known already long ago for diamond (see Figure 1 in earlier studies [17], Figure 3 in [18], and cf. Figure 2 in Section 4 of the present paper). A similar plateau behavior of the ${\Theta _D}\left( T \right)$ curve has also been observed for cubic carbon nitride (c-BN) above 150 K (see Figure 2 in [19]). These two examples of wide-gap materials might already be sufficient to justify the present aim of deriving efficient formulae for the conventional Debye function integral (governed by fixed ${\Theta _D}$), which should hence enable us to perform careful least-mean-square fits of the duly limited (partial) ${C_{V/p}}(T > {\Theta _D}\left( 0 \right)/12)$ datasets (similar to those available for diamond or c-BN) without involving numerical calculation procedures for the original Debye function integrals.

In Section 2, we give a brief overview of the derivation of the conventional low-and high-temperature approximations [11] for the Debye function integral, which are well known already from Debye’s original paper [1]. We confirm that Debye’s low-temperature formula [1,11,20,21] is in principle capable of providing quite adequate numerical values for the Debye function integral even at lattice temperatures comparable with the Debye temperature, provided the respective summation procedure is not truncated at a very low order of summation terms. In contrast to the latter, we find again (in accordance with [1,3,11,22]) that Debye’s original high-temperature Taylor series expansion is of little practical use in view of its rather bad convergence properties, which are caused by the alteration of signs of the respective expansion coefficients, i. e. $c_{2m}^D = {( - 1)^m}\left| {c_{2m}^D} \right|$ for the subsequent even-ordered ${\left( {\frac{{{\Theta _D}}}{T}} \right)^{2n}}$ terms ($n$ = 1, 2, 3, …).

Significant improvements of convergence properties have already been shown in previous studies [11,22] to be realizable by means of suitable transformations [23] of the given Taylor series expansions into structurally different versions that are embedded into conveniently chosen alternative (non-linear) functions. In the appendix, we briefly present two previous sample applications of this method to the high-temperature behavior of the Debye function, which resulted in the construction of an exponential series representation [22], and a derivation of an alternative Taylor series representation for reciprocal Debye function values, ${\left( {{\kappa _{Dh}}\left( x \right)} \right)^{ - 1}}$ [11].

Of special usefulness within the present context is an unprecedented Taylor series representation for the squares of the reciprocal Debye function values, ${\left( {{\kappa _{Dh}}\left( x \right)} \right)^{ - 2}}$, which has been presented in the appendix. In Section 3, we qualitatively consider different combinations of truncated versions of the novel high-capacity Debye function formula, (Eq. (A.12)), with various truncated versions of Debye’s conventional low-capacity Debye function expansion, ${\kappa }_{Dl}^{\left(n_D^{}\right)}\left(x\right)$(Eq. (2.2)).

In Section 4, we visualize the actual temperature dependence of the effective (caloric) Debye temperature for diamond, ${\Theta _D}\left( T \right)$, which is implied by the isobaric heat capacities, ${C_p}\left( T \right)$, given by Desnoyers and Morrison [17], Victor [24], and Dinsdale [25]. In addition, a least-mean-square fitting has been performed for the whole combination of the respective ${C_p}\left( T \right)$ datasets, from approximately 200 K up to 5000 K. This fit is based on an appropriate low-to-high temperature interpolation formula in combination with additional anharmonicity-related terms, which come strongly into play in the high-temperature region. As the main result of this fit, we found the effective Debye temperature to be ${\Theta _D} \cong \;$1855 K.

In Section 5, we give a brief sketch of the basic equations of our previously developed representative hybrid model [11,26], which consists of two continuous low-$T$ curve sections in combination with three discrete (Einstein) phonon energy peaks. We have also performed the corresponding high-accuracy fit of the whole ${C_p}\left( T \right)$ dataset, from approximately 20 K up to 5000 K.

The two theoretical alternative curves for the ${C_p}\left( T \right)$ dependence of diamond, resulting from the fit via the unprecedented high-to-low-$T$ interpolation formulae (in Section 4), on the one hand, and from the high-accuracy fit via the previously developed hybrid model [11,26] (in Section 5), on the other hand, have been found to be almost indistinguishable throughout the interval from 200 K to 5000 K. This good agreement illustrates the potential usefulness of the computational framework of approximate Debye function formulae (displayed in Sec. 3) developed in this work and its possible application to the limited low-to-high-temperature heat capacity datasets, which are available for certain wide-band-gap materials characterized by sufficiently low degrees of phonon dispersion.

Various other aspects and additional quantitative information resulting from the comparable alternative fittings of the ${C_p}\left( T \right)$ datasets, for diamond, under study, by using alternative models have been discussed in Section 6.

2. Conventional Low-and High-temperature Expressions for the Debye Function

According to Debye’s original heat capacity model [1], the isochoric lattice heat capacity is assumed to be given by the product, ${C_V}\left( T \right) = {C_{Vh}}\left( \infty \right){\kappa _D}\left( T \right)$ of the familiar Dulong-Petit limit, ${C_{Vh}}\left( \infty \right) = \;3{n_A}R$, of the harmonic lattice heat capacity (for a material containing ${n_A}$ atoms per molecule) and a normalized (to unity) heat capacity shape function, $0 < {\kappa _D} \le 1$. The latter is well known to be defined by an integral of the form [1,4,5,10,11,12,21]

\[ \kappa_{D}(x) \equiv \frac{3}{x^{3}} \int_{0}^{x} d z \frac{z^{4} e^{z}}{\left(e^{z}-1\right)^{2}}=\frac{12}{x^{3}} \int_{0}^{x} d z \frac{z^{3}}{e^{z}-1}-\frac{3 x}{\exp (x)-1} \tag{2.1} \]

where the ratio between the (presumably material-specific) Debye temperature, ${\Theta _D}$, and the lattice temperature, $T$, is denoted as usual by $x = {\Theta _D}/T$.

Yet, with respect to the eventual numerical applications of this conventional model to the interpretation of certain experimental data (with respect to which this simple model happens to be applicable), it is desirable to avoid the numerical integrations (in Eq. (2.1)) by using the appropriate analytical or algebraic expressions for more or less extended regions of low and/or high temperatures.

For low-temperature regions ($x = {\Theta _D}/T > > 1$), it was found [1] to be convenient to expand the denominator in the second integral of Eq. (2.1) into an infinite Taylor series, ${\left( {{e^z} - 1} \right)^{ - 1}} = {e^{ - z}}{\left( {1 - {e^{ - z}}} \right)^{ - 1}} = {e^{ - z}} + {e^{ - 2z}} + {e^{ - 3z}} + \ldots $. By performing the corresponding integrations (subsequently, by parts), within the second integral of Eq. (2.1), for the low-temperature behavior of the Debye function, ${\kappa _D}\left( x \right)$, an infinite series expression of the familiar form [1,11,12,20,21] is obtained:

$${\kappa }_{Dl}^{\left(n_D^{ }\right)}(x)=\frac{4{\pi }^4}{5x^3}-\frac{3x}{\exp (x)-1}-12x\sum _{n=1}^{n_D^{}(\to \mathrm{\infty })}\exp (-nx)[\frac{1}{nx}+\frac{3}{(nx{)}^2}+\frac{6}{(nx{)}^3}+\frac{6}{(nx{)}^4}].\mathrm{}(2.2)$$

Concerning the actual range of validity of this basic Debye function series expansion (2.2), we would like to point out that it is, in principle, quite exact within the frame of a duly extended summation procedure (i. e., for the ${n_D^{ }} \to \infty \;$limit, at least). It may thus be used, if necessary, for performing high-accuracy calculations of the ${\kappa _D}\left( x \right)$ function (Eq. (2.1)) even for regions of relatively high temperatures ($x = {\Theta _D}/T \approx 1$; see below) provided that the summation procedure is extended to sufficiently high order (e. g., up to $n \approx 30$). At the same time, it is obvious that this infinite series representation (Eq. (2.2)) is not applicable to the $T \to \infty $ limit.

Concerning the opposite regime of the high-temperature behavior of the ${\kappa _D}\left( x \right)$ integral, it is useful to rewrite Eq. (2.1), in accordance with $\left( {{e^{z/2}} - {e^{ - z/2}}} \right)/2\; = \sinh (z/2)$, in the equivalent form [3,10,11,12]

\[ \kappa_{D}(x)=\frac{3}{4 x^{3}} \int_{0}^{x} d z \frac{z^{4}}{\left(\sinh \left(\frac{z}{2}\right)\right)^{2}} \leq 1 \tag{2.3} \]

For the $\left( {\sinh (z/2} \right){)^{ - 2}}$ function one can use, in analogy to previous works [3,11,22], the Taylor series expansion

\[ \left(\sinh \left(\frac{z}{2}\right)\right)^{-2}=4 z^{-2}+4 \sum_{m=1}^{\infty} \frac{(-1)^{m}(2 m-1)\left|B_{2 m}\right| z^{2 m-2}}{(2 m) !} \tag{2.4} \]

where $\left| {{B_{2m}}} \right|$ denotes the absolute values of the respective Bernoulli numbers (i. e., $\left| {{B_2}} \right|$= 1/6, $\left| {{B_4}} \right|$= 1/30, $\left| {{B_6}} \right|$= 1/42, $\left| {{B_8}} \right|$= 1/30, $\left| {{B_{10}}} \right|$= 5/66, $\left| {{B_{12}}} \right|$ = 691/2730, $\left| {{B_{14}}} \right|$= 7/6, $\left| {{B_{16}}} \right|$= 3617/510,…). Using the expansion in Eq. (2.4) in Eq. (2.3) and performing the respective integrations of the individual ${z^{2m + 2}}$ power terms, for the high-temperature behavior of the ${\kappa _D}\left( x \right)$ function, one obtains a Taylor series expansion of the conventional form [1,11,20,21,22]

\[ \kappa_{D h}^{\left(m_{D}^{ }\right)}(x)=1+\sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} c_{2 m}^{D} x^{2 m} \tag{2.5} \]

where the expansion coefficients are given (in accordance with previous studies [1,11,20]) by

\[ c_{2 m}^{D}=(-1)^{m} \frac{3 \cdot(2 m-1)\left|B_{2 m}\right|}{(2 m+3) \cdot(2 m) !} \tag{2.6} \]

i. e. explicitly [1,11,20,21,22]

\[ c_{2}^{D}=\frac{-1}{20}, c_{4}^{D}=\frac{+1}{560}, c_{6}^{D}=\frac{-1}{18144}, c_{8}^{D}=\frac{+1}{633600}, c_{10}^{D}=\frac{-1}{23063040} \tag{2.7} \]

(Note that the values of Debye’s conventional expansion coefficients, $c_{2m}^D$ (2.6), are listed up to an order of $2m \le 16$ in the left part of Table 1).

However, due to the alternating signs of Debye’s original Taylor series expansion coefficients (2.7), a high-temperature series expansion of this conventional type (Eq. (2.5)) does not show good convergence properties (as already pointed out in Debye’s original paper [1] and illustrated in detail in previous studies [11,22]). Thus, one has to look for possible alternative forms of high-temperature representations of the Debye integral.

3. Alternative High-temperature Representations for the Debye Function

A very useful analytical tool for constructing alternative versions of the Taylor series representations for the high-temperature behavior, ${\kappa _{Dh}}\left( x \right)$ (Eq. (2.5)), of the Debye function integral, ${\kappa _D}\left( x \right)$ (Eq. (2.1)), has been found to be provided by the method of performing suitable transformations [23] of a given Taylor series expansion into structurally different counterparts that are embedded into conveniently chosen alternative (non-linear) functions. This method had been used for the first time in [22] for deriving an equivalent Taylor series for the logarithm, $\log ({\kappa _{Dh}}\left( x \right))$ (Eq. (A.2)), of the conventional ${\kappa _{Dh}}\left( x \right)$ representation (Eq. (2.5)). Inverting this relationship and observing that $\exp (\log ({\kappa _{Dh}}\left( x \right))) = {\kappa _{Dh}}\left( x \right)$, in [22], we came to a corresponding exponential series representation (Eq. (A.5)) for the ${\kappa _{Dh}}\left( x \right)$ function which showed markedly better convergence properties (see Figure 6 in [22]) than Debye’s original Taylor series representation (Eq. (2.5); for additional details, see [22] and subsection A.1 of the appendix of the present paper).

In an analogous manner, we have derived an equivalent Taylor series expansion in [11] for the reciprocal Debye function value, $1/{\kappa _{Dh}}\left( x \right)$ (Eq. (A.7/8)). Inverting the latter relationship, so that ${\left[ {1/{\kappa _{Dh}}\left( x \right)} \right]^{ - 1}} = {\kappa _{Dh}}\left( x \right)$, we came to an alternative analytical expression (Eq. (A.9)) for the ${\kappa _{Dh}}\left( x \right)$ function which was found to show still significantly better convergence properties than the preceding exponential representation (for more details, see [11] and subsection A.2 of the appendix).

Table 1 The eight lowest-order expansion coefficients, $r_{2 \le 2m \le 16}^D$, implied by the unprecedented high-capacity (high-temperature) algebraic formula, (Eq. (A.12)), in comparison with their conventional counterparts, $c_{2 \le 2m \le 16}^D$, due to Debye’s original version (see Eq. (12’) in [1]) of the high-temperature Taylor series expansion, (Eq. (25)).

However, with respect to the present purpose of simulating experimentally measured heat capacities of a low-dispersion wide-band-gap material such as diamond (see Section 4), it is useful to display still another version of a rapidly converging series expansion for the ${\kappa _{Dh}}\left( x \right)$ function. The corresponding investigation has been done below (in Subsection A.3 of the Appendix) by considering the squares of the reciprocal Debye function values ${(1/{\kappa _{Dh}}\left( x \right))^2}$ (Eq. (A.10)). Inverting the latter relationship, so that ${\left[ {{{(1/{\kappa _{Dh}}\left( x \right))}^2}} \right]^{ - \frac{1}{2}}} = {\kappa _{Dh}}\left( x \right)$, In the present study (in Appendix A.3), we came to a novel (unprecedented) expression (Eq. (A.12)) for the ${\kappa _{Dh}}\left( x \right)$ function having the algebraic form,

\[ \kappa_{D h}^{\left(m_{P}^{ }\right)}(x)=\left[1+\sum_{m=1}^{m_{P}^{ }(\rightarrow \infty)} r_{2 m}^{D} x^{2 m}\right]^{-\frac{1}{2}} \tag{3.1} \]

This form has been found to be extraordinarily well suited for the present purpose. The actual values of the respective Taylor series expansion coefficients, $r_{2m \le 16}^D$, are listed up to an order of $2m \le 16$ in the right column of Table 1.

Figure 1 Visualization of the rapid approach of a series of truncated versions of the high-temperature algebraic representation, developed in the present work (-----curves, obtained using Eq. (3.1)), for ${m_P^{ }}$ = 1 to 4), toward the exact ${\kappa _D}\left( x \right)$ function (—— curve, obtained using Eq. (2.1)). The upper inset shows that the increasing deviations of the truncated curve from the exact ${\kappa _D}\left( x \right)$ curve remain limited to less than 1% throughout an unusually large high-temperature interval (of $0 \le x < 8$). The “snaky” solid and dotted curves show that the actual deviations of the approximate ${\tilde \kappa _D}\left( x \right)$ dependences due to the relatively simple interpolation formulae given by Eqs. (3.5) and (3.6), are limited to ±0.5% or ±0.25%, respectively. From the middle inset, we see, e.g., that a significant reduction of deviations to only about ±0.007% or ±0.0006% can be achieved by using the interpolation formulae of the type ${\tilde{\kappa }}_D^{\left(n_D^{}\ge 1;m_P^{}=4\right)}\left(x\right)$ (Eq. (3.7)) involving truncated versions of Debye’s conventional low-capacity representation, (Eq. (2.2)), of the two lowest orders, ${n_D^{ }}$= 1 or 2, respectively.

3.1 Alternate Usage of High-and Low-temperature Debye Function Formulae

From Figure 1, we see (in analogy to Figure 8 of an earlier study [11]) that the deviations, ${\kappa }_{Dl}^{\left(n_D^{ }\ge 1\right)}\left(x\right)-{\kappa }_D\left(x\right)>0$, of Debye’s conventional low-temperature curves (Eq. (2.2)) are decreasing with the increasing number of summation steps, ${n_D^{ }}$= 1, 2, …, and they are increasing with decreasing magnitude of the argument $x\left( { = {\Theta _D}/T} \right)$. In contrast to the latter, the deviations, ${\kappa }_{Dh}^{\left(m_P^{ }=7\right)}\left(x\right)-{\kappa }_D\left(x\right)>0$, of the high-temperature curve (Eq. (3.1)) constructed in this study, is seen to monotonically increase with increasing magnitude of the argument $x\left( { = {\Theta _D}/T} \right)$. Consequently, the individual pairs of these two complementary low-and high-temperature curves are mutually crossing, ${\kappa }_{Dl}^{\left(n_D^{ }\ge 1\right)}\left(x_c\right)={\kappa }_{Dh}^{\left(m_P^{ }=7\right)}\left(x_c\right)$, at certain (${n_D^{ }}$-dependent) crossing points, ${x_c}$. We see from a corresponding list for a series of ${n_D^{ }}$ values (in Table 2) that the individual crossing point positions, ${x_c}$, and the corresponding deviations from the exact ${\kappa _D}\left( {{x_c}} \right)$ values,

$\begin{array}{c}{\kappa }_{Dl}^{\left(n_D^{}\ge 1\right)}\left(x_c\right)-{\kappa }_D\left(x_c\right)={\kappa }_{Dl}^{\left(m_P^{}= 7\right)}\left(x_c\right)-{\kappa }_D\left(x_c\right)>0,\end{array}$

$(3.2)$

are rapidly decreasing with increasing ${n_D^{ }}$ (cf. the positions of the solid circles, for the cases of ${n_D^{ }}$ = 1 and 2, in the middle inset to Figure 2).

Table 2 Positions of the crossings, ${x_c}$, of the unprecedented high-temperature curve, ${\kappa }_{Dh}^{\left(m_P^{ }= 7\right)}(x<7)$ (Eq. (3.1);----curve shown in the middle inset of Figure 1), with various (truncated) samples of Debye’s conventional low-capacity (low-temperature) curves, ${\kappa }_{Dl}^{\left(n_D^{ }= 1,2,\mathrm{...}\right)}\left(x\right)$ (Eq. (2.2); ·-·-·-· curves shown, for ${n_D^{ }}$ = 1 and 2, in the middle inset of Figure 1). In the third column, we have quoted the maximum values of the relative deviations, $\delta {\hat \kappa _D}\left( {{x_c}} \right)/{\kappa _D}\left( {{x_c}} \right)$, between the exact ${\kappa _D}\left( x \right)$ curve (Eq. (2.1)) and the compound ${\widehat{\kappa }}_D^{\left(n_D^{ }\ge 1;m_P^{ }=7\right)}\left(x\right)$ curves (Eq. (3.4)).

In this way, we come to rather good approximations, ${\widehat{\kappa }}_D^{\left(n_D^{ } ; m_P^{ }=7\right)}\left(x\right)$, for the Debye integral in terms of ${n_D^{ }}$-specific combinations of complementary high-and low-temperature curve sections,

\[ \hat{\kappa}_{D}^{\left(n_{D}^{ } \ ; \ m_{P}^{ }=7\right)}(x)=\kappa_{D h}^{\left(m_{P}^{ }=7\right)}(x),\ \text {for} \ x \leq x_{c}\left(n_{D}^{ }\right),\ \text {and} \ \hat{\kappa}_{D}^{\left(n_{D}^{ } \ ; \ m_{P}^{ }=7\right)}(x)=\kappa_{D l}^{\left(n_{D}^{ }\right)}(x), \ \text { for } \ x \geq x_{c}\left(n_{D}^{ }\right).\tag{3.3} \]

Equivalently, we can also express the latter combinations in the more compact form,

\[ \hat{\kappa}_{D}^{\left(n_{D}^{ } ; m_{P}^{ }=7\right)}(x)=\min \left\{\kappa_{D h}^{\left(m_{P}^{ }=7\right)}(x) ; \kappa_{D l}^{\left(n_{D}^{ }\right)}(x)\right\} \text { for } 0 \leq x<7 \tag{3.4} \]

(and ${\widehat{\kappa }}_D^{\left(n_D^{ };m_P^{ }=7\right)}\left(x\right)={\kappa }_{Dl}^{\left(m_D^{ }\right)}\left(x\right)$, for $x > 7$). In this way, owing to the requirement (Eq. (3.4)) of choosing just the lower value of the two competitive (approximate) high-and low-temperature functions, ${\kappa }_{Dh}^{\left(m_P^{ }= 7\right)}\left(x\right)$ versus ${\kappa }_{Dl}^{\left(n_D^{ }\right)}\left(x\right)$, it is automatically assured that the switch between the two complementary curve sections occurs just at the respective crossing point, $x_c\left(n_D^{ }\right)$ .

From Table 2, we see that by limiting the summation procedure for ${\kappa }_{Dl}^{\left(n_D^{ }\right)}\left(x\right)$ (in Eq. (2.2)) e.g. to only 4 steps (${n_D^{ }}$= 4), one nevertheless comes to a rather gentle approximation for the Debye function, ${\kappa }_D\left(x\right)\cong {\kappa }_D^{\left(n_D^{ }= 4;m_P^{ }=7\right)}\left(x\right)$, the maximum deviations of which from the exact function are limited to an order of 10⁻⁷. Such a high degree of accuracy is largely sufficient with respect to practical applications.

3.2 Interpolation Formulae

It turns out that it is possible to find yet another way of constructing good analytical approximations for the Debye function integral, ${\kappa _D}\left( x \right)$ (Eq. (2.1)), which does not involve a sudden switch between the high-and low-temperature curve sections at their crossing points. This alternative way is based on the construction of suitable interpolation formulae that bridge the gap between a truncated version of Debye’s low-temperature expression (Eq. (2.2)), ${\kappa }_{Dl}^{\left(n_D^{}\right)}\left(x\right)>{\kappa }_D\left(x\right)$, and a suitably chosen sample of the novel analytical high-temperature expressions, ${\kappa }_{Dh}^{\left(m_P^{}\right)}\left(x\right)$ (Eq. (3.1)), the curve of which should throughout be running below the exact ${\kappa _D}\left( x \right)$ curve. We see from the upper and middle insets of Figure 1 that this feature is found to be realized in particular for the ${m_P^{ }} = \;4$ case.

Let us first consider the simplest case of an interpolation formula, a combination of the high-temperature expression, ${\kappa }_{Dh}^{\left(m_P^{}= 4\right)}\left(x\right)$ (Eq. (3.1), with the asymptotic ($x \to \infty $ limiting) term of the conventional low-temperature expression (Eq. (2.2)), i.e., $\kappa _{Dl}^{(\lim .)}\left( x \right) = \;4{\pi ^4}/5{x^3}$ [1,11,12,20]. An appropriate interpolation formula for connecting this pair of non-crossing high-and low-temperature curves can be chosen to be of the form

\[ \tilde{\kappa}_{D}(x)=\sqrt[2]{\frac{\exp \left(\frac{x_{o}-x}{\Delta x}\right)+1}{\left(1+\sum_{m=1}^{m_{P}^{ }=4} r_{2 m}^{D} x^{2 m}\right) \cdot \exp \left(\frac{x_{o}-x}{\Delta x}\right)+\left(\frac{5}{4 \pi^{4}}\right)^{2} x^{6}}}\tag{3.5} \]

where ${x_o}$ and $\Delta x$ are formula-specific interpolation parameters.

An optimal approach of the interpolation curve, ${\tilde \kappa _D}\left( x \right)$ (Eq. (3.5)), to the exact ${\kappa _D}\left( x \right)$ dependence could be achieved by performing a corresponding least-mean-square fitting, which leads to the adjusted values of ${x_o}$= 9.775 and $\Delta x$ = 0.84 for the respective interpolation coefficients. The remaining deviations $\left( {{{\tilde \kappa }_D}\left( x \right) - {\kappa _D}\left( x \right)} \right)/{\kappa _D}\left( x \right)$, of Eq. (3.5) from the exact Debye function, ${\kappa _D}\left( x \right)$, are visualized by the solid (“snaky”) curve in the upper inset of Figure 1. From the latter, we see that these deviations from the exact Debye function (Eq. (2.1)), which are implied by this primary interpolation formula (Eq. 3.5)), are limited to an order of ±0.5%.

A moderate reduction of deviations from the exact Debye function (Eq. (2.1)) can be readily achieved by including the subsequent term, namely, the term $ - 3x/\left( {{e^x} - 1} \right)$) occurring in Debye’s low-temperature expression, ${\kappa _{Dl}}\left( x \right)$ (Eq. (2.2)), into a conveniently constructed interpolation formula. Accordingly, we consider a combination of a truncated low-temperature expression of the type ${\kappa _{Dl}}\left( x \right) \to 4{\pi ^4}/5{x^3} - 3x/\left( {{e^x} - 1} \right)$ with the same high-temperature expression as above, ${\kappa }_{Dh}^{\left(m_P^{}= 4\right)}\left(x\right)$ (Eq. (3.1). We have found an appropriate version of an interpolation formula for this combination to be given by the expression

\[ \tilde{\kappa}_{D}^{(0)}(x)=\frac{\left[\left(1+\sum_{m=1}^{m_{P}^{ }=4} r_{2 m}^{D} x^{2 m}\right)^{-\frac{1}{2}} \cdot \exp \left(\frac{x_{o}-x}{\Delta x}\right)+\frac{4 \pi^{4}}{5 x^{3}}-\frac{3 x}{e^{x}-1}\right]}{\left[\exp \left(\frac{x_{o}-x}{\Delta x}\right)+1\right]} \tag{3.6} \]

involving adjusted values of ${x_o}$= 9.05 and $\Delta x$ = 0.60 for the respective interpolation coefficients. The remaining deviations $\left( {\tilde \kappa _D^{\left( 0 \right)}\left( x \right) - {\kappa _D}\left( x \right)} \right)/{\kappa _D}\left( x \right)$, of Eq. (3.6) from the exact Debye function, ${\kappa _D}\left( x \right)$, are visualized by the dotted (“snaky”) curve in the upper inset in Figure 1. From the latter, we can see that the respective deviations of Eq. (3.6) from the exact Debye function (Eq. (2.1)), are limited to an order of ±0.25%. Thus, alternative use of Eq. (3.6), instead of Eq. (3.5), involves a reduction of deviations from the exact Debye function by a factor of about 2.

For the sake of completeness, we would like to mention that an even more pronounced (order-of-magnitude) reduction of deviations from the exact ${\kappa _D}\left( x \right)$ function, which is inherent to interpolation formulas of the type given by Eq. (3.6), can be achieved when one envisages combinations of the same high-temperature expression, ${\kappa }_{Dh}^{\left(m_P^{}= 4\right)}\left(x\right)$ (Eq. (2.1)), with low-temperature expressions of Debye’s type, ${\kappa }_{Dl}^{\left(n_D^{}\ge 1\right)}\left(x\right)$ (Eq. (2.2)), that is retaining at least one summation step, i.e.,

\[ \tilde{\kappa}_{D}^{\left(n_{D}^{ } \geq 1 ; m_{P}^{ }=4\right)}(x)=\frac{\left[\left(1+\sum_{m=1}^{m_{p}^{ }=4} r_{2 m}^{D} x^{2 m}\right)^{-\frac{1}{2}} \cdot \exp \left(\frac{x_{o}-x}{\Delta x}\right)+\kappa_{D l}^{\left(n_{D}^{ } \geq 1\right)}(x)\right]}{\left[\exp \left(\frac{x_{o}-x}{\Delta x}\right)+1\right]} \tag{3.7} \]

In this connection, we succeeded in realizing rather good simulations of the true ${\kappa _D}\left( x \right)$ function by an interpolation formula of the latter type (Eq. (3.7)) even for the two lowest-order versions, where Debye’s low-temperature expansion, ${\kappa }_{Dl}^{\left(n_D^{}\right)}\left(x\right)$ (Eq. (2.2)), has been truncated just after the first or second summation step (i.e., $n_D^{}$= 1 or 2).

Appropriate interpolation coefficients for $n_D^{}= 1$ were found to be ${x_o}$= 5.393 in combination with $\Delta x$ = 0.24, involving a limitation of deviations to less than ±0.007% (cf. the solid (“snaky”) curve in the middle inset in Figure 1). For $n_D^{}= 2$, the adjusted interpolation coefficients turned out to be ${x_o}$= 4.146 in combination with $\Delta x$ = 0.15, involving an order-of-magnitude reduction of deviations to less than ±0.0006% (cf. the corresponding solid (“snaky”) curve in the middle inset in Figure 1).

Thus, if necessary, one can also use an interpolation formula of the type given by Eq. (3.7) in order to come to a high-accuracy approximation for the Debye function integral provided that at least two summation steps are retained in the corresponding low-capacity expression, ${\kappa }_{Dl}^{\left(n_D^{}\ge 2\right)}\left(x\right)$ (Eq. (2.2)). An advantageous feature of any one of the interpolation formulae displayed above (i.e., Eq. (3.5) to (3.7)) is due to the fact that the slopes (as well as the higher-order derivatives) of the approximate ${\tilde \kappa _D}\left( x \right)$ curves are continuous functions of their argument, $x = {\Theta _D}/T$. At variance to this, the slopes (including the higher-order derivatives) of approximate functions of the type ${\widehat{\kappa }}_D^{\left(n_D^{}\ge 1;m_P^{}=7\right)}\left(x\right)$ (Eq. (3.4)) are changing discontinuously at the respective points of intersection, ${x_c}$, between the connected high-and low-temperature curve sections, ${\kappa }_{Dh}^{\left(m_P^{}= 7\right)}\left(x\right)$ and ${\kappa }_{Dl}^{\left(n_D^{}\right)}\left(x\right)$.

4. Sample Applications to Diamond

It is a matter of principle that any one of the preceding approximate expressions (displayed in the preceding sections 2 and 3) for the Debye function integral, ${\kappa _D}\left( x \right) = {\kappa _D}\left( {{\Theta _D}/T} \right)$ (Eq. (2.1)), are applicable only to those heat capacity datasets for, which the basic assumption of a constant Debye temperature, ${\Theta _D}\left( T \right) \approx {\Theta _D}$ = constant, is realized to a good approximation (at least within a limited low-to-high temperature region). Such a special state of affairs can hardly be expected to be found for typical semiconductor materials (${E_g}$ < 3 eV), and it is only very rarely encountered even for wide-band-gap materials (${E_g}$ > 3 eV). Nevertheless, with respect to the latter, one can find a few cases from among a large variety of isolators, for which the basic assumption of a nearly constant Debye temperature appears to be approximately fulfilled (within a limited temperature interval, at least). As a typical case of this type, let us consider the actual magnitudes of the effective (caloric) Debye temperature, ${\Theta _D}\left( T \right)$, of diamond (visualized in Figure 2).

Figure 2 Magnitudes of the effective (caloric) Debye temperatures for diamond, ${\Theta _D}\left( T \right)$, which have been derived using two high-accuracy formulae (Eq. (4.1a) and (4.1b)) from the (isobaric) heat capacities, ${C_p}\left( T \right)$, measured in the low-capacity region by Desnoyer and Morrison [17] (○) and in the high-temperature region by Victor [24] (△) and Dinsdale [25] (□). One can see, in particular, that the ${\Theta _D}\left( T \right)$ values are nearly constant (deviating by less than ±12 K from an average value of about 1857 K) within a rather large interval of about 140 K < $T$ < 700 K.

Concerning the individual sets of ${\Theta _D}\left( T \right)$ values shown in Figure 2, we would like to point out that the ${\Theta _D}\left( T \right)$ data points represented by open circles (○) for the low-temperature region (18 K < $T$ < 278 K) correspond to their original ones (○) shown by Desnoyer and Morrison [17] (in their Figure 1). At variance to these well-known low-temperature ${\Theta _D}\left( T \right)$ data, the different sets of high-temperature ${\Theta _D}\left( T \right)$ data points shown here (in Figure 2) have been derived from the available sets of isobaric heat capacity values, ${C_p}\left( {T \ge 300K} \right)$, that were given by Victor [24] (△) and Dinsdale [25] (□). The corresponding point-by-point transformations have been performed here for the individual ${C_p}\left( T \right)$ vs. ${C_{Vh}}\left( \infty \right)$ ratios, 0 < $\kappa \left( T \right) = {C_p}\left( T \right)/{C_{Vh}}\left( \infty \right)$ < 1, by means of two analytical transformation equations [11] (cf. Eq. (B.5) and (B.3) in [11]), i.e.,

\[ \Theta_{D}^{ }(T) \cong\left[\left(\frac{4 \pi^{4}}{5 \kappa(T)}\right)^{\frac{1}{3}}-g \cdot \kappa(T) \exp \left(-\frac{p}{\kappa(T)}\right)\right] \times T, \ \text { for } \ 0 < \kappa(T) \leq 0.5975, \tag{4.1a} \]

(involving the parameters $g$ = 3.145972 and $p$ = 0.07037359) in combination with

$${\mathrm{\Theta }}_D(T)\cong \sqrt[2]{-35+5\cdot \sqrt[2]{49+56\cdot ((\kappa (T){)}^{-1}-1)}}\times T,$$ $for$ $0.5975\le \kappa (T)\le 1.$

$\left(4.1b\right)$

Note that the maximum inaccuracy resulting from these two of ${\Theta _D}\left( T \right)$ expressions amounts to only 0.13% (at the crossing point, ${\kappa _c}$= 0.5975, between Eqs. (4.1a) and (4.1b)).

Viewing the whole ${\Theta _D}\left( T \right)$ dependence shown in Figure 2, we see, on the one hand, that the sequence of low-temperature data points due to Desnoyers and Morrison’s cryogenic data show a rapid fall from a $T \to 0$ limiting magnitude of about ${\Theta _D}\left( 0 \right) \approx $ 2223 K to magnitudes below 1900 K (in the vicinity of $T \approx $ 150 K). In view of such a rapid change in the $T$-dependent Debye temperature, it is obvious that it is impossible to describe the temperature dependence of the underlying cryogenic heat capacity of diamond, ${C_p}(T < 150\;{\rm{K}})$ (cf. the ${C_p}\left( T \right)$ list presented in the Appendix of [17]), in terms of a Debye function integral, ${\kappa _D}\left( {{\Theta _D}/T} \right)$ (Eq. (2.1)), with a fixed Debye temperature (whichever would be the magnitude chosen for ${\Theta _D}$). At the same time, we see also from Figure 2 that the magnitudes of ${\Theta _D}\left( T \right)$ within a rather large interval (of about 140 K < $T$ < 700 K) are within a relative narrow interval of ${\Theta _D}\left( T \right) \approx $ (1857 ±12) K (corresponding to the variations of ${\Theta _D}\left( T \right)$ values by less than 1.3%). From this observation one can expect that it should be possible to perform a reasonable fitting of the temperature dependence of the heat capacity, within this limited low-to-high temperature region, by a Debye function parameterized by a fixed Debye temperature, ${\Theta _D}$ = const. ≈ 1857 K.

In performing the calculations, for the isochoric (harmonic) part of the heat capacity,

\[ C_{V h}(T) \cong C_{V h}(\infty) \cdot \tilde{\kappa}_{D}\left(\frac{\Theta_{D}}{T}\right) \tag{4.2} \]

and the associated isobaric heat capacities, ${C_p}\left( T \right)$ (see below), we choose, within the present study, the primary interpolation formula proposed above, ${\tilde \kappa _D}\left( x \right)$ (Eq. 3.5), i.e., explicitly

\[ \tilde{\kappa}_{D}\left(\frac{\Theta_{D}}{T}\right) \cong \sqrt[2]{\frac{\exp \left(\frac{x_{o}-\frac{\Theta_{D}}{T}}{\Delta x}\right)+1}{\left(1+\sum_{n=1}^{4} r_{2 n}^{D}\left(\frac{\Theta_{D}}{T}\right)^{2 n}\right) \cdot \exp \left(\frac{x_{o}-\frac{\Theta_{D}}{T}}{\Delta x}\right)+\left(\frac{5}{4 \pi^{4}}\right)^{2}\left(\frac{\Theta_{D}}{T}\right)^{6}}} \tag{4.3} \]

which appears to be directly applicable to a region of 170 K < $T $< 400 K (cf. the corresponding ${\Theta _{Dh}}\left( T \right)$ behavior shown by the dashed curve in Figure 2). At the same time, we see from Figure 2 that the actual (effective) ${\Theta _D}\left( T \right)$ values reach their local maximum in the region between approximately 400– 500 K and subsequently decrease with increasing temperature. Such behavior is well known as a characteristic feature of non-negligible contributions of lattice expansion and anharmonicities to the observable (isobaric) heat capacities, ${C_p}\left( T \right)$. In accordance with our frequently adopted analytical relationship [3,6,10,22,26,32] between isochoric (harmonic) and the associated isobaric heat capacities in solids, we can hence represent the associated ${C_p}\left( T \right)$ dependence by an analytical expression of the form

\[ C_{p}^{ }(T)=C_{V h}(\infty) \cdot \tilde{\kappa}_{D}\left(\frac{\Theta_{D}}{T}\right)\left[1+\tilde{\kappa}_{D}\left(\frac{\Theta_{D}}{T}\right) \cdot \sum_{n=1}^{(\infty)} A_{n} T^{n}\right] \tag{4.4} \]

where the material-specific expansion coefficients ${A_n}$, $n$ = 1,2,…, are quantifying the combined effects of lattice expansion and anharmonicities.

On the basis of Eq. (4.4) (in combination with Eq. (4.3)), we performed a least-mean-square fitting process of the combined set of isobaric heat capacity data points, ${C_p}\left( T \right)$ [17,24,25], that are ranging within an unusually large temperature interval (from 180 K up to 5000 K; cf. Figure 3). As the adjustable parameters, we have obtained an effective (fixed) Debye temperature of ${\Theta _D}$= 1854.8 K in combination with the ${A_n}$ coefficients of ${A_1}$= 2.079 × 10⁻⁵ K⁻¹ and ${A_2}$= 2.421 × 10⁻⁹ K⁻². The respective theoretical ${C_p}\left( T \right)$ and $\rho \left( T \right) = {C_p}\left( T \right)/{T^3}$ dependencies are represented by the dash-dotted curves (-·-·-·) in Figure 3.

Figure 3 Comparison of the alternative fittings (up to 5000 K) of the total set of isobaric heat capacity data, ${C_p}\left( T \right)$, measured by Desnoyer and Morrison [17] (○), Victor [24] (△), and Dinsdale [25] (□), using qualitatively different analytical ${C_p}\left( T \right)$ representations. Good partial (intermediate-to-high-temperature) fittings, at a fixed Debye temperature values of about ${\Theta _D} \approx $1853 K, can be seen to have been realized on the basis of the unprecedented interpolation formula ${\tilde \kappa _D}\left( x \right)$ (Eq. (3.5)), in combination with Eq. (4.4) (for the region 180 K < $T$ < 5000 K), as well as via the truncated high-temperature (high-capacity) algebraic formula (Eq. (A.12)), in combination with Eq. (4.4) (for the region 230 K < $T$ < 5000 K, at least). An exemplary high-accuracy fit for the whole ${C_p}\left( T \right)$ dataset (solid curve, due to Eq. (5.10)) had been obtained on the basis of the representative hybrid model displayed in Sec. 5. Also shown in the figure is the corresponding high-temperature behavior of the associated isochoric (harmonic) heat capacity, ${C_{Vh}}\left( T \right) = \;3R \cdot {\kappa _P}\left( T \right)$ (dashed curve, according to Eq. (5.4)).

It is still instructive to consider a somewhat rougher, alternative approximation for the isochoric part of the heat capacity, ${C_{Vh}}\left( T \right)$, and the associated isobaric heat capacity, ${C_p}\left( T \right)$. This is obtained by replacing the approximate (interpolated) Debye function, ${\tilde \kappa _D}\left( {{\Theta _D}/T} \right)$ (in Eq. (4.2) and (4.4)), simply by the truncated high-temperature (algebraic) formula, ${\kappa }_{Dh}^{\left(m_P^{}= 4\right)}\left({\Theta }_D/T\right)$ (Eq. (3.1)), for the Debye function, i.e., explicitly

$$C_{Vh}(T)\to C_{Vh}(\mathrm{\infty })\cdot {\kappa }_{Dh}^{(m_P^{}=4)}\left(\frac{{\Theta }_D}{T}\right)=\frac{C_{Vh}(\mathrm{\infty })}{\sqrt[2]{1+\sum _{n=1}^4{r}_{2n}^D{\left(\frac{{\mathrm{\Theta }}_D}{T}\right)}^{2n}}}.(4.5)$$

The corresponding theoretical $\rho \left( T \right) \to {C_{Vh}}\left( T \right)/{T^3}$ dependence is represented by the dash-double-dot curve (-··-··-··) in the inset of Figure 3.

5. Alternative High-accuracy Fit due to the Representative Hybrid Model

The actual variation of the Debye temperature, ${\Theta _D}\left( T \right)$, of the wide-band-gap as well as semiconductor materials is well-known to usually undergo rapid falls from their limiting ($T \to 0$) magnitudes, ${\Theta _D}\left( 0 \right)$, toward significantly lower minimum values, ${\Theta _{D\;min}}$, which used to be located at temperatures of the order ${T_{min}} \approx {\Theta _D}\left( 0 \right)/12$, (in accordance with Figure 2 for diamond; see also the cryogenic ${\Theta _D}\left( T \right)$ curves shown for Si an Ge in [2,3], for III-V semiconductors in [4,5,6], and for II-VI semiconductors in [7,8,9,10,11]). Thus, it is impossible to simulate the temperature dependences of the measured ${C_{p/Vh}}\left( T \right)$ dependences in cryogenic (liquid-hydrogen) regions, $0 < T < {\Theta _D}\left( 0 \right)/12$, by means of a Debye function integral (2.1) with a constant Debye temperature (whichever magnitude would have been chosen for a fixed ${\Theta _D}$ value). A conventional way of reasonable fittings for at least partial (strongly truncated) sets of cryogenic ${C_{Vh}}\left( T \right)$ data is well known to be given in terms of an odd-order Taylor series expansion [6,12,14,17,27,28,29,30,31,32],

\[ C_{V h}(T)=c_{3} T^{3}+c_{5} T^{5}+c_{7} T^{7}+\ldots \tag{5.1} \]

(cf., in particular, the numerous applications of this simple analytical model to Group-IV, III-V, and II-VI materials in [12]). However, this truncated odd-order series expansion is continually found to be applicable only to narrow intervals of about $0 < T < {\Theta _D}\left( 0 \right)/\left( {25\; \pm 10} \right)$ [12]. (Note that, in cases of non-negligible contributions of an electronic system to the heat capacity in the liquid-helium-region, an additional linear term [11,28,29], ${C_{el}}\left( T \right) \to {c_1}T$, still needs to be included in Eq. (5.1)).

Suitable incorporation of this commonly observable cryogenic ${C_{Vh}}\left( T \right)$ behavior (Eq. (5.1)) into theoretical heat capacity models thus represents the necessary condition for their eventual applicability to the whole $T$-regions of practical interest, i.e., from the $T \to 0$ limit up to very high temperatures (up to melting points, if necessary). This condition is well-fulfilled by two qualitatively different types of duly general analytical models, namely, the representative hybrid model [3,10,11,26], on the one hand, and one of the two (essentially equivalent) versions of an analytical Non-Debye formula [6,32], on the other hand. By comparing the alternative fittings of various comprehensive ${C_p}\left( T \right)$ datasets by these two types of theoretical heat capacity models, we have found that the representative hybrid model [3,10,11,26] is capable of providing the finest numerical simulations of given ${C_p}\left( T \right)$ datasets.

This model is based on a hybrid ansatz for the phonon density of states (PDOS) function, $g_P^{ }\left(\epsilon \right)$, which consists of three discrete (Einstein) peaks, ${w_{Ei}}\delta \left( {\varepsilon - {\varepsilon _i}} \right)$ (where $0 < {\varepsilon _1} < {\varepsilon _2} < {\varepsilon _3}$), in combination with two continuous (quadratic and quartic) power function sections, ${w_{C1}} \cdot 3{\varepsilon ^2}/\varepsilon _1^3$ and ${w_{C2}} \cdot 5{\varepsilon ^4}/\varepsilon _1^5$ (for the low energy region from 0 up to the first Einstein peak, ${\varepsilon _1}$). This hybrid model thus involves a global normalization condition of the form

\[ w_{C 1}+w_{C 2}+w_{E 1}+w_{E 2}+w_{E 3}=1 \tag{5.2} \]

for the total set of weighting factors. Accordingly, the moments for this theoretical hybrid PDOS spectrum are given by [10,11,26]

\[ \mu_{P}^{(n)}=\left(\frac{3}{3+n} w_{C 1}+\frac{5}{5+n} w_{C 2}\right) \varepsilon_{1}^{n}+\sum_{i=1}^{3} w_{E i} \varepsilon_{i}^{n} \text { for }-3 < n < 0 \text { and } n>0 \tag{5.3} \]

Within the framework of this model, the temperature dependence of the isochoric (harmonic) part of the heat capacity, ${C_{Vh}}\left( T \right) = {C_{Vh}}\left( \infty \right) \cdot {\kappa _P}\left( T \right)$, is described by the respective (duly normalized) heat capacity shape function, ${\kappa _P}\left( T \right) \le 1$ [26],

\[\kappa_{P}(T)=w_{C 1} \cdot \kappa_{C 1}\left(\frac{\Theta_{1}}{T}\right)+w_{C 2} \cdot \kappa_{C 2}\left(\frac{\Theta_{1}}{T}\right)+\sum_{i=1}^{3} w_{E i} \cdot\left(\frac{\Theta_{i}}{2 T \sinh \left(\frac{\Theta_{i}}{2 T}\right)}\right)^{2} \tag{5.4} \]

Here, ${\Theta _i} = {\varepsilon _i}/{k_B}$ denotes the respective Einstein temperatures pertaining to the individual discrete peaks. The normalized contributions of the quadratic and quartic power function components are given in terms of the ratio ${x_1} = {\Theta _1}/T$, by the integrals [11]

\[ \kappa_{C 1}\left(x_{1}^{}\right)=\frac{3}{4 x_{1}^{3}} \int_{0}^{x_{1}^{}} d z \frac{z^{4}}{\left(\sinh \left(\frac{Z}{2}\right)\right)^{2}} \leq 1 \text { and } \kappa_{C 2}\left(x_{1}^{}\right)=\frac{5}{4 x_{1}^{5}} \int_{0}^{x_{1}^{}} d z \frac{z^{6}}{\left(\sinh \left(\frac{Z}{2}\right)\right)^{2}} \leq 1. \tag{5.5} \]

Comparing the integral representation for the ${\kappa _{C1}}\left( {{x_1}} \right)$ function (Eq. (5.5)) with Eq. (2.3) for the conventional Debye function integral, ${\kappa _D}\left( x \right)$, we see that, that the dependence of ${\kappa _{C1}}\left( {{x_1}} \right)$ on the ratio ${x_1} = {\Theta _1}/T$ is identical to the dependence of ${\kappa _D}\left( x \right)$ on $x = {\Theta _D}/T$. This means that the asymptotic ($T \to 0$) temperature dependence of ${\kappa _{C1}}\left( {{x_1}} \right)$ for sufficiently low temperatures, ${x_1} = {\Theta _1}/T > 12$, tends to the well-known limiting cubic power function, ${\kappa _{C1}}\left( {{\Theta _1}/T} \right) \to \left( {4{\pi ^4}/5} \right) \cdot {\left( {T/{\Theta _1}} \right)^3}$ (in accordance with [11,26]). Analogously, for the asymptotic behavior of the second component, ${\kappa _{C2}}\left( {{x_1}} \right)$, an asymptotic ($T \to 0$) temperature dependence is obtained in the form of a quintic power function, ${\kappa _{C2}}\left( {{\Theta _1}/T} \right) \to \left( {80{\pi ^6}/21} \right) \cdot {\left( {T/{\Theta _1}} \right)^5}$ [11,26]. For regions of very low temperatures, $T < {\Theta _1}/12$, this hybrid model reduces to combinations of unambiguously fixed cubic and quintic $T$-power function terms,

\[ C_{V h}(T) \rightarrow C_{V h}(\infty)\left(w_{C 1} \cdot \frac{4 \pi^{4}}{5 \Theta_{1}^{3}} \cdot T^{3}+w_{C 2} \cdot \frac{80 \pi^{6}}{21 \Theta_{1}^{5}} \cdot T^{5}\right) \tag{5.6} \]

The latter asymptote is in qualitative accordance with the occurrence of the respective two lowest-order power terms occurring in the well-known empirical low-temperature heat capacity formula given by Eq. (5.1), ${C_{Vh}}\left( T \right) \to {c_3}{T^3} + {c_5}{T^5}$. Comparing the latter with Eq. (5.6), we see that the empirical expansion coefficients, ${c_3}$ and ${c_5}$, can be represented in terms of the weighting factors ${w_{C1}}$, ${w_{C2}}$, and the characteristic phonon temperature, ${\Theta _1}$, due to the first (lowest) Einstein phonon peak, by the following relations:

\[ c_{3}=C_{V h}(\infty) \cdot w_{C 1} \frac{4 \pi^{4}}{5 \Theta_{1}^{3}} \text { and } c_{5}=C_{V h}(\infty) \cdot w_{C 2} \frac{80 \pi^{6}}{27 \Theta_{1}^{5}} \tag{5.7} \]

Further, we observe that the coefficient ${c_3}$ is well-known to be unambiguously connected with the $T \to 0$ limiting magnitude of the Debye temperature, ${\Theta _D}\left( 0 \right)$, by the relation [12]

\[ \Theta_{D}(0)=\left(\frac{\left(\frac{4}{5}\right) \pi^{4} C_{V h}(\infty)}{c_{3}}\right)^{\frac{1}{3}}, \text { which is equivalent to } c_{3}=C_{V h}(\infty) \cdot \frac{4 \pi^{4}}{5\left(\Theta_{D}(0)\right)^{3}}. \tag{5.8} \]

Comparing the latter expression for ${c_3}$ (in (5.8)) with the preceding one (in (5.7)), we see that the ratio between ${\Theta _1}$ and ${\Theta _D}\left( 0 \right)$ is unambiguously determined by the weighting factor, ${w_{C1}}$, for the cubic component [11], i.e.,

\[ \left(\frac{\Theta_{1}}{\Theta_{D}(0)}\right)^{3}=w_{C 1}, \text { which is equivalent to } \Theta_{D}(0)=\frac{\Theta_{1}}{\sqrt[3]{w_{C 1}}} \tag{5.9} \]

Finally, we consider the measured (isobaric) temperature dependences, ${C_p}\left( T \right)$, for higher temperatures, which use to be co-determined by anharmonicity effects. The entire ${C_p}\left( T \right)$ curve, from 0 up to very high temperatures, can be represented (in analogy to Eq. (4.4)) by an expression of the form [11,26]

\[ C_{p}(T)=C_{V h}(\infty) \cdot \kappa_{P}(T) \cdot\left[1+\kappa_{P}(T) \cdot\left(A_{1} T+A_{2} T^{2}+\ldots\right)\right] \tag{5.10} \]

We have performed a complete high-accuracy fit of the combined ${C_p}\left( T \right)$ datasets under study [17,24,25] for diamond (as shown by the solid curve in Figure 3). As a result of the multi-parameter least-mean-square fitting process, we have obtained the following values for the effective Einstein temperatures of the three discrete phonon peaks under consideration:

\[ \Theta_{1}=778.5 \mathrm{~K}, \Theta_{2}=1108.5 \mathrm{~K}, \text { and } \Theta_{3}=1733.6 \mathrm{~K} \tag{5.11a} \]

in combination with the weighting factors

\[ w_{C 1}=0.04294, w_{C 2}=0.01425, w_{1}=0.06514, w_{2}=0.30894, w_{3}=0.56872, \tag{5.11b} \]

and the associated ${A_n}$ coefficients

\[ A_{1}=2.158 \times 10^{-5} \mathrm{~K}^{-1} \text { and } A_{2}=2.451 \times 10^{-9} \mathrm{~K}^{-2} \tag{5.11c} \]

The corresponding theoretical ${C_p}\left( T \right)$ curve (up to 5000 K) and the associated $\rho \left( T \right) = {C_p}\left( T \right)/{T^3}$ dependence is represented by the solid curves (——) in Figure 3.

Furthermore, we have shown above (in Figure 2) the corresponding theoretical ${\Theta _D}\left( T \right)$ and ${\Theta _{Dh}}\left( T \right)$ dependences, which have been obtained via the given transformation equations ((Eqs. (4.1a) and (4.1b)) for the calculated ratios $\kappa \left( T \right) \to {\kappa _P}\left( T \right) = {C_{Vh}}\left( T \right)/{C_{Vh}}\left( \infty \right)$ (Eq. (5.4)) and $\kappa \left( T \right) \to {C_p}\left( T \right)/{C_{Vh}}\left( \infty \right)$ (Eq. (5.10)), respectively.

6. Discussion

Above all, in this paper, we have developed an unprecedented algebraic formula, ${\kappa }_{Dh}^{\left(m_P^{ }\right)}\left(x\right)$ (Eq. (A.12)), for the high-temperature behavior of the Debye function. Using, in particular, a suitably truncated version of the latter, ${\kappa }_{Dh}^{\left(m_P^{}=4\right)}\left(x\right)$ (in Eq. (4.5)), we were able to realize a good approach to the true ${C_{Vh}}(T > {T_c})$ dependence for diamond (at least within an interval from about 230 K to about 600 K (cf. the dash-double-dot curve for the $\rho \left( T \right) \approx {C_{Vh}}\left( T \right)/{T^3}$-dependence, in the inset of Figure 3). A somewhat closer approach to the true $\rho \left( T \right)$ behavior within the region 180 K to 230 K could be achieved when we used, alternatively, the high-to-low-temperature interpolation formula, ${\tilde \kappa _D}\left( {{\Theta _D}/T} \right)$ (Eq. (4.3)), for the fitting. However, it is a matter of principle that the corresponding ${\rho _D}\left( T \right)$ dependence (dash-dot curve) exhibits a plateau (approximate constancy) from 0 up to approximately 150 K. This typical Debye function feature is, of course, in striking contrast to the well-known non-Debye behavior (see Eq. (5.1)) of the measured ${C_p}\left( T \right)$ dependences in the corresponding cryogenic regions, $0 < T < {T_c} \approx {\Theta _D}\left( 0 \right)/12$.

6.1 Additional Quantitative Information Implied by the High-accuracy Fit via the Representative Hybrid Model

As described in Section 5, we have performed a very accurate fit, from approximately 20 K up to 5000 K (cf. the solid curves in Figure 3), for the whole combination of the partial ${C_p}\left( T \right)$ datasets given for low and/or high temperatures [17,24,25]. This least-mean-square fit has been realized on the basis of the representative hybrid model [10,11,26], which is represented here by Eq. (5.10) in combination with Eq. (5.4). From the adjusted set of the model-specific parameter values (5.11), it follows that the average deviations of the theoretical ${C_p}\left( T \right)$ values from the respective experimental points amount to approximately ±0.45%. The unusually small size of these average deviations can be considered to be a confirmation of the high quality of the analytical apparatus used (in Sec. 5) and also of the experimental datasets [17,24,25] considered in this study.

On the basis of the resulting parameter values (Eq. (5.11)), one can readily evaluate the characteristic quantities associated with the $T \to 0$ limiting behavior of the fitted $\rho \left( T \right)$ and ${\Theta _D}\left( T \right)$ curves. Using the first relation in Eq. (5.7), a limiting $\rho \left( {T \to 0} \right)$ value of ${c_3} = \rho \left( 0 \right) = $ 0.17694 µJ K⁻⁴ mol⁻¹ (as indicated by a horizontal $\rho \left( {T \to 0} \right)$ line in the inset of Figure 3) is obtained. According to the first relation in Eq. (5.8), this ${c_3}$-value corresponds to the $T \to 0$ limiting value of ${\Theta _D}\left( 0 \right) \cong $ 2223 K (as indicated by a horizontal line in Figure 2). (Note that the same ${\Theta _D}\left( 0 \right)$ value can be obtained using the second relation in Eq. (5.9)). This value is found to be in good agreement with the frequently quoted caloric ${\Theta _D}\left( 0 \right)$ value of about 2220 K [33,34,35] for diamond.

The maximum of the $\rho \left( T \right)$ curve, ${\rho _{max}} = \rho \left( {{T_c}} \right) = $ 0.2995 µJ K⁻⁴ mol⁻¹, has been found to be located at the characteristic point of ${T_c} = \;$174 K, in the vicinity of which the corresponding ${C_p}\left( T \right)$ the curve shows a certain cubic dependence, ${C_p}\left( {T \sim {T_c}} \right) \to \rho \left( {{T_c}} \right) \cdot {T^3}$. However, the latter is to be clearly distinguished from the well-known $T \to 0$ (limiting) cubic dependence based on Debye’s theory [1], ${C_V}\left( {T \to 0} \right) \to {c_3}{T^3}$ (cf. Eq. (5.1)). The qualitative and quantitative differences between these two (local) cubic $T$-dependences in the cryogenic region are due to the circumstance that the magnitudes of the respective proportionality factors, i.e., $\rho \left( {{T_c}} \right)$ in comparison with $\rho \left( 0 \right) = {c_3}$, are markedly different, $\rho \left( {{T_c}} \right)/\rho \left( 0 \right) = $ 1.693.

Other interesting by-products of the fit via the hybrid model are the moments, $\mu _P^{\left( n \right)}$ (Eq. (5.3)), of the underlying theoretical (hybrid) PDOS model spectrum under consideration (see Section 5). Of particular interest for an assessment of the degree of phonon dispersion are the first-and second-order moments, the magnitudes of which follow from the fitted parameter set (Eq. (5.11)) to be $\mu _P^{\left( 1 \right)} = \;$121.795 meV and $\mu _P^{\left( 2 \right)} = $ (126.35 meV)². Thus, the characteristic dispersion coefficient, ${\Delta _P}$ [36,37] has a magnitude of

$${\mathrm{\Delta }}_P\equiv \frac{\sqrt[2]{{\mu }_P^{(2)}-{\left({\mu }_P^{(1)}\right)}^2}}{{\mu }_P^{(1)}}=0.276.$$

$\left(6.1\right)$

This value is in very good agreement with the earlier values of ${\Delta _P} = \;$0.275 ±0.005, which had been estimated in previous works [36,37] on the basis of the calculated PDOS spectra published for diamond by four different authors (see Table 1 in [36] and Table 1 in [37]). This ${\Delta _P}$-value (Eq. (6.1)) indicates an unusually low degree of dispersion since it is only slightly higher than its theoretical counterpart implied by Debye’s original (purely quadratic) PDOS model function [1], ${g_D}\left( \varepsilon \right) = \;3{\varepsilon ^2}/\varepsilon _D^3$ (for the interval $0 < \varepsilon < {\varepsilon _D}$). For the first and second-order moments, this simple model invoked by Debye gives the values $\mu _D^{\left( 1 \right)} = \left( {3/4} \right) \cdot {\varepsilon _D}$ and $\mu _D^{\left( 2 \right)} = \left( {3/5} \right) \cdot \varepsilon _D^2$. Inserting these special $\mu _D^{\left( 1 \right)}$ and $\mu _D^{\left( 2 \right)}$ values into the general expression (Eq. (6.1)) for the dispersion coefficients gives for this conventional model, a magnitude of [26,37]

$${\mathrm{\Delta }}_D\equiv \frac{\sqrt[2]{{\mu }_D^{(2)}-{\left({\mu }_D^{(1)}\right)}^2}}{{\mu }_D^{(1)}}=\sqrt[2]{\frac{1}{15}}=0.258199.$$

$\left(6.2\right)$

Thus, we see that the material-specific magnitude of the dispersion coefficient, ${\Delta _P}$ (Eq. (6.1)), for diamond is only approximately 7% higher than its model-specific counterpart, ${\Delta _D}$ (Eq. (6.2)), implied by Debye’s original model [1]. The relatively similar values of ${\Delta _P}$ to ${\Delta _D}$ is the obvious reason of the relatively good functioning of Debye’s very special heat capacity shape function, ${\kappa _D}\left( x \right)$ (Eq. (2.1)), for an approximate description of the ${C_{Vh}}(T > {T_c})$ dependence (as shown in Figure 3).

Another important theoretical quantity is the high-temperature limiting $\left( {T \to \infty } \right)$ value of the Debye temperature for harmonic lattice oscillations, ${\Theta _{Dh}}\left( {T \to \infty } \right)$. This quantity is well-known to be connected unambiguously with the second moment, $\mu _P^{\left( 2 \right)}$, by the equation [26]

$${\Theta }_{Dh}(\mathrm{\infty })=\frac{\sqrt[2]{\left(\frac{5}{3}\right){\mu }_P^{(2)}}}{k_B}=1893.0\mathrm{K}.$$

$\left(6.3\right)$

The magnitude of this limiting (theoretical) Debye temperature value, ${\Theta _{Dh}}\left( \infty \right)$, is indicated in Figure 2 by an arrow.

6.2 Assessment of the Incomplete Fittings by Novel Formulae Involving a Fixed Debye Temperature

Taking the high-accuracy fit performed for a diamond based on the representative hybrid model (in Section 5) into consideration, let us now discuss in more detail the possible practical usability of two versions of the unprecedented Debye function formulas displayed in Section 3.

Firstly, concerning the limited fit of the experimental ${C_p}\left( T \right)$ datasets, from 180 K to 5000 K, done using the interpolation formula, ${\tilde \kappa _D}\left( {{\Theta _D}/T} \right)$ (Eq. (4.3)), in combination with Eq. (4.4), we are concerned with average deviations of ±0.54% between the theoretical and experimental ${C_p}\left( T \right)$ values (which are only slightly higher than those associated with the abovementioned fit by the hybrid model). Accordingly, the approximate theoretical ${C_p}\left( T \right)$ and $\rho \left( T \right) = {C_p}\left( T \right)/{T^3}$ dependences, which are represented by the dash-dot curves in Figure 3, are almost coinciding, for any $T$ > 180 K, with their counterparts due to the hybrid model (solid curves). The simple model, based on the primary interpolation formula (Eq. (3.5)), might thus be of some practical use with respect to the materials having similarly low ${\Delta _P}$-values, such as diamond (Eq. (6.1)). This interpolation model could, in particular, be of use within the framework of thermo-chemistry, where one is mainly interested in fine numerical descriptions of ${C_p}\left( T \right)$ dependences within the high-$T$ region, $T $> 298.15 K. In addition, from the inset of Figure 3, it can be seen that this model, which exhibits a low-temperature plateau behavior of the respective, approximate ${C_p}\left( T \right)/{T^3}$ function (dash-dot-curve), is absolutely inapplicable to $T < {T_c}$ = 174 K.

With regards to the alternative fit of the experimental ${C_p}\left( T \right)$ datasets, from 230 K to 5000 K, via the truncated high-temperature Debye function formula, ${\kappa }_{Dh}^{\left(m_P^{}=4\right)}\left(x\right)$ (Eq. (3.1)), average deviations of ±0.74% between the theoretical and experimental ${C_p}\left( T \right)$ values (being higher by a factor of approximately 1.6 than those due to the hybrid model) are observed within this limited region. in addition, we see from the respective, approximate ${C_p}\left( T \right)/{T^3}$ function (represented by a dash-double-dot curve, in the inset of Figure 3) that this simplified (algebraic) model nevertheless provides a relatively close fit (involving underestimation of the experimental values by less than 8%) in the region from approximately 80 K to 230 K.

6.3 Sub-quartic Behavior of the Low-temperature Heat Capacity Dependence

Such a relatively close fit to the actual (experimentally observed) $T$-dependence, within a region of approximately ${T_c}/2 < T < {T_c}$, is not surprising in view of the circumstance that the $T$-dependence of the truncated Debye function formula, ${\kappa }_{Dh}^{\left(m_P^{}=4\right)}\left({\Theta }_D/T\right)$ (Eq. (3.1)), is tending in the $T \to 0$ limit to a certain (fictive) quartic ${C_{Vh}}\left( T \right)$-dependence (cf. (Eq. (4.5)), $C_{Vh}\left(T\right)\to \left(C_{Vh}\left(\infty \right)/\sqrt[2]{r_8^D}\right)\cdot {\left(T/{\Theta }_D\right)}^4$. This peculiar feature is at least in qualitative accordance with the global expectation [38], according to which, in the cryogenic region below ${T_c}$, the effective exponents, $\eta $, of the power function tangents, ${C_p}\left( T \right) \propto {T^\eta }$, should range somewhere between 3 and 4 (for a wide-band-gap material such as diamond; cf. Figure 7 in [38]). In order to estimate the maximum value of this exponent, ${\eta _m}$, let us consider for the vicinity of the point of inflexion, ${T_i}$ = 99.3 K, of the fitted $\rho \left( T \right)$ curve, an approximation by a linear tangent [38],

\[ \rho\left(T \sim T_{i}\right) \rightarrow \rho\left(T_{i}\right)+\rho^{(1)}\left(T_{i}\right) \cdot\left(T-T_{i}\right) \equiv d_{0}+d_{1} \cdot T \tag{6.4} \]

where ${d_0} \equiv \rho \left( {{T_i}} \right) - {\rho ^{\left( 1 \right)}}\left( {{T_i}} \right) \cdot {T_i} = 0.104745$ µJ K⁻⁴ mol⁻¹ and ${d_1} = {\rho ^{\left( 1 \right)}}\left( {{T_i}} \right)$ = 0.13911×10⁻² µJ K⁻⁵ mol⁻¹ (see the dotted line in the inset of Figure 3). This tangent (Eq. (6.4)) corresponds to the local heat capacity dependence in the form of a combination of a cubic with a quartic term [38],

\[ C_{p}\left(T \sim T_{i}\right)=\rho\left(T \sim T_{i}\right) \cdot T^{3} \rightarrow d_{0} \cdot T^{3}+d_{1} \cdot T^{4} \tag{6.5} \]

Using the general expression [38] for the $T$-dependence of the effective power function exponents,

\[ \eta(T)=\frac{d \log \left(C_{p}(T)\right)}{d \log (T)}=\left(\frac{d C_{p}(T)}{d T}\right) \cdot\left(\frac{T}{C_{p}(T)}\right) \tag{6.6} \]

from the estimated ${C_p}\left( {T \sim {T_i}} \right)$ dependence (Eq. (6.5)) we obtain a magnitude of

\[ \eta_{m}=\eta\left(T_{i}\right)=\frac{3 d_{0} T_{i}^{3}+4 d_{1} T_{i}^{4}}{d_{0} T_{i}^{3}+d_{1} T_{i}^{4}}=3+\frac{d_{1} T_{i}}{d_{0}+d_{1} T_{i}}=3.5687 \tag{6.7} \]

for the effective (maximum) power function exponent in the vicinity of the point of inflection (as indicated by a solid triangle in the inset of Figure 3). Thus, in the case of diamond, we are concerned with a typical sub-quartic behavior of the cryogenic ${C_p}\left( T \right)$ dependence. This is in accordance with our preliminary assessment shown for diamond in Figure 7 of an earlier study [38].

6.4 Comparison with a Recently Proposed Multi-Debye-function Model

Of interest might be a comparison with a recently published model [39] for heat capacity fitting, from cryogenic to high temperatures, on the basis of a certain interpolation model consisting of a combination of three Debye function components, for low-temperature regions, with the popular (thermo-chemical) Maier-Kelley approximation [40], for high temperatures. Concerning the case of diamond, one can see from a comparison, e.g., of the ${C_{Vh}}\left( T \right)$ curve shown in Figure 7 from an earlier study [39] with the present counterpart (shown by the dashed curve in Figure 3) that their behavior appears to be nearly in mutual agreement (within a limited range of 200 K to 2000 K, at least). At first sight, such approximate equality might appear somewhat surprising since one could rather expect more pronounced qualitative and/or quantitative differences between the present interpretation in terms of a single Debye function (Eq. 4.3)), for the region from 200 K to approximately 600 K, and the numerical simulation by Vassiliev and Taldrik [39], who employed a combination of three Debye functions for the same purpose. Yet, an explanation of this peculiarity can be readily given by observing that there are only relatively small differences between the individual ${\Theta _{Di}}$ values (${\Theta _{D1}}$, ${\Theta _{D2}}$, and ${\Theta _{D3}}$) estimated by Vassiliev and Taldrik [39] (which are listed, together with the respective weighting factors, ${A_i}$, e.g., under number 1b in Table 6 of their article). The average magnitude, $\overline\Theta _{D}=A_1\Theta_{D1}+A_2\Theta_{D2}+A_3\Theta_{D3}$, of this set of ${\Theta _{Di}}$ values amounts to $\overline\Theta_{D}$ = 1896 K. The latter is by only 2.2% higher than our single Debye temperature value of ${\Theta _D}$= 1855 K, which resulted from the present fitting using Eq. (4.2) and/or Eq. (4.5). Furthermore, we see that the individual ${\Theta _{Di}}$ values estimated by Vassiliev and Taldrik [39] differ from their average value, $\overline\Theta_{D}$, by less than 4%. The relatively small differences between the individual ${\Theta _{Di}}$ values provide an explanation for the similar success of our alternative models, which are involving only a single effective Debye temperature for the two comparable numerical simulations of the ${C_p}(T > {T_c})$ data presented in Section 4.

6.5 Possible Usefulness of a Properly Adopted Non-Debye Interpolation Formula

Satisfactory numerical fitting to the comprehensive ${C_p}\left( T \right)$ datasets for the high-dispersion wide-band-gap materials GaN and ZnO [36,37], from very low up to high temperatures, had been performed for the first time in a previous study [32] on the basis of a novel algebraic non-Debye formula. A suitable alternative version of the latter had been displayed and used in another work [6] for comprehensive fitting (from liquid-helium temperatures up to the melting points of the materials) of the combined ${C_p}\left( T \right)$ datasets, that are available for a series of Group III-V (high-dispersion [37]) semiconductor materials.

We have found that an appropriate specialization of the preceding analytical non-Debye model [6] with respect to its possible applications to low-dispersion materials, such as diamond, can readily be established by adopting an algebraic expression of the interpolation type

\[ \kappa_{D ; N D}(T)=\frac{1+\frac{c_{5}^{ }}{c_{7}^{ } T^{2}}+\frac{c_{3}^{ }}{c_{7}^{ } T^{4}}}{ \sqrt[2]{1+\sum_{n=1}^{3} r_{2 n}^{D}\left(\frac{\Theta_{D}}{T}\right)^{2 n}+r_{8}^{N D}\left(\frac{\Theta_{D}}{T}\right)^{8}+\left(\frac{C_{V h}(\infty)}{c_{7}^{ } T^{7}}\right)^{2}}} \tag{6.8} \]

for the isochoric heat capacity shape function (the structure of which is analogous to the structure of the original non-Debye heat capacity shape function (A3) considered in [6]). Envisaging Eq. (6.8) for calculations of the respective isochoric heat capacities,

\[ C_{V h}(T)=C_{V h}(\infty) \cdot \kappa_{D ; N D}(T) \tag{6.9} \]

we see, on the one hand, that at sufficiently low temperatures, where the contributions of the first four $ \propto {T^{ - 2n}}$ terms ($n$ = 1 to 4, in the denominator of Eq. (6.8)) are small compared to the contribution of the $ \propto {T^{ - 14}}$ term, the ${C_{Vh}}\left( T \right)$ expression in Eq. (6.9) approaches asymptotically to the familiar odd-order Taylor series expansion in Eq. (5.1). On the other hand, at sufficiently high temperatures, the numerator tends to unity and the denominator approaches asymptotically to the truncated high-temperature expression given by Eq. (4.5). However, a characteristic deviation in the latter case is due to the circumstance that the fourth expansion coefficient, $r_8^{ND}$, for the $ \propto {T^{ - 8}}$ term, must be admitted to deviate more or less significantly from its (originally fixed) Debye-function-related predecessor, $r_8^D$ (Table 1). This generalization is inevitable for assuring a smooth change between the respective asymptotic low-and high-temperature behavior.

Within the framework of this unprecedented interpolation model, the isobaric heat capacities, ${C_p}\left( T \right)$, are represented by an analytical function of the usual type:

\[ C_{p}(T)=C_{V h}(\infty) \cdot \kappa_{D ; N D}(T) \cdot\left[1+\kappa_{D ; N D}(T) \cdot\left(A_{1} T+A_{2} T^{2}+\ldots\right)\right] \tag{6.10} \]

(in analogy to Eqs. (4.4) and (5.10)). Performing, tentatively, an alternative fit of the whole ${C_p}\left( T \right)$ dataset under study for diamond using Eq. (6.10) (in combination with Eq. (6.8)) we have obtained the cryogenic parameter values ${c_3}$= 0.17907 µJ K⁻⁴ mol⁻¹, ${c_5}$ = 0.32631×10⁻⁵ µJ K⁻⁶ mol⁻¹, and ${c_7}$= 0.40863×10⁻⁹ µJ K⁻⁸ mol⁻¹ (cf. [12]) in combination with a fitted non-Debye expansion coefficient of $r_8^{ND}$= 1.862 × 10⁻⁶, an effective Debye temperature of ${\Theta _D}$= 1904.8 K, and the associated anharmonicity coefficients of ${A_1}$= 1.799×10⁻⁵ K⁻¹ and ${A_2}$= 3.278 × 10⁻⁹ K⁻². In particular, assessing the presently obtained ${\Theta _D}$ value, we see that it is by only approximately 0.6% higher than the theoretical high-temperature limiting $\left( {T \to \infty } \right)$ value of the Debye temperature for harmonic lattice oscillations, ${\Theta _{Dh}}\left( {T \to \infty } \right)$ = 1893 K. (cf. Eq. (6.3)). Thus, the unprecedented algebraic model function of the interpolation type (Eq. (6.10)) appears to offer us a chance of giving, as an important by-product of fitting, even the correct order-of-magnitude of the material-specific ${\Theta _{Dh}}\left( {T \to \infty } \right)$ values. Furthermore, we see that the fitted magnitude of the Non-Debye expansion coefficient, $r_8^{ND}$, is a factor of ~3 higher than its preceding counterpart, $r_8^D$, due to the original Debye function expansion (cf. Table 1). This considerable difference is obviously due to strongly increasing non-Debye effects in the vicinity of ${T_c} \approx {\Theta _D}\left( 0 \right)/12$.

Concerning the associated estimation of the $T \to 0$ limiting value of the Debye temperature, we see that the alternatively fitted ${c_3}$-value of ${c_3}$= 0.17907 µJ K⁻⁴ mol⁻¹ corresponds (according to Eq. (5.8)) to a ${\Theta _D}\left( 0 \right)$ value of about 2214 K. The latter is approximately 0.3% lower than the apparently more exact value of 2223 K (as resulting from the refined fit by the hybrid model used for the calculations presented in Section 5). Further, we would like to mention that the latter fit via the unprecedented algebraic interpolation model function (Eq. (6.10)) implies the maximum value of the $\rho \left( T \right)$ curve to be ${\rho _{max}} = \rho \left( {{T_c}} \right) = $ 0.3042 µJ K⁻⁴ mol⁻¹, (which is by ~1.5% higher than its counterpart due to the hybrid model used for the calculations presented in Section 5). The respective ${T_c}$-position of the maximum, ${\rho _{max}} = \rho \left( {{T_c}} \right)$, resulting from Eq. (6.10) is located at ${T_c} \approx $ 164 K (which is by ~6% lower than the ${T_c}$-position estimated via the hybrid model used in Section 5). Finally, we would like to mention that the average deviations of the theoretical ${C_p}\left( T \right)$ values resulting from the fit of the experimental points under consideration using Eq. (6.10) (in combination with Eq. (6.8)), amount to approximately ±1%, which corresponds to the same order of magnitude as the typical experimental uncertainties in the low-temperature region. In view of these still moderate deviations between the theoretically fitted and experimental ${C_p}\left( T \right)$ values, one can conclude that the latter version (Eq. (6.10), in combination with (6.8)) of an algebraic non-Debye interpolation formula can be successfully applied to numerical fittings of comprehensive ${C_p}\left( T \right)$ datasets, from the liquid helium region ($T \to 0$) up to very high temperatures, provided that the dispersion coefficient, ${\Delta _P}$ (Eq. (6.1)), of a wide-band-gap material under consideration is not significantly higher than those estimated for diamond and c-BN [36,37].

7. Concluding Remarks

We have displayed within this study a variety of structurally relatively simple interpolation formulae that are convenient for fitting the measured low-to-high temperature heat capacity datasets of certain wide-band-gap materials that are characterized by a sufficiently low degree of phonon dispersion (0.258 < ${\Delta _P}$ < 1/3; cf. Eq. (6.1) and (6.2)).

A characteristic empirical feature of such low-dispersion materials is the quasi-plateau-behavior of the effective (caloric) Debye temperature, ${\Theta _D}\left( T \right)$, which is observed in $T$-regions from ${\Theta _D}\left( 0 \right)/12$ up to an order of ${\Theta _D}\left( 0 \right)/3$ (cf. Figure 2 for diamond; see also Figure 2 in [19] for c-BN).

Looking among a larger variety of papers on thermal properties of wide-band gap-materials for further potential targets for the applications of the presently developed analytical apparatus, the case of MgO [41,42,43,44] is the most prominent. Indeed, one can see from the corresponding ${\Theta _D}\left( T \right)$ and/or ${\Theta _{Dh}}\left( T \right)$ curves shown in these papers (cf. Figure 2 in [41], Figure 5 in [42], Figure 4 in [43], and Figure 3 in [44]) that, after a rapid fall from a $T \to 0$ limiting level of ${\Theta _D}\left( 0 \right) \cong $ 946 K to a minimum of ${\Theta _{D\;min}} \approx $ 750 K (±10 K), the further run from 100 K up to ~700 K takes place at a nearly constant Debye temperature of ${\Theta _D} \approx $ 760 K. Thus, it should be possible to describe the corresponding ${C_p}\left( T \right)$ data for MgO, from 100 K up to high temperatures, e.g., by the interpolation formula given by Eq. (4.3) parameterized by an effective Debye temperature of ${\Theta _D} \approx $ 760 K. In contrast, the obviously unconsidered choice of the limiting $\left( {T \to 0} \right)$ Debye temperature of ${\Theta _D}\left( 0 \right)$ = 946 K reported in previous studies [45,46] can hardly be considered as an adequate ${\Theta _D}$-value for a proper description of the actual ${C_p}\left( T \right)$ dependence of MgO (from the cryogenic region up to 2100 K).

Clearly, the arbitrary choice of ${\Theta _D} \to $ 920 K for ZnO [45,46] is also inadequate. This is due to the circumstance that the ${\Theta _D}\left( T \right)$ dependence of ZnO shows a very pronounced “snaky” shape (in analogy to the typical ${\Theta _D}\left( T \right)$ behavior of semiconductor materials [2,3,4,5,6,7,8,9,10,11]). Thus, it is impossible to find a certain range in which the Debye temperature would have a chance to be considered as approximately constant. Furthermore, one can readily estimate (via Eqs. (4.1a) and (4.1b)) that the whole range of possible ${\Theta _{Dh}}\left( T \right)$ values for ZnO extends from a lower boundary of ${\Theta _{D\;min}} \approx $ 350 K up to an upper boundary of ${\Theta _{Dh}}\left( \infty \right) \approx $ 690 K. A considerably higher ${\Theta _D}$-value, such as ${\Theta _D} \to $ 920 K, suggested in previous studies [45,46], for ZnO, ranges clearly outside the realm of reason.

One can make a pre-selection of groups of materials, the heat capacities of which should actually have (or $not$ have) a chance to be reasonably described by the presently developed Debye function interpolation formulae, on the basis of their material-specific dispersion coefficients, ${\Delta _P}$ (Eq. (6.2)). On viewing, e.g., a list of the corresponding ${\Delta _P}$ values given for a large variety of wide-band-gap materials (in Table 1 of an earlier study [36]), which had been estimated on the basis of various theoretically calculated PDOS spectra, we see that diamond and c-BN are the only two materials for which their dispersion coefficients, ${\Delta _P} \cong $ 0.275 ±0.005, are closely approaching the Debye’s limiting value of ${\Delta _D} \cong 0.258$. This is the reason of the actual applicability of the presently developed analytical apparatus (see subsections 3.2 and 6.5) to the available ${C_p}\left( T \right)$ datasets for these two materials. In contrast, the ${\Delta _P}$ values (> 0.35) listed in Table 1 (of Ref. [36]) for all the other wide-band-gap-materials under consideration, are throughout > 40% higher than ${\Delta _D}$. In particular, it should be noted that the ${\Delta _P}$ values estimated for GaN and ZnO are even a factor of ~2 higher than ${\Delta _D}$. Such large qualitative differences in the actual PDOS spectra versus Debye’s oversimplified model spectrum are the reason for the strict inapplicability of Debye’s original model (based on an invocation of a fixed Debye temperature, ${\Theta _D}$ = const.).

Analogously, we find for a larger variety of typical semiconductor (Group-IV, Group-III-V, and Group-II-VI) materials considered [37] that their material-specific ${\Delta _P}$ values are throughout > 60% higher than ${\Delta _D}$ (cf. Table I in [37]). Such large differences between Debye’s fictive and the actual degrees of phonon dispersion automatically exclude the eventual applicability of Debye’s original model (invoking a fixed ${\Theta _D}$) to these materials.

On the other hand, by assessing the wealth of theoretically calculated PDOS spectra that are shown for a large variety of insulators [47], one can find at least four binary materials having ${\Delta _P}$ values very similar to Debye’s value, ${\Delta _D} \cong 0.258$. This concerns the two alkaline earth oxides MgO (${\Delta _P} \approx \;$0.29) and CaO (${\Delta _P} \approx \;$0.30) as well as the two alkali halides NaF (${\Delta _P} \approx $ 0.28) and NaCl ${\Delta _P} \approx \;$0.30). Thus, it should be possible to directly apply the presently developed analytical apparatus, if necessary, to the latter wide-band-gap materials.

Finally, we would like to point out that of primary importance for numerical simulations of the ${C_p}\left( T \right)$ datasets available for such low-dispersion materials is the unprecedented interpolation formula of the mixed type, ${\kappa _{D;ND}}\left( T \right)$ (Eq. (6.8)), which smoothly combines the characteristic high-temperature Debye model behavior (Eq. (3.1)) with the qualitatively different limiting low-temperature non-Debye behavior (Eq. (5.1)). The actual use of this ambiguous heat capacity shape function, ${\kappa _{D;ND}}\left( T \right)$, within the frame of the usually considered analytical expression for isobaric heat capacities (Eq. (6.10)) provides, as a rule, good numerical simulations for the comprehensive ${C_p}\left( T \right)$ datasets for any temperature of practical interest, i.e., from the liquid-helium region (including the $T \to 0$ limit) up to the melting points of the low-dispersion materials.

The crucial problem of reasonable theoretical interpretation of the measured ${C_p}\left( T \right)$ datasets for a large variety of wide-band-gap materials, in combination with an adequate analytical description of the qualitatively different Debye temperature behavior, is intended to be treated in due detail in a forthcoming study.

Appendix A Alternative Versions of the Infinite Taylor Series Representations for the High-$T$ Behavior of the Debye Function Integral

A useful analytical tool for deriving alternative versions of the Taylor series representations for the Debye integral represented by Eq. (2.1) was provided by a general method [23] of infinite Taylor series transformations. Within the present context, the starting point for such transformations is given by Debye’s conventional high-$T$ Taylor series expansion, Eq. (2.5), for the Debye integral, Eq. (2.1),

\[ \kappa_{D h}^{\left(m_{D}^{ }\right)}(x)=1+\sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} c_{2 m}^{D} x^{2 m} \tag{A.1} \]

where the dimensionless ($T$-dependent) quantity, $x\left( T \right)$, represents the ratio of the Debye temperature versus the lattice temperature, $x = {\Theta _D}/T$.

A.1 Exponential Series Representation

The first application of this method to the Debye function had been made in a previous study [22] in the form of the ansatz

\[ \log \left[\kappa_{D h}(x)\right]=\log \left[1+\sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} c_{2 m}^{D} \cdot x^{2 m}\right] \rightleftarrows \sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} d_{2 m}^{D} \cdot x^{2 m} \tag{A.2} \]

The actual values of the corresponding expansion coefficients, $d_{2m}^D$, for the $\log \left[ {{\kappa _{Dh}}\left( x \right)} \right]$ function, are obtained by using the known series expansion of the logarithmic function, $\log (1 + z) = z - \frac{{{z^2}}}{2} + \frac{{{z^3}}}{3} - \frac{{{z^4}}}{4} \ldots $, with respect to the compound argument

\[ z=\sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} c_{2 m}^{D} x^{2 m} \tag{A.3} \]

Evaluating the respective expansion coefficients associated with the individual ${x^{2m}}$ power terms (occurring in the argument of the logarithmic function), and identifying them with the novel (transformed) expansion coefficients, $d_{2m}^D$ , we have obtained in previous work [22] the following values for the latter:

\[ d_{2}^{D}=\frac{-1}{20}, d_{4}^{D}=\frac{3}{5600}, d_{6}^{D}=\frac{-17}{2268000}, d_{8}^{D}=\frac{3631}{27941760000}, d_{10}^{D}=\frac{-11521}{4540536000000}. \tag{A.4} \]

(cf. the appendix of [22]). Based on these transformed expansion coefficients (Eq. (A.4)), it was possible to represent (in accordance with Eq. (A.2)) the high-$T$ curve section of the Debye function, ${\kappa _{Dh}}\left( x \right)$ (A.1), in terms of an equivalent (exponential) power series expansion [22],

\[ \kappa_{D h}^{\left(m_{D}^{ }\right)}(x)=\exp \left[\sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} d_{2 m}^{D} x^{2 m}\right] \tag{A.5} \]

Comparing the sequence of transformed expansion coefficients, $d_{2m}^D$ (A.5), with the sequence of Debye’s original ones, $c_{2m}^D$ (Eq. (2.7)), we see that the signs of both sequences of expansion coefficients are alternating, i.e., $d_{2m}^D = {( - 1)^m}\left| {d_{2m}^D} \right|$ and $c_{2m}^D = {( - 1)^m}\left| {c_{2m}^D} \right|$. In addition, we find that the magnitudes of the novel expansion coefficients, $\left| {d_{2m}^D} \right|$, are decreasing more rapidly (with increasing order, $2m$ = 4, 6, 8,…) than the magnitudes, $\left| {c_{2m}^D} \right|$, of Debye’s original high-$T$ expansion coefficients. This important feature is the reason for the markedly improved convergence properties of the novel (exponential) series expansion (A.5) in comparison with Debye’s original one (Eq. (2.5)).

We have shown in Figure 6 of the appendix of an earlier work [22] that, e.g., the deviations of a truncated ${\kappa _{Dh}}\left( x \right)$ expression of this exponential type (Eq. (A.5)),

\[ \kappa_{D h}^{\left(m_{D}^{ }=5\right)}(x)=\exp \left[\sum_{m=1}^{m_{D}^{ }=5} d_{2 m}^{D} x^{2 m}\right] \tag{A.6} \]

from the exact Debye function, ${\kappa _D}\left( x \right)$ (Eq. (2.1)), are smaller than 1% throughout an interval of $x = {\Theta _D}/T$ from 0 to ~5. This advantageous feature of Eq. (A.6) already represented significant progress in comparison with Debye’s original high-$T$ approximation (Eq. (2.5)).

A.2 High-Temperature Representation for Reciprocal Debye Function Values

In the course of the forthcoming analytical studies, we succeeded in establishing a few more (structurally different) versions of the high-$T$ representations for the Debye function (Eq. (2.1)), which turned out to show even markedly better convergence properties than the exponential one sketched above. An interesting example had been presented, firstly, in [11] with respect to the reciprocal values, ${\left[ {{\kappa _{Dh}}\left( x \right)} \right]^{ - 1}}$, of the ${\kappa _{Dh}}\left( x \right)$ function (Eq. (2.5)). Accordingly we have made, in a previous study [11], the alternative ansatz

\[ \left[\kappa_{D h}(x)\right]^{-1}=\left[1+\sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} c_{2 m}^{D} \cdot x^{2 m}\right]^{-1} \rightleftarrows 1+\sum_{m=1}^{m_{R}^{ }(\rightarrow \infty)} c_{2 m}^{D R} \cdot x^{2 m} \tag{A.7} \]

In this case, the actual values of the corresponding expansion coefficients, $c_{2m}^{DR}$, for the ${\left[ {{\kappa _{Dh}}\left( x \right)} \right]^{ - 1}}$ function have been obtained using the familiar series expansion ${\left( {1 + z} \right)^{ - 1}} = 1 - z + {z^2} - {z^3} + ...$ (with respect to the compound argument $z$ (Eq. (A.3)). Evaluating the respective expansion coefficients associated with the individual ${x^{2m}}$ power terms (occurring in the argument of the ${(1 + z)^{ - 1}}$ function, Eq. (A.7)), and identifying them with the novel (transformed) expansion coefficients, $c_{2m}^{DR}$, we have obtained, in [11], the following values for the latter:

\[ c_{2}^{D R}=+\frac{1}{20}, c_{4}^{D R}=+\frac{1}{1400}, c_{6}^{D R}=+\frac{1}{648000}, c_{8}^{D R}=-\frac{73}{3492720000}, \ldots \tag{A.8} \]

(concerning the values of the subsequent four expansion coefficients, up to an order of $2{m_R} = \;16$, see Table 6 in [11]). Consequently, on the basis of the latter, it was possible to represent (in accordance with Eq. (A.7)) the high-$T$ curve section, ${\kappa _{Dh}}\left( x \right)$, of the Debye function (2.1)), in terms of the reciprocal values of a corresponding Taylor series expansion,

\[ \kappa_{D h}^{\left(m_{R}^{ }\right)}(x)=\left[1+\sum_{m=1}^{m_{R}^{ }(\rightarrow \infty)} c_{2 m}^{D R} x^{2 m}\right]^{-1} \tag{A.9} \]

(in agreement with Eq. (A.15) given in [11]).

It has been shown in Figure 8 of Ref. [11] that the deviations, e.g., of a truncated ${\kappa }_{Dh}^{\left(m_R^{ }= 5\right)}\left(x\right)$ expansion of the type given by Eq. (A.9) from the exact Debye function, ${\kappa _D}\left( x \right)$ (Eq. (2.1)), are smaller than 1% throughout the interval of $x = {\Theta _D}/T$ values from 0 to ~7. This widening of the range of the possible applicability of Eq. (A.9) already represented significant progress over the somewhat worse convergence properties of the preceding exponential series representation (Eq. (A.5)). This improvement of convergence properties is because the signs of the first three expansion coefficients (A.8) are the same, $c_{2m = \;2,4,6}^{DR} > 0$ (in contrast to the alternating signs of the expansion coefficients for the exponential series representation, $d_{2m}^D = {( - 1)^m}\left| {d_{2m}^D} \right|$).

It should be noted that the ${\kappa }_{Dh}^{\left(m_R^{ }= 8\right)}\left(x\right)$ version of Eq. (A.9) had been used in [11] for performing high-precision calculations (accurate to 10 decimal places, within the interval $0 \le {\Theta _D}/T \le 1.6$) in combination with the usage of the conventional low-$T$ formula ${\kappa }_{Dl}^{\left(n_D^{}= 15\right)}\left(x\right)$ (Eq. (2.2)) for the interval $1.6 \le {\Theta _D}/T \le 40$ (see the corresponding results for selected ${\kappa _D}\left( x \right)$ values listed in Table 7 of Ref. [11]).

A.3 Unprecedented Series Representation for Reciprocal Square-root Debye Function Values

Another crucial example of a rapidly converging Taylor series expansion which turned out to play an especially useful role within the present context, concerns the square of reciprocal values, ${(1/{\kappa _{Dh}}\left( x \right))^2}$, of the ${\kappa _{Dh}}\left( x \right)$ function (A.1).

Accordingly, we make the alternative ansatz

\[ \left(\frac{1}{\kappa_{D h}^{\left(m_{D}^{ }\right)}(x)}\right)^{2}=\left[1+\sum_{m=1}^{m_{D}^{ }(\rightarrow \infty)} c_{2 m}^{D} \cdot x^{2 m}\right]^{-2} \rightleftarrows 1+\sum_{m=1}^{m_{P}^{ }(\rightarrow \infty)} r_{2 m}^{D} \cdot x^{2 m} \tag{A.10} \]

In this case, the actual values of the corresponding expansion coefficients, $r_{2m}^D$, for the ${\left[{\kappa }_{Dh}^{\left(m_D^{}\right)}\left(x\right)\right]}^{-2}$ function have been obtained via the usage of the familiar Taylor series expansion, ${\left( {1 + z} \right)^{ - 2}} = \;1 - 2z + 3{z^2} - 4{z^3} + \ldots $ (with respect to the compound argument $z$ (Eq. (A.3)). Evaluating the respective expansion coefficients associated with the individual ${x^{2m}}$ power terms (occurring in the argument of the ${(1 + z)^{ - 2}}$ function, Eq. (A.10)), and identifying them with the novel (transformed) expansion coefficients, $r_{2m}^D$, we have obtained the following values for the latter:

\[ r_{2}^{D}=\frac{1}{10}, r_{4}^{D}=\frac{11}{2800}, r_{6}^{D}=\frac{169}{2268000}, r_{8}^{D}=\frac{29}{46569600} \tag{A.11} \]

(The values of the subsequent four expansion coefficients, up to an order of $2{m_P^{ }} = \;16$, are given in Table 1).

Consequently, on the basis of the latter coefficients (A.11), it is possible to represent (in accordance with Eq. (A.10)) the high-$T$ curve section, ${\kappa _{Dh}}\left( x \right)$, of the Debye function (2.1)), in the form of the following algebraic expression:

\[ \kappa_{D h}^{\left(m_{P}^{ }\right)}(x)=\left[1+\sum_{m=1}^{m_{P}^{ }(\rightarrow \infty)} r_{2 m}^{D} x^{2 m}\right]^{-\frac{1}{2}} \tag{A.12} \]

The advantageous properties of the latter representation have been visualized and discussed in detail in Section 3.

Here, we have used the ${\kappa }_{Dh}^{\left(m_P^{}= 8\right)}\left(x\right)$ version of Eq. (A.12) for performing high-precision calculations (accurate to 13 decimal places) within the interval $0 \le {\Theta _D}/T \le 1.6$ (as indicated in the lower inset of Figure 1), in combination with the usage of the conventional low-$T$ formula ${\kappa }_{Dl}^{\left(n_D^{}= 20\right)}\left(x\right)$ (Eq. (2.2)), for the interval $1.6 \le {\Theta _D}/T \le 40$. The corresponding high-accuracy results for selected couples of $x$ and ${\kappa _D}\left( x \right)$ values are listed in the first and second columns of Table A.1.

Table A.1 Exemplary high-accuracy Debye function values, ${\kappa _D}\left( x \right)$, for selected magnitudes of the ratios ${\Theta _D}/T = x$. The respective calculations have been performed due to a combination of the unprecedented algebraic ${\kappa }_{Dh}^{\left(m_P^{}= 8\right)}\left(x\right)$ expression (Eq. (A.12)), for the interval $0 \le {\Theta _D}/T \le 1.6$, with Debye’s conventional ${\kappa }_{Dl}^{\left(n_D^{}= 20\right)}\left(x\right)$ expression (Eq. (2.2)), for the interval $1.6 \le {\Theta _D}/T \le 45$.

Finally, we would like to mention that the presently performed Taylor series transformation indicated in Eq. (A.10), as well as the preceding ones [11,22] (indicated in Eq. (A.7) and (A.2)), could be readily performed (in both directions) within the framework of a representative mathematical standard program collection like MAPLE.

Author Contributions

The author did all the research work of this study.

Competing Interests

The author has declared that no competing interests exist.