Efficient Debye Function Interpolation Formulae: Sample Applications to Diamond
Technische Universität Chemnitz, Institut für Physik, D09107 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.
Abstract
The wellknown classical heat capacity model developed by Debye proposed an approximate description of the temperaturedependence of heat capacities of solids in terms of a characteristic integral, the Tdependent values of which are parameterized by the Debye temperature, $\Theta _D$. However, numerous tests of this simple model have shown that within Debye’s original supposition of approximately constant, materialspecific Debye temperature, it has little chance to be applicable to a larger variety of nonmetals, except for a few widebandgap materials such as diamond or cubic boron nitride, which are characterized by an unusually low degree of phonon dispersion. In this study, we present a variety of structurally simple, unprecedented algebraic expressions for the hightemperature behavior of Debye’s conventional heatcapacity integral, which provide fine numerical descriptions of the isochoric (harmonic) heat capacity dependences parameterized by a fixed Debye temperature. The present sample application of an appropriate hightolow temperature interpolation formula to the isobaric heat capacity data for diamond measured by Desnoyers and Morrison [17], Victor [24], and Dinsdale [25] provided a fine numerical simulation of data within a range of 200 to 600 K, involving a fixed Debye temperature of about 1855 K. Representing the monotonically increasing difference of the isobaric versus isochoric heat capacities by two associated anharmonicity coefficients, we were able to extend the accurate fit of the given heat capacity ($C_p \left( T \right)$) data up to 5000 K. Furthermore, we have performed a highaccuracy fit of the whole $C_p \left( T \right)$ dataset, from approximately 20 K to 5000 K, on the basis of a previously developed hybrid model, which is based on two continuous lowT curve sections in combination with three discrete (Einstein) phonon energy peaks. The two theoretical alternative curves for the $C_p \left( T \right)$ dependence of diamond were found to be almost indistinguishable throughout the interval from 200 K to 5000 K.
Keywords
Heat capacities; Debye function; nonDebye behavior; widebandgap materials; diamond; phonon dispersion; Debye temperature
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 phonondensitiesofstates (PDOS) spectra of metals as well as nonmetals, 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 cutoff (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 modelspecific parameter, the socalled 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 nonmetals (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 IIIV materials [4,5,6] and IIVI materials [7,8,9,10,11].
These materialspecific ${\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) DulongPetit 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, materialspecific Debye temperatures, ${\Theta _D}$.
On the other hand, among a larger variety of widebandgap 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 quasiplateau 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 (cBN) above 150 K (see Figure 2 in [19]). These two examples of widegap 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 leastmeansquare fits of the duly limited (partial) ${C_{V/p}}(T > {\Theta _D}\left( 0 \right)/12)$ datasets (similar to those available for diamond or cBN) 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 lowand hightemperature approximations [11] for the Debye function integral, which are well known already from Debye’s original paper [1]. We confirm that Debye’s lowtemperature 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 hightemperature 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 evenordered ${\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 (nonlinear) functions. In the appendix, we briefly present two previous sample applications of this method to the hightemperature 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 highcapacity Debye function formula, ${\kappa}_{Dh}^{\left({m}_{P}^{}\right)}\mathrm{\left(x\right)}$ (Eq. (A.12)), with various truncated versions of Debye’s conventional lowcapacity 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 leastmeansquare 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 lowtohigh temperature interpolation formula in combination with additional anharmonicityrelated terms, which come strongly into play in the hightemperature 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 highaccuracy 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 hightolow$T$ interpolation formulae (in Section 4), on the one hand, and from the highaccuracy 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 lowtohightemperature heat capacity datasets, which are available for certain widebandgap 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 Lowand Hightemperature 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 DulongPetit 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 materialspecific) 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 lowtemperature 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 lowtemperature 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}}{5{x}^{3}}\frac{3x}{\mathrm{exp}(x)1}12x\sum _{n=1}^{{n}_{D}^{}(\to \mathrm{\infty})}\mathrm{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 highaccuracy 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 hightemperature 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 m1)\leftB_{2 m}\right z^{2 m2}}{(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 hightemperature 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 m1)\leftB_{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 hightemperature 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 hightemperature representations of the Debye integral.
3. Alternative Hightemperature Representations for the Debye Function
A very useful analytical tool for constructing alternative versions of the Taylor series representations for the hightemperature 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 (nonlinear) 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 lowestorder expansion coefficients, $r_{2 \le 2m \le 16}^D$, implied by the unprecedented highcapacity (hightemperature) algebraic formula, ${\kappa}_{Dh}^{({m}_{P}^{}=8)}\mathrm{\left(x\right)}$ (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 hightemperature Taylor series expansion, ${\kappa}_{Dh}^{({m}_{D}^{}=8)}\mathrm{\left(x\right)}$ (Eq. (25)).
However, with respect to the present purpose of simulating experimentally measured heat capacities of a lowdispersion widebandgap 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 hightemperature algebraic representation, ${\kappa}_{Dh}^{\left({m}_{P}^{}\right)}\mathrm{\left(x\right)}$ 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 ${\kappa}_{Dh}^{({m}_{P}^{}=4)}\mathrm{\left(x\right)}$ curve from the exact ${\kappa _D}\left( x \right)$ curve remain limited to less than 1% throughout an unusually large hightemperature 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 ${\stackrel{~}{\kappa}}_{D}^{\left({n}_{D}^{}\ge 1;{m}_{P}^{}=4\right)}\left(x\right)$ (Eq. (3.7)) involving truncated versions of Debye’s conventional lowcapacity representation, ${\kappa}_{Dl}^{\left({n}_{D}^{}\right)}\mathrm{\left(x\right)}$ (Eq. (2.2)), of the two lowest orders, ${n_D^{ }}$= 1 or 2, respectively.
3.1 Alternate Usage of Highand Lowtemperature 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 lowtemperature 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 hightemperature 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 lowand hightemperature 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 hightemperature curve, ${\kappa}_{Dh}^{\left({m}_{P}^{}=7\right)}(x7)$ (Eq. (3.1);curve shown in the middle inset of Figure 1), with various (truncated) samples of Debye’s conventional lowcapacity (lowtemperature) 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 highand lowtemperature 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) highand lowtemperature 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^{−7}. 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 highand lowtemperature 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 lowtemperature 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 hightemperature 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 hightemperature 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 lowtemperature 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 noncrossing highand lowtemperature 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 formulaspecific 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 leastmeansquare 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 lowtemperature expression, ${\kappa _{Dl}}\left( x \right)$ (Eq. (2.2)), into a conveniently constructed interpolation formula. Accordingly, we consider a combination of a truncated lowtemperature 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 hightemperature 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 (orderofmagnitude) 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 hightemperature expression, ${\kappa}_{Dh}^{\left({m}_{P}^{}=4\right)}\left(x\right)$ (Eq. (2.1)), with lowtemperature 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 lowestorder versions, where Debye’s lowtemperature 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 orderofmagnitude 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 highaccuracy approximation for the Debye function integral provided that at least two summation steps are retained in the corresponding lowcapacity 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 higherorder 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 higherorder 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 highand lowtemperature 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 lowtohigh 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 widebandgap 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 highaccuracy formulae (Eq. (4.1a) and (4.1b)) from the (isobaric) heat capacities, ${C_p}\left( T \right)$, measured in the lowcapacity region by Desnoyer and Morrison [17] (○) and in the hightemperature 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 lowtemperature 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 wellknown lowtemperature ${\Theta _D}\left( T \right)$ data, the different sets of hightemperature ${\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 pointbypoint 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 lowtemperature 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 lowtohigh 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 nonnegligible 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 materialspecific 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 leastmeansquare 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^{−5} K^{−1} and ${A_2}$= 2.421 × 10^{−9} K^{−2}. 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 dashdotted 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 (intermediatetohightemperature) 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 hightemperature (highcapacity) algebraic formula ${\kappa}_{Dh}^{({m}_{P}^{}=4)}\mathrm{\left(x\right)}$ (Eq. (A.12)), in combination with Eq. (4.4) (for the region 230 K < $T$ < 5000 K, at least). An exemplary highaccuracy 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 hightemperature 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 hightemperature (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 dashdoubledot curve (······) in the inset of Figure 3.
5. Alternative Highaccuracy Fit due to the Representative Hybrid Model
The actual variation of the Debye temperature, ${\Theta _D}\left( T \right)$, of the widebandgap as well as semiconductor materials is wellknown 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 IIIV semiconductors in [4,5,6], and for IIVI 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 (liquidhydrogen) 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 oddorder 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 GroupIV, IIIV, and IIVI materials in [12]). However, this truncated oddorder 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 nonnegligible contributions of an electronic system to the heat capacity in the liquidheliumregion, 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 wellfulfilled 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 NonDebye 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 wellknown 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 lowestorder power terms occurring in the wellknown empirical lowtemperature 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 wellknown 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 codetermined 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 highaccuracy 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 multiparameter leastmeansquare 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 hightemperature 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 dashdoubledot 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 hightolowtemperature 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 (dashdot 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 wellknown nonDebye 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 Highaccuracy 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 leastmeansquare 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 modelspecific 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^{−4 }mol^{−1} (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^{−4 }mol^{−1}, 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 wellknown $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 byproducts 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 firstand secondorder 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)^{2}. 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 secondorder 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 materialspecific magnitude of the dispersion coefficient, ${\Delta _P}$ (Eq. (6.1)), for diamond is only approximately 7% higher than its modelspecific 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 hightemperature 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 wellknown 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\text{}\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 highaccuracy 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 dashdot 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 thermochemistry, 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 lowtemperature plateau behavior of the respective, approximate ${C_p}\left( T \right)/{T^3}$ function (dashdotcurve), 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 hightemperature 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 dashdoubledot 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 Subquartic Behavior of the Lowtemperature 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 widebandgap 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(TT_{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^{−4 }mol^{−1} and ${d_1} = {\rho ^{\left( 1 \right)}}\left( {{T_i}} \right)$ = 0.13911×10^{−2} µJ K^{−5 }mol^{−1} (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 subquartic 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 MultiDebyefunction 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 lowtemperature regions, with the popular (thermochemical) MaierKelley 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 NonDebye Interpolation Formula
Satisfactory numerical fitting to the comprehensive ${C_p}\left( T \right)$ datasets for the highdispersion widebandgap 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 nonDebye formula. A suitable alternative version of the latter had been displayed and used in another work [6] for comprehensive fitting (from liquidhelium 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 IIIV (highdispersion [37]) semiconductor materials.
We have found that an appropriate specialization of the preceding analytical nonDebye model [6] with respect to its possible applications to lowdispersion 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 nonDebye 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 oddorder 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 hightemperature 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) Debyefunctionrelated predecessor, $r_8^D$ (Table 1). This generalization is inevitable for assuring a smooth change between the respective asymptotic lowand hightemperature 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^{−4 }mol^{−1}, ${c_5}$ = 0.32631×10^{−5 }µJ K^{−6 }mol^{−1}, and ${c_7}$= 0.40863×10^{−9 }µJ K^{−8 }mol^{−1} (cf. [12]) in combination with a fitted nonDebye expansion coefficient of $r_8^{ND}$= 1.862 × 10^{−6}, an effective Debye temperature of ${\Theta _D}$= 1904.8 K, and the associated anharmonicity coefficients of ${A_1}$= 1.799×10^{−5} K^{−1} and ${A_2}$= 3.278 × 10^{−9} K^{−2}. In particular, assessing the presently obtained ${\Theta _D}$ value, we see that it is by only approximately 0.6% higher than the theoretical hightemperature 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 byproduct of fitting, even the correct orderofmagnitude of the materialspecific ${\Theta _{Dh}}\left( {T \to \infty } \right)$ values. Furthermore, we see that the fitted magnitude of the NonDebye 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 nonDebye 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^{−4 }mol^{−1} 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^{−4 }mol^{−1}, (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 lowtemperature 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 nonDebye 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 widebandgap material under consideration is not significantly higher than those estimated for diamond and cBN [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 lowtohigh temperature heat capacity datasets of certain widebandgap 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 lowdispersion materials is the quasiplateaubehavior 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 cBN).
Looking among a larger variety of papers on thermal properties of wideband gapmaterials 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 preselection 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 materialspecific 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 widebandgap 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 cBN 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 widebandgapmaterials 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 (GroupIV, GroupIIIV, and GroupIIVI) materials considered [37] that their materialspecific ${\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 widebandgap 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 lowdispersion materials is the unprecedented interpolation formula of the mixed type, ${\kappa _{D;ND}}\left( T \right)$ (Eq. (6.8)), which smoothly combines the characteristic hightemperature Debye model behavior (Eq. (3.1)) with the qualitatively different limiting lowtemperature nonDebye 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 liquidhelium region (including the $T \to 0$ limit) up to the melting points of the lowdispersion materials.
The crucial problem of reasonable theoretical interpretation of the measured ${C_p}\left( T \right)$ datasets for a large variety of widebandgap 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 HighTemperature 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 highprecision 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 Squareroot 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 highprecision 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 highaccuracy 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 highaccuracy 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.
References
 Debye P. The theory of specific warmth. Ann Phys. 1912; 39: 789839. [CrossRef]
 Flubacher P, Leadbetter AJ, Morrison JA. The heat capacity of pure Silicon and Germanium and properties of their vibrational frequency spectra. Philos Mag. 1959; 4: 273294. [CrossRef]
 Pässler R. Dispersionrelated theory for heat capacities of semiconductors. Phys Stat Sol B. 2007; 244: 46054623. [CrossRef]
 Piesbergen U. Heat capacity and Debye temperatures. In: Semiconductors and semimetals. New York: Academic Press; 1966. [CrossRef]
 Blakemore JS. Semiconducting and other mayor properties of gallium arsenide. J Appl Phys. 1982; 53: R123. [CrossRef]
 Pässler R. NonDebye heat capacity formula refined and applied to GaP, GaAs, GaSb, InP, InAs, and InSb. AIP Adv. 2013; 3. doi: 10.1063/1.4818273. [CrossRef]
 Rusakov AP, Vekilov YK, Kadyshevich AE. Specific heats of CdTe and HgTe and their vibrational frequency spectra. Fiz Tverd Tela. 1970; 12: 32383243.
 Vegalatos N, Wehe D, King JS. Phonon dispersion and phonon densities of states for ZnS and ZnTe. J Chem Phys. 1974; 60: 3613. [CrossRef]
 Talwar DN, Vandevyver M, Kunc K, Zigone M. Lattice dynamics of zinc chalcogenides under compression: Phonon dispersion, mode Grüneisen, and thermal expansion. Phys Rev B. 1981; 24: 741. [CrossRef]
 Pässler R. NonDebye behaviours of heat capacities of cubic II–VI materials. J Phys Chem Solids. 2011; 72: 12961311. [CrossRef]
 Pässler R. Unprecedented integralfree debye temperature formulas: Sample applications to heat capacities of ZnSe and ZnTe. Adv Condens Matter Phys. 2017; 2017. doi: 10.1155/2017/9321439. [CrossRef]
 Pässler R. Limiting Debye temperature behavior following from cryogenic heat capacity data for groupIV, III–V, and II–VI materials. Phys Stat Sol B. 2010; 247: 7792. [CrossRef]
 Berg WT, Morrison JA. The thermal properties of alkali halide crystals. I. The heat capacity of potassium chloride, potassium bromide, potassium iodide and sodium iodide between 2.8 and 270˚K. Proc Roy Soc A. 1957; 242: 467477. [CrossRef]
 Barron TH, Berg WT, Morrison JA. The thermal properties of alkali halide crystals II. Analysis of experimental results. Proc Roy Soc A. 1957; 242: 478492. [CrossRef]
 Tosi MP, Fumi FG. Temperature dependence of the Debye temperatures for the thermodynamic functions of alkali halide crystals. Phys Rev. 1963; 131: 14581465. [CrossRef]
 Gopal ES. Specific heats at low temperatures. New York: Plenum Press; 1966. [CrossRef]
 Desnoyers JE, Morrison JA. The heat capacity of diamond between 12·8° and 277°k. Phil Mag. 1958; 3: 4248. [CrossRef]
 Dolling G, Cowley RA. The thermodynamic and optical properties of germanium, silicon, diamond and gallium arsenide. Proc Roy Soc. 1966; 88: 463. [CrossRef]
 Inaba A,Yoshiasa A. Lowtemperature heat capacity of wurtzitetype boron nitride. Jpn J Appl Phys. 1997; 36: 5644. [CrossRef]
 Schrödinger E. Spezifische wärme. In: Thermische eigenschaften der stoffe. Berlin: SpringerVerlag; 1926. p. 275320. [CrossRef]
 Blackman M. The specific heat of solids. In: Encyclopedia of physics. Berlin, Heidelberg: SpringerVerlag; 1955. p. 325385. [CrossRef]
 Pässler R. Exponential series representation for heat capacities of semiconductors and widebandgap materials. Phys Stat Sol B. 2008; 245: 11331146. [CrossRef]
 Schell HJ. Unendliche Reihen. Leipzig: BSB Teubner; 1974. [CrossRef]
 Victor AC. Heat capacity of diamond at high temperatures. J Chem Phys. 1962; 36: 1903. [CrossRef]
 Dinsdale AT. SGTE data for pure elements. Calphad. 1991; 15: 317425. [CrossRef]
 Pässler R. Representative hybrid model used for analyses of heat capacities of groupIV, III–V, and II–VI materials. Phys Stat Sol B. 2011; 248: 904920. [CrossRef]
 Barron TH, Morrison JA. On the specific heat of solids at low temperatures. Can J Phys. 1957; 35: 799810. [CrossRef]
 Cetas TC, Tilford CR, Swenson CA. Specific heats of Cu, GaAs, GaSb, InAs, and InSb from 1 to 30°K. Phys Rev. 1968; 174: 835844. [CrossRef]
 Holste JC. Specific heats of GaSb, GaAs, InSb, InAs, Bi, Cd, Sn, and Zn below 30 K. Phys Rev B. 1972; 6: 24952497. [CrossRef]
 Berg WT. Lowtemperature heat capacities of silver chloride and lithium iodide. Phys Rev B. 1976; 13: 26412645. [CrossRef]
 Barron TH, White GK. Heat capacity and thermal expansion at low temperatures. Boston, MA: Springer; 1999. [CrossRef]
 Pässler R. Characteristic nonDebye heat capacity formula applied to GaN and ZnO. J Appl Phys. 2011; 110. doi: 10.1063/1.3622668. [CrossRef]
 Siethoff H, Ahlborn K. The dependence of the Debye temperature on the elastic constants. Phys Stat Sol B. 1995; 190: 179191. [CrossRef]
 Adachi S. Properties of groupIV, IIIV, and IIVI semiconductors. Hoboken, New Jersey: John Wiley & Sons, Ltd; 2005. [CrossRef]
 LandoltBörnstein: Numerical data and functional relationships in science and technology. Berlin, Germany: SpringerVerla; 1982, 1989, and 19992002.
 Pässler R. Moments of phonon density of states spectra and characteristic phonon temperatures of wide band gap materials. Phys Stat Sol B. 2006; 243: 27192727. [CrossRef]
 Pässler R. Basic moments of phonon density of states spectra and characteristic phonon temperatures of group IV, III–V, and II–VI materials. J Appl Phys. 2007; 101. doi: 10.1063/1.2721749. [CrossRef]
 Pässler R. Qualitatively different nonquartic lowtemperature heat capacity behaviors due to various nonmetallic solids. Phys Stat Sol. 2019; 256. doi: 10.1002/pssb.201900154. [CrossRef]
 Vassiliev VP, Taldrik AF. Description of the heat capacity of solid phases by a multiparameter family of functions. J Alloys Compd. 2021; 872: 159682. [CrossRef]
 Maier CG, Kelley KK. An equation for the representation of hightemperature heat content data. J Am Chem Soc. 1932; 54: 32433246. [CrossRef]
 Barron TH, Berg WT, Morrison JA. On the heat capacity of crystalline magnesium oxide. Proc Roy Soc. 1959; 250: 7083. [CrossRef]
 Gmelin E. Thermal properties of alcalineearthoxides I. Specific heat measurements. Z Naturforsch.1969; 24: 17941800. [CrossRef]
 Sangster MJ, Peckham G, Saunderson DH. Lattice dynamics of magnesium oxide. J Phys C Solid St Phys. 1970; 3: 1026. [CrossRef]
 Singh RK, Upadhyaya KS. Crystal dynamics of magnesium oxide. Phys Rev B. 1972; 6: 15891596. [CrossRef]
 Mamedov BA, Eser E, Koç H, Askerov IM. Accurate evaluation of the specific heat capacity of solids and its application to MgO and ZnO crystals. Int J Thermophys. 2009; 30: 10481054. [CrossRef]
 Bouafia M, Sadat H. An approximate model for the heat capacity of solids. Int J Thermophys. 2020; 41. doi: 10.1007/s1076502002658z. [CrossRef]
 Bilz H, Kress W. Phonon dispersion relations in insulators. New York: Springer; 1979. [CrossRef]