Recent Progress in Materials Efficient Debye Function Interpolation Formulae: Sample Applications to Diamond

The well-known classical heat capacity model developed by Debye proposed an approximate description of the temperature-dependence of heat capacities of solids in terms of a characteristic integral, the T -dependent values of which are parameterized by the Debye temperature, 𝛩 𝐷 . However, numerous tests of this simple model have shown that within Debye’s original supposition of approximately constant , material-specific Debye temperature, it has little chance to be applicable to a larger variety of non-metals, except for a few wide -band-gap 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 high- temperature behavior of Debye’s conventional heat-capacity 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 high-to-low temperature interpolation formula to the isobaric heat capacity data for diamond measured by Desnoyers and Morrison [17], Victor [24], and Dinsdale two associated anharmonicity coefficients, we were able to extend the accurate fit of the given heat capacity ( 𝐶 𝑝 (𝑇) ) data up to 5000 K. Furthermore, we have performed a high-accuracy fit of the whole 𝐶 𝑝 (𝑇) dataset, from approximately 20 K to 5000 K, on the basis of a previously developed hybrid model, which is based on two continuous low- T curve sections in combination with three discrete (Einstein) phonon energy peaks. The two theoretical alternative curves for the 𝐶 𝑝 (𝑇) dependence of diamond were found to be almost indistinguishable throughout the interval from 200 K to 5000 K.


Introduction
A basic model for approximate numerical simulations of measured heat capacities had been constructed more than 100 years ago by Debye [1]. This familiar model was based on the assumption that the relevant phonon-densities-of-states (PDOS) spectra of metals as well as non-metals, should be describable to a good approximation in the form of purely quadratic phonon energy dependences, ( ) ∝ 2 , from 0 up to a certain cut-off (Debye) energy, . 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 ( ( )), which are generally defined as being given by the partial derivatives, ( ) ≡ ( ( , )/ ) , of internal energies, ( , ), 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 ( ) in terms of certain Debye function integrals [1], the T-dependent magnitudes of which are governed by a unique model-specific parameter, the so-called Debye temperature, = / (where represents the Boltzmann constant).
However, since the middle of the past century, a wealth of experimental studies on the thermal properties of non-metals (semiconductors as well as isolators) has been performed. From the numerous results of these studies, one could conclude that, as a rule, it becomes possible to simulate the measured heat capacity data using the Debye function integrals only under the assumption that the Debye temperature depends (more or less strongly) on the lattice temperature, T. The corresponding ( ) dependencies have been found to be very pronounced for typical semiconductor materials, in particular, for Si and Ge [2,3] as well as for numerous III-V materials [4][5][6] and II-VI materials [7][8][9][10][11].
These material-specific ( ) dependences show "snaky" shapes, which exhibit the following: (i) Rapid fall from their → 0 limiting magnitudes [12], (0), to significantly lower minima, , which are usually located at temperatures of the order (0)/12 (ii) Subsequent rise up to a certain maxima, , which are frequently located at temperatures of the order (0)/2 (iii) Final drop to zero, ( ) = 0, at the characteristic points on the T-scale, , where the measured ( ) curves are crossing the classical (limiting) Dulong-Petit value for harmonic lattice heat capacities, ( ) = ℎ (∞) = 3 (where represents the gas constant and 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" ( ) curves are significantly higher than the associated minima, > . 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, / ( ), in terms of the Debye function integrals governed by fixed, material-specific Debye temperatures, .
On the other hand, among a larger variety of wide-band-gap materials ( > 3 eV), one can find at least a few cases where the differences between and are relatively small (being limited to a few %). This indicates that the ( ) values within the respective T-intervals, from approximately (0)/12 up to approximately (0)/2, are more or less constant. Such a limited quasi-plateau behavior, ( ) ≈ = ., of a ( ) 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 ( ) curve has also been observed for cubic carbon nitride (c-BN) above 150 K (see Figure 2 in [19]). These two examples of wide-gap materials might already be sufficient to justify the present aim of deriving efficient formulae for the conventional Debye function integral (governed by fixed ), which should hence enable us to perform careful least-mean-square fits of the duly limited (partial) / ( > (0)/12) datasets (similar to those available for diamond or c-BN) without involving numerical calculation procedures for the original Debye function integrals.
In Section 2, we give a brief overview of the derivation of the conventional low-and hightemperature approximations [11] for the Debye function integral, which are well known already from Debye's original paper [1]. We confirm that Debye's low-temperature formula [1,11,20,21] is in principle capable of providing quite adequate numerical values for the Debye function integral even at lattice temperatures comparable with the Debye temperature, provided the respective summation procedure is not truncated at a very low order of summation terms. In contrast to the latter, we find again (in accordance with [1,3,11,22]) that Debye's original high-temperature Taylor series expansion is of little practical use in view of its rather bad convergence properties, which are caused by the alteration of signs of the respective expansion coefficients, i. e. 2 = (−1) | 2 | for the subsequent even-ordered ( ) 2 terms ( = 1, 2, 3, …).
Significant improvements of convergence properties have already been shown in previous studies [11,22] to be realizable by means of suitable transformations [23] of the given Taylor series expansions into structurally different versions that are embedded into conveniently chosen alternative (non-linear) functions. In the appendix, we briefly present two previous sample applications of this method to the high-temperature behavior of the Debye function, which resulted in the construction of an exponential series representation [22], and a derivation of an alternative Taylor series representation for reciprocal Debye function values, ( ℎ ( )) −1 [11].
Of special usefulness within the present context is an unprecedented Taylor series representation for the squares of the reciprocal Debye function values, ( ℎ ( )) −2 , which has been presented in the appendix. In Section 3, we qualitatively consider different combinations of truncated versions of the novel high-capacity Debye function formula, ℎ ( ) ( ) (Eq. (A.12)), with various truncated versions of Debye's conventional low-capacity Debye function expansion, ( ) ( )(Eq. (2.2)).
In Section 4, we visualize the actual temperature dependence of the effective (caloric) Debye temperature for diamond, ( ), which is implied by the isobaric heat capacities, ( ), given by Desnoyers and Morrison [17], Victor [24], and Dinsdale [25]. In addition, a least-mean-square fitting has been performed for the whole combination of the respective ( ) datasets, from approximately 200 K up to 5000 K. This fit is based on an appropriate low-to-high temperature interpolation formula in combination with additional anharmonicity-related terms, which come strongly into play in the high-temperature region. As the main result of this fit, we found the effective Debye temperature to be ≅ 1855 K. In Section 5, we give a brief sketch of the basic equations of our previously developed representative hybrid model [11,26], which consists of two continuous low-T curve sections in combination with three discrete (Einstein) phonon energy peaks. We have also performed the corresponding high-accuracy fit of the whole ( ) dataset, from approximately 20 K up to 5000 K.
The two theoretical alternative curves for the ( ) dependence of diamond, resulting from the fit via the unprecedented high-to-low-T interpolation formulae (in Section 4), on the one hand, and from the high-accuracy fit via the previously developed hybrid model [11,26] (in Section 5), on the other hand, have been found to be almost indistinguishable throughout the interval from 200 K to 5000 K. This good agreement illustrates the potential usefulness of the computational framework of approximate Debye function formulae (displayed in Sec. 3) developed in this work and its possible application to the limited low-to-high-temperature heat capacity datasets, which are available for certain wide-band-gap materials characterized by sufficiently low degrees of phonon dispersion.
Various other aspects and additional quantitative information resulting from the comparable alternative fittings of the ( ) datasets, for diamond, under study, by using alternative models have been discussed in Section 6.

Conventional Low-and High-temperature Expressions for the Debye Function
According to Debye's original heat capacity model [1], the isochoric lattice heat capacity is assumed to be given by the product, ( ) = ℎ (∞) ( ) of the familiar Dulong-Petit limit, , of the harmonic lattice heat capacity (for a material containing atoms per molecule) and a normalized (to unity) heat capacity shape function, 0 < ≤ 1. The latter is well known to be defined by an integral of the form [1,4,5,[10][11][12]21] where the ratio between the (presumably material-specific) Debye temperature, , and the lattice temperature, , is denoted as usual by = / . Yet, with respect to the eventual numerical applications of this conventional model to the interpretation of certain experimental data (with respect to which this simple model happens to be applicable), it is desirable to avoid the numerical integrations (in Eq. (2.1)) by using the appropriate analytical or algebraic expressions for more or less extended regions of low and/or high temperatures.
For low-temperature regions ( = / >> 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, ( − 1) −1 = − (1 − − ) −1 = − + −2 + −3 + ⋯. By performing the corresponding integrations (subsequently, by parts), within the second integral of Eq. (2.1), for the low-temperature behavior of the Debye function, ( ), an infinite series expression of the familiar form [1,11,12,20,21] is obtained: 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 → ∞ limit, at least). It may thus be used, if necessary, for performing high-accuracy calculations of the ( ) function (Eq. (2.1)) even for regions of relatively high temperatures ( = / ≈ 1; see below) provided that the summation procedure is extended to sufficiently high order (e. g., up to ≈ 30). At the same time, it is obvious that this infinite series representation (Eq. (2.2)) is not applicable to the → ∞ limit.

Alternative High-temperature Representations for the Debye Function
A very useful analytical tool for constructing alternative versions of the Taylor series representations for the high-temperature behavior, ℎ ( ) (Eq. (2.5)), of the Debye function integral, ( ) (Eq. (2.1)), has been found to be provided by the method of performing suitable transformations [23] of a given Taylor series expansion into structurally different counterparts that are embedded into conveniently chosen alternative (non-linear) functions. This method had been used for the first time in [22] for deriving an equivalent Taylor series for the logarithm, ( ℎ ( )) (Eq. (A.2)), of the conventional ℎ ( ) representation (Eq. (2.5)). Inverting this relationship and observing that ( ( ℎ ( ))) = ℎ ( ), in [22], we came to a corresponding exponential series representation (Eq. (A.5)) for the ℎ ( ) 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/ ℎ ( ) (Eq. (A.7/8)). Inverting the latter relationship, so that [1/ ℎ ( )] −1 = ℎ ( ), we came to an alternative analytical expression (Eq. (A.9)) for the ℎ ( ) function which was found to show still significantly better convergence properties than the preceding exponential representation (for more details, see [11] and subsection A.2 of the appendix). Table 1 The eight lowest-order expansion coefficients, 2≤2 ≤16 , implied by the unprecedented high-capacity (high-temperature) algebraic formula, ℎ  However, with respect to the present purpose of simulating experimentally measured heat capacities of a low-dispersion wide-band-gap material such as diamond (see Section 4), it is useful to display still another version of a rapidly converging series expansion for the ℎ ( ) 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/ ℎ ( )) 2 (Eq. (A.10)). Inverting the latter relationship, so that [(1/ ℎ ( )) 2 ] − 1 2 = ℎ ( ), In the present study (in Appendix A.3), we came to a novel (unprecedented) expression (Eq. (A.12)) for the ℎ ( ) function having the algebraic form, This form has been found to be extraordinarily well suited for the present purpose. The actual values of the respective Taylor series expansion coefficients, 2 ≤16 , are listed up to an order of 2 ≤ 16 in the right column of Table 1. to less than 1% throughout an unusually large high-temperature interval (of 0 ≤ < 8). The "snaky" solid and dotted curves show that the actual deviations of the approximate ̃( ) 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 ̃( ≥1 ; =4) ( ) (Eq. (3.7)) involving truncated versions of Debye's conventional lowcapacity representation, ( ) ( ) (Eq. (2.2)), of the two lowest orders, = 1 or 2, respectively.

Alternate Usage of High-and Low-temperature Debye Function Formulae
From Figure 1, we see (in analogy to Figure 8 of an earlier study [11]) that the deviations, ( ≥1) ( ) − ( ) > 0, of Debye's conventional low-temperature curves (Eq. (2.2)) are decreasing with the increasing number of summation steps, = 1, 2, …, and they are increasing with decreasing magnitude of the argument (= / ). In contrast to the latter, the deviations,  Table 2) that the individual crossing point positions, , and the corresponding deviations from the exact ( ) values, are rapidly decreasing with increasing (cf. the positions of the solid circles, for the cases of = 1 and 2, in the middle inset to Figure 2). Table 2 Positions of the crossings, , of the unprecedented high-temperature curve,  In this way, we come to rather good approximations, ̂( ; =7) ( ), for the Debye integral in terms of -specific combinations of complementary high-and low-temperature curve sections, Equivalently, we can also express the latter combinations in the more compact form, From Table 2, we see that by limiting the summation procedure for to only 4 steps ( = 4), one nevertheless comes to a rather gentle approximation for the Debye function, ( ) ≅ ( = 4 ; =7) ( ), 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.

Interpolation Formulae
It turns out that it is possible to find yet another way of constructing good analytical approximations for the Debye function integral, ( ) (Eq. (2.1)), which does not involve a sudden switch between the high-and low-temperature curve sections at their crossing points. This alternative way is based on the construction of suitable interpolation formulae that bridge the gap between a truncated version of Debye's low-temperature expression (Eq.  [1,11,12,20]. An appropriate interpolation formula for connecting this pair of non-crossing high-and lowtemperature curves can be chosen to be of the form where and are formula-specific interpolation parameters. An optimal approach of the interpolation curve, ̃( ) (Eq. (3.5)), to the exact ( ) dependence could be achieved by performing a corresponding least-mean-square fitting, which leads to the adjusted values of = 9.775 and = 0.84 for the respective interpolation coefficients. The remaining deviations (̃( ) − ( ))/ ( ), of Eq. (3.5) from the exact Debye function, ( ), 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 −3 /( − 1)) occurring in Debye's low-temperature expression, ( ) (Eq. (2.2)), into a conveniently constructed interpolation formula. Accordingly, we consider a combination of a truncated low-temperature expression of the type ( ) → 4 4 /5 3 − 3 /( − 1) with the same high-temperature expression as above, We have found an appropriate version of an interpolation formula for this combination to be given by the expression For the sake of completeness, we would like to mention that an even more pronounced (orderof-magnitude) reduction of deviations from the exact ( ) function, which is inherent to interpolation formulas of the type given by Eq. 2)), that is retaining at least one summation step, i.e., In this connection, we succeeded in realizing rather good simulations of the true ( ) function by an interpolation formula of the latter type (Eq. (3.7)) even for the two lowest-order versions, where Debye's low-temperature expansion, ( ) ( ) (Eq. (2.2)), has been truncated just after the first or second summation step (i.e., = 1 or 2). Appropriate interpolation coefficients for = 1 were found to be = 5.393 in combination with = 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 = 2, the adjusted interpolation coefficients turned out to be = 4.146 in combination with = 0.15, involving an order-of-magnitude reduction of deviations to less than ±0.0006% (cf. the corresponding solid ("snaky") curve in the middle inset in Figure 1). Thus, if necessary, one can also use an interpolation formula of the type given by Eq. (3.7) in order to come to a high-accuracy approximation for the Debye function integral provided that at least two summation steps are retained in the corresponding low-capacity expression, ( ≥2) ( ) (Eq. (2.2)). An advantageous feature of any one of the interpolation formulae displayed above (i.e., Eq. (3.5) to (3.7)) is due to the fact that the slopes (as well as the higher-order derivatives) of the approximate ̃( ) curves are continuous functions of their argument, = / . At variance to this, the slopes (including the higher-order derivatives) of approximate functions of the type ̂(

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, ( ) = ( / ) (Eq. (2.1)), are applicable only to those heat capacity datasets for, which the basic assumption of a constant Debye temperature, ( ) ≈ = constant, is realized to a good approximation (at least within a limited low-to-high temperature region). Such a special state of affairs can hardly be expected to be found for typical semiconductor materials ( < 3 eV), and it is only very rarely encountered even for wide-band-gap materials ( > 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, ( ), of diamond (visualized in Figure 2).

Figure 2
Magnitudes of the effective (caloric) Debye temperatures for diamond, ( ), which have been derived using two high-accuracy formulae (Eq. (4.1a) and (4.1b)) from the (isobaric) heat capacities, ( ), measured in the low-capacity region by Desnoyer and Morrison [17] () and in the high-temperature region by Victor [24] ( △) and Dinsdale [25] (). One can see, in particular, that the ( ) 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 < < 700 K.
Concerning the individual sets of ( ) values shown in Figure 2, we would like to point out that the ( ) data points represented by open circles () for the low-temperature region (18 K < < 278 K) correspond to their original ones () shown by Desnoyer and Morrison [17] (in their Figure 1). At variance to these well-known low-temperature ( ) data, the different sets of hightemperature ( ) data points shown here (in Figure 2) have been derived from the available sets of isobaric heat capacity values, ( ≥ 300 ), that were given by Victor [24] (△) and Dinsdale [25] (). The corresponding point-by-point transformations have been performed here for the individual ( ) vs. ℎ (∞) ratios, 0 < ( ) = ( )/ ℎ (∞) < 1, by means of two analytical transformation equations [11] (cf. Eq. (B.5) and (B.3) in [11]), i.e., (involving the parameters = 3.145972 and = 0.07037359) in combination with Note that the maximum inaccuracy resulting from these two of ( ) expressions amounts to only 0.13% (at the crossing point, = 0.5975, between Eqs. (4.1a) and (4.1b)).
Viewing the whole ( ) dependence shown in Figure 2, we see, on the one hand, that the sequence of low-temperature data points due to Desnoyers and Morrison's cryogenic data show a rapid fall from a → 0 limiting magnitude of about (0) ≈ 2223 K to magnitudes below 1900 K (in the vicinity of ≈ 150 K). In view of such a rapid change in the -dependent Debye temperature, it is obvious that it is impossible to describe the temperature dependence of the underlying cryogenic heat capacity of diamond, ( < 150 K) (cf. the ( ) list presented in the Appendix of [17]), in terms of a Debye function integral, ( / ) (Eq. (2.1)), with a fixed Debye temperature (whichever would be the magnitude chosen for ). At the same time, we see also from Figure 2 that the magnitudes of ( ) within a rather large interval (of about 140 K < < 700 K) are within a relative narrow interval of ( ) ≈ (1857 ±12) K (corresponding to the variations of ( ) values by less than 1.3%). From this observation one can expect that it should be possible to perform a reasonable fitting of the temperature dependence of the heat capacity, within this limited low-to-high temperature region, by a Debye function parameterized by a fixed Debye temperature, = const. ≈ 1857 K. In performing the calculations, for the isochoric (harmonic) part of the heat capacity, ℎ ( ) ≅ ℎ (∞) ⋅̃( ) , (4.2) and the associated isobaric heat capacities, ( ) (see below), we choose, within the present study, the primary interpolation formula proposed above, ̃( ) (Eq. 3.5), i.e., explicitly which appears to be directly applicable to a region of 170 K < T < 400 K (cf. the corresponding ℎ ( ) behavior shown by the dashed curve in Figure 2). At the same time, we see from Figure 2 that the actual (effective) ( ) values reach their local maximum in the region between approximately 400-500 K and subsequently decrease with increasing temperature. Such behavior is well known as a characteristic feature of non-negligible contributions of lattice expansion and anharmonicities to the observable (isobaric) heat capacities, ( ). 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 ( ) dependence by an analytical expression of the form where the material-specific expansion coefficients , = 1,2,…, are quantifying the combined effects of lattice expansion and anharmonicities. On the basis of Eq. (4.4) (in combination with Eq. (4.3)), we performed a least-mean-square fitting process of the combined set of isobaric heat capacity data points, ( ) [17,24,25], that are ranging within an unusually large temperature interval (from 180 K up to 5000 K; cf.  (for the region 230 K < < 5000 K, at least). An exemplary high-accuracy fit for the whole ( ) dataset (solid curve, due to Eq. (5.10)) had been obtained on the basis of the representative hybrid model displayed in Sec. 5. Also shown in the figure is the corresponding high-temperature behavior of the associated isochoric (harmonic) heat capacity, ℎ ( ) = 3 ⋅ ( ) (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, ℎ ( ), and the associated isobaric heat capacity, ( ). This is obtained by replacing the approximate (interpolated) Debye function, ̃( / ) (in Eq. (4.5) The corresponding theoretical ( ) → ℎ ( )/ 3 dependence is represented by the dashdouble-dot curve (---) in the inset of Figure 3.

Alternative High-accuracy Fit due to the Representative Hybrid Model
The actual variation of the Debye temperature, ( ), of the wide-band-gap as well as semiconductor materials is well-known to usually undergo rapid falls from their limiting ( → 0) magnitudes, (0), toward significantly lower minimum values, , which used to be located at temperatures of the order ≈ (0)/12, (in accordance with Figure 2 for diamond; see also the cryogenic ( ) curves shown for Si an Ge in [2,3], for III-V semiconductors in [4][5][6], and for II-VI semiconductors in [7][8][9][10][11]). Thus, it is impossible to simulate the temperature dependences of the measured / ℎ ( ) dependences in cryogenic (liquid-hydrogen) regions, 0 < < (0)/12, by means of a Debye function integral (2.1) with a constant Debye temperature (whichever magnitude would have been chosen for a fixed value). A conventional way of reasonable fittings for at least partial (strongly truncated) sets of cryogenic ℎ ( ) data is well known to be given in terms of an odd-order Taylor series expansion [6,12,14,17,[27][28][29][30][31][32], ℎ ( ) = 3 3 + 5 5 + 7 7 + ⋯ (5.1) (cf., in particular, the numerous applications of this simple analytical model to Group-IV, III-V, and II-VI materials in [12]). However, this truncated odd-order series expansion is continually found to be applicable only to narrow intervals of about 0 < < (0)/(25 ± 10) [12]. (Note that, in cases of non-negligible contributions of an electronic system to the heat capacity in the liquidhelium-region, an additional linear term [11,28,29], ( ) → 1 , still needs to be included in Eq. (5.1)). Suitable incorporation of this commonly observable cryogenic ℎ ( ) behavior (Eq. (5.1)) into theoretical heat capacity models thus represents the necessary condition for their eventual applicability to the whole -regions of practical interest, i.e., from the → 0 limit up to very high temperatures (up to melting points, if necessary). This condition is well-fulfilled by two qualitatively different types of duly general analytical models, namely, the representative hybrid model [3,10,11,26], on the one hand, and one of the two (essentially equivalent) versions of an analytical Non-Debye formula [6,32], on the other hand. By comparing the alternative fittings of various comprehensive ( ) 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 ( ) datasets.

(5.8)
Comparing the latter expression for 3 (in (5.8)) with the preceding one (in (5.7)), we see that the ratio between 1 and (0) is unambiguously determined by the weighting factor, 1 , for the cubic component [11], i.e., Finally, we consider the measured (isobaric) temperature dependences, ( ), for higher temperatures, which use to be co-determined by anharmonicity effects. The entire ( ) 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]  We have performed a complete high-accuracy fit of the combined ( ) datasets under study [17,24,25] for diamond (as shown by the solid curve in Figure 3). As a result of the multiparameter least-mean-square fitting process, we have obtained the following values for the effective Einstein temperatures of the three discrete phonon peaks under consideration: The corresponding theoretical ( ) curve (up to 5000 K) and the associated ( ) = ( )/ 3 dependence is represented by the solid curves (⎯⎯) in Figure 3.

Discussion
Above all, in this paper, we have developed an unprecedented algebraic formula, ℎ ( ) ( ) (Eq. (A.12)), for the high-temperature behavior of the Debye function. Using, in particular, a suitably truncated version of the latter, ℎ ( = 4) ( ) (in Eq. (4.5)), we were able to realize a good approach to the true ℎ ( > ) dependence for diamond (at least within an interval from about 230 K to about 600 K (cf. the dash-double-dot curve for the ( ) ≈ ℎ ( )/ 3 -dependence, in the inset of Figure 3). A somewhat closer approach to the true ( ) behavior within the region 180 K to 230 K could be achieved when we used, alternatively, the high-to-low-temperature interpolation formula, ̃( / ) (Eq. (4.3)), for the fitting. However, it is a matter of principle that the corresponding ( ) dependence (dash-dot curve) exhibits a plateau (approximate constancy) from 0 up to approximately 150 K. This typical Debye function feature is, of course, in striking contrast to the wellknown non-Debye behavior (see Eq. (5.1)) of the measured ( ) dependences in the corresponding cryogenic regions, 0 < < ≈ (0)/12.

Additional Quantitative Information Implied by the High-accuracy Fit via the Representative Hybrid Model
As described in Section 5, we have performed a very accurate fit, from approximately 20 K up to 5000 K (cf. the solid curves in Figure 3), for the whole combination of the partial ( ) datasets given for low and/or high temperatures [17,24,25]. This least-mean-square fit has been realized on the basis of the representative hybrid model [10,11,26], which is represented here by Eq. (5.10) in combination with Eq. (5.4). From the adjusted set of the model-specific parameter values (5.11), it follows that the average deviations of the theoretical ( ) 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 → 0 limiting behavior of the fitted ( ) and ( ) curves. Using the first relation in Eq. (5.7), a limiting ( → 0) value of 3 = (0) = 0.17694 µJK −4 mol −1 (as indicated by a horizontal ( → 0) line in the inset of Figure 3) is obtained. According to the first relation in Eq. (5.8), this 3 -value corresponds to the → 0 limiting value of (0) ≅ 2223 K (as indicated by a horizontal line in Figure 2). (Note that the same (0) 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 (0) value of about 2220 K [33][34][35] for diamond.
The maximum of the ( ) curve, = ( ) = 0.2995 µJK −4 mol −1 , has been found to be located at the characteristic point of = 174 K, in the vicinity of which the corresponding ( ) the curve shows a certain cubic dependence, ( ∼ ) → ( ) ⋅ 3 . However, the latter is to be clearly distinguished from the well-known → 0 (limiting) cubic dependence based on Debye's theory [1], ( → 0) → 3 3 (cf. Eq. (5.1)). The qualitative and quantitative differences between these two (local) cubic -dependences in the cryogenic region are due to the circumstance that the magnitudes of the respective proportionality factors, i.e., This value is in very good agreement with the earlier values of = 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 -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], ( ) = 3 2 / 3 (for the interval 0 < < ). For the first and second-order moments, this simple model invoked by Debye gives the values (1) = (3/4) ⋅ and (2) = (3/5) ⋅ 2 . Inserting these special (1) and (2) values into the general expression (Eq. (6.1)) for the dispersion coefficients gives for this conventional model, a magnitude of [26,37]  Thus, we see that the material-specific magnitude of the dispersion coefficient, (Eq. (6.1)), for diamond is only approximately 7% higher than its model-specific counterpart, (Eq. (6.2)), implied by Debye's original model [1]. The relatively similar values of to is the obvious reason of the relatively good functioning of Debye's very special heat capacity shape function, ( ) (Eq. (2.1)), for an approximate description of the ℎ ( > ) dependence (as shown in Figure 3).

Assessment of the Incomplete Fittings by Novel Formulae Involving a Fixed Debye Temperature
Taking the high-accuracy fit performed for a diamond based on the representative hybrid model (in Section 5) into consideration, let us now discuss in more detail the possible practical usability of two versions of the unprecedented Debye function formulas displayed in Section 3.
Firstly, concerning the limited fit of the experimental ( ) datasets, from 180 K to 5000 K, done using the interpolation formula, ̃( / ) (Eq. (4.3)), in combination with Eq. (4.4), we are concerned with average deviations of ±0.54% between the theoretical and experimental ( ) values (which are only slightly higher than those associated with the abovementioned fit by the hybrid model). Accordingly, the approximate theoretical ( ) and ( ) = ( )/ 3 dependences, which are represented by the dash-dot curves in Figure 3, are almost coinciding, for any > 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 -values, such as diamond (Eq. (6.1)). This interpolation model could, in particular, be of use within the framework of thermo-chemistry, where one is mainly interested in fine numerical descriptions of ( ) dependences within the high-region, T > 298.15 K. In addition, from the inset of Figure 3, it can be seen that this model, which exhibits a low-temperature plateau behavior of the respective, approximate ( )/ 3 function (dash-dotcurve), is absolutely inapplicable to < = 174 K.
With regards to the alternative fit of the experimental ( ) datasets, from 230 K to 5000 K, via the truncated high-temperature Debye function formula, ℎ ( = 4) ( ) (Eq. (3.1)), average deviations of ±0.74% between the theoretical and experimental ( ) 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 ( )/ 3 function (represented by a dash-double-dot curve, in the inset of Figure 3) that this simplified (algebraic) model nevertheless provides a relatively close fit (involving underestimation of the experimental values by less than 8%) in the region from approximately 80 K to 230 K.

Sub-quartic Behavior of the Low-temperature Heat Capacity Dependence
Such a relatively close fit to the actual (experimentally observed) -dependence, within a region of approximately /2 < < , is not surprising in view of the circumstance that thedependence of the truncated Debye function formula, ℎ ( = 4) ( / ) (Eq. (3.1)), is tending in the → 0 limit to a certain (fictive) quartic ℎ ( ) -dependence (cf. (Eq. (4.5)), ℎ ( ) → ( ℎ (∞)/√ 8 2 ) ⋅ ( / ) 4 . This peculiar feature is at least in qualitative accordance with the global expectation [38], according to which, in the cryogenic region below , the effective exponents, , of the power function tangents, ( ) ∝ , should range somewhere between 3 and 4 (for a wide-band-gap material such as diamond; cf. Figure 7 in [38]). In order to estimate the maximum value of this exponent, , let us consider for the vicinity of the point of inflexion, = 99.3 K, of the fitted ( ) curve, an approximation by a linear tangent [38],  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], Using the general expression [38] for the -dependence of the effective power function exponents, for the effective (maximum) power function exponent in the vicinity of the point of inflection (as indicated by a solid triangle in the inset of Figure 3). Thus, in the case of diamond, we are concerned with a typical sub-quartic behavior of the cryogenic ( ) dependence. This is in accordance with our preliminary assessment shown for diamond in Figure 7 of an earlier study [38].

Comparison with a Recently Proposed Multi-Debye-function Model
Of interest might be a comparison with a recently published model [39] for heat capacity fitting, from cryogenic to high temperatures, on the basis of a certain interpolation model consisting of a combination of three Debye function components, for low-temperature regions, with the popular (thermo-chemical) Maier-Kelley approximation [40], for high temperatures. Concerning the case of diamond, one can see from a comparison, e.g., of the ℎ ( ) 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 values ( 1 , 2 , and 3 ) estimated by Vassiliev and Taldrik [39] (which are listed, together with the respective weighting factors, , e.g., under number 1b in Table 6 of their article). The average magnitude, ̄= 1 1 + 2 2 + 3 3 , of this set of values amounts to ̄= 1896 K. The latter is by only 2.2% higher than our single Debye temperature value of = 1855 K, which resulted from the present fitting using Eq. (4.2) and/or Eq. (4.5). Furthermore, we see that the individual values estimated by Vassiliev and Taldrik [39] differ from their average value, ̄, by less than 4%. The relatively small differences between the individual 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 ( > ) data presented in Section 4.

Possible Usefulness of a Properly Adopted Non-Debye Interpolation Formula
Satisfactory numerical fitting to the comprehensive ( ) datasets for the high-dispersion wide-band-gap materials GaN and ZnO [36,37], from very low up to high temperatures, had been performed for the first time in a previous study [32] on the basis of a novel algebraic non-Debye formula. A suitable alternative version of the latter had been displayed and used in another work [6] for comprehensive fitting (from liquid-helium temperatures up to the melting points of the materials) of the combined ( ) datasets, that are available for a series of Group III-V (highdispersion [37]) semiconductor materials.
We have found that an appropriate specialization of the preceding analytical non-Debye model [6] with respect to its possible applications to low-dispersion materials, such as diamond, can readily be established by adopting an algebraic expression of the interpolation type for the isochoric heat capacity shape function (the structure of which is analogous to the structure of the original non-Debye heat capacity shape function (A3) considered in [6]). Envisaging Eq. (6.8) for calculations of the respective isochoric heat capacities, ℎ ( ) = ℎ (∞) ⋅ ; ( ), (6.9) we see, on the one hand, that at sufficiently low temperatures, where the contributions of the first four ∝ −2 terms ( = 1 to 4, in the denominator of Eq. (6.8)) are small compared to the contribution of the ∝ −14 term, the ℎ ( ) expression in Eq. (6.9) approaches asymptotically to the familiar odd-order Taylor series expansion in Eq. (5.1). On the other hand, at sufficiently high temperatures, the numerator tends to unity and the denominator approaches asymptotically to the truncated high-temperature expression given by Eq. (4.5). However, a characteristic deviation in the latter case is due to the circumstance that the fourth expansion coefficient, 8 , for the ∝ −8 term, must be admitted to deviate more or less significantly from its (originally fixed) Debyefunction-related predecessor, 8 (Table 1). This generalization is inevitable for assuring a smooth change between the respective asymptotic low-and high-temperature behavior.
Within the framework of this unprecedented interpolation model, the isobaric heat capacities, ( ), are represented by an analytical function of the usual type: values. Furthermore, we see that the fitted magnitude of the Non-Debye expansion coefficient, 8 , is a factor of ~3 higher than its preceding counterpart, 8 , due to the original Debye function expansion (cf. Table 1). This considerable difference is obviously due to strongly increasing non-Debye effects in the vicinity of ≈ (0)/12.
Concerning the associated estimation of the → 0 limiting value of the Debye temperature, we see that the alternatively fitted 3 -value of 3 = 0.17907 µJK −4 mol −1 corresponds (according to Eq. (5.8)) to a (0) value of about 2214 K. The latter is approximately 0.3% lower than the ap-parently 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 ( ) curve to be = ( ) = 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 -position of the maximum, = ( ), resulting from Eq. (6.10) is located at ≈ 164 K (which is by ~6% lower than the -position estimated via the hybrid model used in Section 5). Finally, we would like to mention that the average deviations of the theoretical ( ) values resulting from the fit of the experimental points under consideration using Eq. (6.10) (in combination with Eq. (6.8)), amount to approximately ±1%, which corresponds to the same order of magnitude as the typical experimental uncertainties in the low-temperature region. In view of these still moderate deviations between the theoretically fitted and experimental ( ) values, one can conclude that the latter version (Eq. (6.10), in combination with (6.8)) of an algebraic non-Debye interpolation formula can be successfully applied to numerical fittings of comprehensive ( ) datasets, from the liquid helium region ( → 0) up to very high temperatures, provided that the dispersion coefficient, (Eq. (6.1)), of a wide-band-gap material under consideration is not significantly higher than those estimated for diamond and c-BN [36,37].

Concluding Remarks
We have displayed within this study a variety of structurally relatively simple interpolation formulae that are convenient for fitting the measured low-to-high temperature heat capacity datasets of certain wide-band-gap materials that are characterized by a sufficiently low degree of phonon dispersion (0.258 < < 1/3; cf. Eq. (6.1) and (6.2)). A characteristic empirical feature of such low-dispersion materials is the quasi-plateaubehavior of the effective (caloric) Debye temperature, ( ), which is observed in T-regions from (0)/12 up to an order of (0)/3 (cf. Figure 2 for diamond; see also Figure 2 in [19] for c-BN).
Looking among a larger variety of papers on thermal properties of wide-band gap-materials for further potential targets for the applications of the presently developed analytical apparatus, the case of MgO [41][42][43][44] is the most prominent. Indeed, one can see from the corresponding ( ) and/or ℎ ( ) 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 → 0 limiting level of (0) ≅ 946 K to a minimum of ≈ 750 K (±10 K), the further run from 100 K up to ~700 K takes place at a nearly constant Debye temperature of ≈ 760 K. Thus, it should be possible to describe the corresponding ( ) 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 ≈ 760 K. In contrast, the obviously unconsidered choice of the limiting ( → 0) Debye temperature of (0) = 946 K reported in previous studies [45,46] can hardly be considered as an adequate -value for a proper description of the actual ( ) dependence of MgO (from the cryogenic region up to 2100 K).
Clearly, the arbitrary choice of → 920 K for ZnO [45,46] is also inadequate. This is due to the circumstance that the ( ) dependence of ZnO shows a very pronounced "snaky" shape (in analogy to the typical ( ) 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 ℎ ( ) values for ZnO extends from a lower boundary of ≈ 350 K up to an upper boundary of ℎ (∞) ≈ 690 K. A considerably higher -value, such as → 920 K, suggested in previous studies [45,46], for ZnO, ranges clearly outside the realm of reason.
One can make a pre-selection of groups of materials, the heat capacities of which should actually have (or not have) a chance to be reasonably described by the presently developed Debye function interpolation formulae, on the basis of their material-specific dispersion coefficients, (Eq. (6.2)). On viewing, e.g., a list of the corresponding values given for a large variety of wideband-gap materials (in Table 1 of an earlier study [36]), which had been estimated on the basis of various theoretically calculated PDOS spectra, we see that diamond and c-BN are the only two materials for which their dispersion coefficients, ≅ 0.275 ±0.005, are closely approaching the Debye's limiting value of ≅ 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 ( ) datasets for these two materials. In contrast, the values (> 0.35) listed in Table 1 (of Ref. [36]) for all the other wide-band-gap-materials under consideration, are throughout > 40% higher than . In particular, it should be noted that the values estimated for GaN and ZnO are even a factor of ~2 higher than . 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, = const.). Analogously, we find for a larger variety of typical semiconductor (Group-IV, Group-III-V, and Group-II-VI) materials considered [37] that their material-specific values are throughout > 60% higher than (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 ) 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 values very similar to Debye's value, ≅ 0.258. This concerns the two alkaline earth oxides MgO ( ≈ 0.29) and CaO ( ≈ 0.30) as well as the two alkali halides NaF ( ≈ 0.28) and NaCl ≈ 0.30). Thus, it should be possible to directly apply the presently developed analytical apparatus, if necessary, to the latter wide-band-gap materials.
Finally, we would like to point out that of primary importance for numerical simulations of the ( ) datasets available for such low-dispersion materials is the unprecedented interpolation formula of the mixed type, ; ( ) (Eq. (6.8)), which smoothly combines the characteristic hightemperature Debye model behavior (Eq. (3.1)) with the qualitatively different limiting lowtemperature non-Debye behavior (Eq. (5.1)). The actual use of this ambiguous heat capacity shape function, ; ( ), 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 ( ) datasets for any temperature of practical interest, i.e., from the liquid-helium region (including the → 0 limit) up to the melting points of the low-dispersion materials.
The crucial problem of reasonable theoretical interpretation of the measured ( ) datasets for a large variety of wide-band-gap materials, in combination with an adequate analytical description of the qualitatively different Debye temperature behavior, is intended to be treated in due detail in a forthcoming study.
We have shown in Figure 6 of the appendix of an earlier work [22] that, e.g., the deviations of a truncated ℎ ( ) expression of this exponential type (Eq. (A.5)),

A.2 High-Temperature Representation for Reciprocal Debye Function Values
In the course of the forthcoming analytical studies, we succeeded in establishing a few more (structurally different) versions of the high-T representations for the Debye function (Eq. (2.1)), which turned out to show even markedly better convergence properties than the exponential one sketched above. An interesting example had been presented, firstly, in [11] with respect to the reciprocal values, [ ℎ ( )] −1 , of the ℎ ( ) function (Eq. (2.5)). Accordingly we have made, in a previous study [11], the alternative ansatz  A.3)). Evaluating the respective expansion coefficients associated with the individual 2 power terms (occurring in the argument of the (1 + ) −1 function, Eq. (A.7)), and identifying them with the novel (transformed) expansion coefficients, 2 , we have obtained, in [11], the following values for the latter:  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, ℎ ( ), of the Debye function (2.1)), in terms of the reciprocal values of a corresponding Taylor series expansion, ( . 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 ℎ  Table  7 of Ref. [11]).

A.3 Unprecedented Series Representation for Reciprocal Square-root Debye Function Values
Another crucial example of a rapidly converging Taylor series expansion which turned out to play an especially useful role within the present context, concerns the square of reciprocal values, (1/ ℎ ( )) 2 , of the ℎ ( ) function (A.1).
Accordingly, we make the alternative ansatz (The values of the subsequent four expansion coefficients, up to an order of 2 = 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, ℎ ( ), of the Debye function (2.1)), in the form of the following algebraic expression: The advantageous properties of the latter representation have been visualized and discussed in detail in Section 3.
Here, we have used the ℎ 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.