Detailed analytical and numerical analyses are performed for combinations of several complementary sets of measured heat capacities, for ZnSe and ZnTe, from the liquid-helium region up to 600 K. The isochoric (harmonic) parts of heat capacities, , are described within the frame of a properly devised four-oscillator hybrid model. Additional anharmonicity-related terms are included for comprehensive numerical fittings of the isobaric heat capacities, . The contributions of Debye and non-Debye type due to the low-energy acoustical phonon sections are represented here for the first time by unprecedented, integral-free formulas. Indications for weak electronic contributions to the cryogenic heat capacities are found for both materials. A novel analytical framework has been constructed for high-accuracy evaluations of Debye function integrals via a couple of integral-free formulas, consisting of Debye’s conventional low-temperature series expansion in combination with an unprecedented high-temperature series representation for reciprocal values of the Debye function. The zero-temperature limits of Debye temperatures have been detected from published low-temperature data sets to be significantly lower than previously estimated, namely, 270 (±3) K for ZnSe and 220 (±2) K for ZnTe. The high-temperature limits of the “true” (harmonic lattice) Debye temperatures are found to be 317 K for ZnSe and 262 K for ZnTe.

1. Introduction

Since the invocation of the concept of apparently characteristic, material-specific temperature parameters, , within Debye’s classical paper [1] on specific heats of solids, one was concerned with a large variety of quotations of corresponding values (so-called “Debye temperatures”) within numerous thermophysical research papers, including various representative review articles [25] and books [610]. A brief inspection of the enormous data material published hitherto, however, showed readily that it is obviously not possible to find unique values for the individual materials in question, which should apply to duly broad temperature regions. In order to save, nevertheless, the apparently rather popular concept of Debye temperatures, it became thus already long ago a widely accepted custom to adopt the physically somewhat problematic notion of variable (-dependent) Debye temperatures, . Accordingly, Debye’s heat capacity model function [1] used to be employed in practice in the general form [3, 4, 6, 1015]where the upper limits of integration are given by the ratios, , of adjustable Debye temperatures versus the respective lattice temperatures and represents the familiar Dulong-Petit limiting value for the isochoric lattice heat capacity in harmonic approximation (i.e., in particular , for binary materials).

In the course of numerical fittings of measured (isobaric) heat capacities, , on the basis of (1), it has continually been found that, in contrast with Debye’s original suggestion [1], proper simulations of such heat capacity curves can in fact be realized only by admitting rather strong -dependencies of the material-specific Debye temperatures, . This statement applies above all to the liquid-helium-hydrogen region, where the adjusted values are as a rule rapidly falling [35, 1217] from their limiting levels, , to certain minimum values, , which used to be located in the vicinities of . Furthermore, towards higher temperatures, many curves show a more or less pronounced increase up to certain material-specific maxima, , the actual magnitudes of which are in many cases significantly higher than the respective levels. Such typical nonmonotonic behaviours of curves in the cryogenic region have been found in particular for Si and Ge [1821] as well as for a large variety of III–V [11, 12, 15, 2230], II–VI [14, 3138], and I–VII [16, 39, 40] materials.

Globally one can thus assess that Debye’s original idea, according to which heat capacities of solids might be represented by Debye function integrals (of type (1)), whose upper limits of integration should be based on fixed, material-specific values, is largely illusionary. Actually, not a single binary or ternary material could be found hitherto, within a wealth of thermophysical researches, for which it would have been possible to simulate in adequate way the temperature dependence of measured heat capacities, from the liquid-helium region up to the vicinity of room temperature, in terms of a unique value. This means, among other things, that there is hardly a chance to come to appropriate quantitative descriptions of the temperature dependencies of heat capacities mainly via constructions of various, highly elaborate analytical models [4144] for integrals of Debye’s type and imputing, at the same time, only some fixed (-independent) values for Debye temperatures (like = 632 K for GaN [43] or = 920 K for ZnO [44] and = 946 K for MgO [44]).

The basic cause of the occurrence of relatively deep minima, , in the cryogenic region had been indicated already many years ago by Schrödinger [2] to be obviously due to the onset of thermal activation of short-wavelength transverse acoustical (TA) phonons, which are quantum-theoretically manifested by pronounced peaks (singularities) in the respective material-specific phonon density of states (PDOS) spectra (see, e.g., the calculated spectra shown in [45, 46] for ZnSe and ZnTe). Such drastic deviations of physically realistic PDOS spectra from Debye’s naïve (quadratic) PDOS model function [1], , are continually confirmed to be indeed the main cause of the obviously typical, nonmonotonic (non-Debye) behaviours of curves.

Of course, when a representative set of properly determined values is actually available for one or the other material in question, one can readily calculate the respective heat capacities, either via direct numerical evaluations of the integral occurring in (1) or by interpolations between discrete values quoted in some of the more or less comprehensive or tables presented by various authors [6, 8, 4749] (see also the list of presently obtained high-accuracy sample values given in Table 7).

Considerably more problematic and difficult, however, is a proper computational solution of the hitherto largely ignored inverse problem, namely, the reliable determination of effective Debye temperatures, , on the basis of measured heat capacities, . In view of the longstanding absence of corresponding integral-free formulas, it was the custom in thermophysical investigations over a long time to resort to some of the tables quoted above and to perform thus approximate estimations of the respective -ratios via interpolations between discrete -ratios quoted for the couples of the two neighbouring (or ) values. Such conventional interpolation procedures were still feasible as long as one used to be concerned only with relative small sets (limited to few dozens) of measured data points. However, in contemporary literature, one is repeatedly concerned with rather large sets (of hundreds) of data points (like the presently analyzed data sets given in [46] for ZnSe and ZnTe). With respect to the latter such clumsy point-by-point interpolation procedures appear to be largely obsolete. A practicable computational alternative, however, requires the availability of a duly elaborate analytical framework for the inverse relationship, which should hence enable performing straightforward evaluations of the actual -ratios directly on the basis of given heat capacity values, . It is the main purpose of the present study to develop here (within Appendix B) such a representative analytical framework of integral-free formulas and to present corresponding sample results for Debye temperatures of ZnSe and ZnTe.

Concerning the successive development of this novel analytical apparatus for Debye temperatures, we would still like to note that some partial results in form of asymptotic (approximate) Debye temperature expressions had already been published in two preceding papers, namely, in [13] for the low-temperature region (),and in [15] for the complementary range of higher temperatures (),

However, due to a lack of space, the somewhat unusual derivation procedure for this highly useful high-temperature formula (i.e., Eq. in [15]) could not be included into the former paper. A duly detailed display of the corresponding analytical derivation procedure is now made up, among other things, in Appendix B.

In Section 2 we briefly sketch the basic arguments and equations of the analytical framework for theoretical calculations of isobaric heat capacities. In Section 3 we display a reformulated, integral-free version of the multioscillator hybrid model [14, 21, 50]. This is based mainly on the usage of novel, unprecedented analytical formulas (derived in Appendix A) for the contributions of the continuous, long-wavelength sections of transverse acoustical (TA) phonons to the resulting lattice heat capacities. In Section 4 we adopt (in analogy to [14]) a special 4-peak Einstein oscillator constellation of the hybrid model, and we use the respective analytical heat capacity expression for careful fittings of compatible sets of low- and high-temperature data available for ZnSe and ZnTe (up to 600 K). In Section 5 we perform, on the basis of a couple of precision formulas derived in Appendix B, transformations of the experimental heat capacity data in question (including the corresponding theoretical and curves), into the respective Debye temperature representations. The results are discussed in Section 6.

2. Basic Equations for Temperature Dependencies of Heat Capacities in Solids

Within the frame of the harmonic lattice regime, the temperature dependencies of the isochoric heat capacities, , are well known to be given by corresponding heat capacity shape functions, , of the general form [4, 6, 45, 5055]where we have denoted by the energies of the individual phonons and represents the material-specific PDOS spectral function [3, 4, 19, 56]. The latter is supposed within the present context to have been normalized to unity [6, 45, 52, 53, 55]; that is, .

Within the frame of numerical analyses of experimental heat capacity data, however, one is continually concerned with a basic theoretical complication due to the inherent differences between the temperature dependencies of theoretical functions, on the one hand, and those of measured (isobaric) heat capacities, , on the other hand. The latter are well known to be throughout somewhat higher, , than their theoretical (isochoric) counterparts pertaining to the harmonic lattice regime [4, 50, 54]. The respective differences, , are usually found to be very small from 0 up to temperatures of order , where the heat capacity amounts to about 50% of the Delong-Petit limiting value, . On the other hand, one is usually concerned with a relatively strong monotonic increase of these differences, , at higher temperatures (). These differences are generally ascribed to cumulative effects of lattice expansion and lattice anharmonicities [3, 11, 18, 21, 23, 39, 50, 54, 56, 57]. We have repeatedly found within a larger series of heat capacity studies [14, 15, 21, 50, 54] that, within regions from up to temperatures of order (at least), the rapid increase of these differences can be simulated in good approximation by a proportionality of type (note that the structure of the latter is analogous to the known Nernst-Lindemann formula [2, 3, 58], , for the differences between isobaric and isochoric heat capacities). Admitting, furthermore, that at sufficiently high temperatures (up to melting points, if necessary) some additional higher-order power terms of may come into play, it appeared reasonable to simulate the -dependencies of isobaric heat capacities by a duly general algebraic expression of the form [14, 15, 21, 50, 54, 59]where the empirical Taylor series expansion coefficients, , are quantifying the material-specific weights of contributions of the individual, anharmonicity-related power terms [15, 21, 54, 57]. Furthermore we have still included into (5) an additional linear term, , which has to represent possible contributions of a degenerate electronic system [4, 6, 13, 23, 6063] to the measured heat capacities. Such electronic contributions (“Sommerfeld terms” [4, 6, 60]) are known to be rather strong (tending to become even the dominant ones, at sufficiently low temperatures) in metals, due to large concentrations of conduction electrons. In contrast to this, such electronic contributions are as a rule rather weak in semiconductors and insulators due to the relatively low concentrations of free electrons associated with lattice imperfections (impurities). Accordingly, the actual strengths of such electronic contributions can vary considerably even between different samples of one and the same semiconductor material [4, 62]. We will find below (in Section 4) indications for some weak electronic contributions to the isobaric heat capacities measured at liquid-helium temperatures [46], for both ZnSe and ZnTe.

3. Multioscillator Hybrid Model and Its Integral-Free Representation

Physically realistic, material-specific phonon density of states (PDOS) functions, (like those shown, e.g., in Figure  3 of [45] or in Figure  4 of [46], for ZnSe and ZnTe), which are determining the -dependencies of the respective harmonic heat capacity shape functions, (see (4)), are known to consist in general of numerous peaks (critical points) in combination with various continuous curve sections and/or gaps between them. Many fine details of such material-specific PDOS spectra, however, are practically not detectible via analyses of experimental heat capacity data [3, 64]. Their possible influence on thermal properties used to be strongly reduced by the thermal averaging process (in (4)), so that they are not discernable from the background of accidental experimental uncertainties (being usually of the order 1%). Consequently, it is as a rule possible to realize rather good fittings of available data on the basis of relatively simple theoretical model functions for PDOS spectra [64], provided that the latter are sufficiently general and flexible for reasonable reproductions of various prominent features, at least. This concerns above all the liquid-helium-hydrogen region [13, 64], where the lattice heat capacities are known to be represented in good approximation by Taylor series of odd-order terms [13, 16, 18, 23, 35, 39, 65]

This limiting behaviour is the automatic consequence of the well-known fact that the low-energy (TA phonon) tails of PDOS spectra for three-dimensional crystals are generally given by Taylor series [13, 16, 18, 23, 35, 39, 64, 65] of exclusively even-order power terms, , of the phonon energy. An incorporation of a corresponding continuous low-energy section into theoretical model functions is therefore indispensable for their actual applicability to the cryogenic region. At the same time it is also clear that such continuous low-energy tail sections are naturally ending (i.e., ) at the positions of the first critical points [14, 21, 50, 64], , which are usually corresponding to the first peaks of short-wavelength TA phonons (note that, for materials with zinc blende structure like ZnSe and ZnTe, the energy positions of the first critical points are corresponding to those of the L-points on the first Brillouin zone boundaries, i.e., [7, 32, 36, 45, 46, 66, 67]).

3.1. General Analytical Framework of Multioscillator Hybrid Models

Less critical than the general ansatz for the region is the modelling of spectra for energies higher than the first critical points, . Numerous analyses have shown that excellent simulations of dependencies (4), from cryogenic to room temperatures, can be realized on the basis of combinations of several (macroscopic) Einstein oscillators [68], ( = 1 to ), the ordered energy positions of which, , including their actual weights, , are to be properly adjusted in the course of least-mean-square fitting processes for heat capacities. Thus one can take the normalized PDOS function in (4) to be generally representable by a multioscillator hybrid model function of the sufficiently general form [14, 21, 50]where we have denoted by ( for and 0 for ) the step function, which realizes the cut of the continuous part of the function just at the first singular point. By integrating (7) over the whole spectrum (from 0 beyond the highest Einstein oscillator), we see that the supposed normalization of the function to unity is equivalent to a normalization condition of the form [14, 21, 50]for the whole set of weighting factors. This condition is to be strictly fulfilled within any fitting process of experimental data on the basis of the model function (7).

Inserting, finally, the function (see (7)) into the integral (4) and denoting by ( = 1 to ) the characteristic Einstein temperatures corresponding to the energy positions of the individual macroscopic oscillators we obtain for the normalized (harmonic) heat capacity shape function the general expression [14, 21, 50]where the functions (being due to the two continuous components, , in (7)) are given in form of integrals of Debye and non-Debye type [14, 21]

Here we have denoted bythe ratio of the characteristic phonon temperature associated with the first (lowest) special point, , versus lattice temperature, (in analogy with the notation used within Debye’s theory, ; cf. (1) and Appendix A).

Informative byproducts of numerical fittings of measured data sets via (5), in combinations with (9) for the harmonic heat capacity shape function, (including (10), for the components), are the moments, , of material-specific PDOS spectral functions [3, 4, 16, 18, 19, 21, 50, 52, 56, 6972], which are generally defined by

Inserting the hybrid model function (see (7)) into the integral for (see (12)), we obtain for these moments the model-specific expression [14, 21]

3.2. Integral-Free Formulas for the Contributions of Long-Wave Acoustical (TA) Phonons

On the basis of the hitherto used original (i.e., integral) expressions (see (10)) for the two low-phonon-energy contributions, (in (9)), the corresponding least-mean-square analysis procedures performed in several preceding studies [14, 21, 50] turned out to be somewhat heavy and time-consuming. This disadvantage of the former hybrid model versions [14, 21, 50] was one of the main motivations for performing here (in Appendix A) a comprehensive analytical study aiming at the construction of more practicable, integral-free formulas for the integrals (10) of Debye and non-Debye type.

Consider first the integral (10) for the low-energy component of Debye type, (see Figure 1). A detailed analyses (in Section A.2 of Appendix A) has shown that the corresponding dependence (see (A.13)) can be represented in good approximation in terms of a proper combination (A.16) of truncated low- and high-temperature expansions, (see (A.14)) and (see (A.15)), assuming a switch between the two complementary branches just at their crossing point, = 6.5135 (cf. the upper inset of Figure 1). This means for the low-temperature branch, ,and for the high-temperature branch, , (the respective expansion coefficients of which, , are listed in Table 6).

Concerning the integral (10) for the second low-energy component, (see Figure 1), we have found within the frame of a detailed study (in Section A.2, of Appendix A) that the corresponding dependence of non-Debye type (see (A.17)) can be represented in good approximation in terms of the combination (A.24) of the truncated low- and high-temperature expansions, (see (A.18)) and (see (A.23)), assuming again a switch between the two complementary branches at the respective crossing point, (cf. the lower inset of Figure 1). This means for the low-temperature branch, , explicitlyand for the high-temperature branch, , (the respective expansion coefficients of which, , are listed in Table 6).

The large scale representations of the and curves in Figure 1 are suggesting apparently perfect coincidences between the exact (solid) curves and the respective approximate curve sections. At the same time we see, from the insets of Figure 1, that the approaches to the exact functions are really excellent only in the regions < 5 and > 9, whereas one is concerned in the vicinities of the crossing points, that is, at = 6.5135, for the Debye-like component, and at = 7.0962, for the non-Debye component, with certain overestimations of the exact values by the approximate formulas (up to a maximum order of 0.008%, for = 1, or of 0.015%, for = 2). Yet, such computational inaccuracies are by 1 to 2 orders of magnitude smaller than the typical uncertainties of measured heat capacities, . Thus, these computational inaccuracies are usually negligible within the frame of practical applications.

3.3. Connection with Conventional Cryogenic Heat Capacity Parameters

At sufficiently low (liquid-helium) temperatures, where , the first terms ( or ) in the respective low-temperature expressions (see (14) and (16)) are strongly dominating. Thus the Debye-like and non-Debye components, ((see (10) for = 1 and 2, resp.), tend here to the known and asymptotes [14, 21]

The asymptote of the component is of Debye type [1], whereas the asymptote of the component is responsible for the onset of the typical non-Debye behaviour. The latter is frequently emerging already at very low temperatures (of order ), where the Einstein oscillator terms in (9) for the dependence are practically negligible (owing to their plateau behaviour, ; cf. also the curve in Figure 1). In this way we satisfy ourselves that the low-temperature dependence of the harmonic lattice heat capacity, (see (4)), reduces, as expected, to a couple of two odd-order terms [13, 16, 18, 23, 35, 39, 65], (see (6)). The respective expansion coefficients, and , are readily seen to be given in terms of the weighting factors and , the characteristic phonon temperature associated with the first Einstein peak, , and the classical Dulong-Petit heat capacity value, , by the relations

Comparing, finally, the analytical expression for the limiting (cubic) asymptote implied by the present formalism, (due to (19)), with its conventional predecessor, (due to Debye’s theory [13, 8, 13]), we see that these two expressions are coincident when the weighting factor for the Debye-like component, , is just given by the cube [14, 21] of the ratio between the characteristic phonon temperature associated with the first Einstein peak, , and the limit value, of the conventional Debye temperature; that is,

4. Sample Analyses of Heat Capacity Data Sets for ZnS and ZnTe

For the sake of numerical analyses of given experimental heat capacity data within the frame of a multioscillator hybrid model like that displayed above, it is still necessary to specify the material-specific constellation (in particular the total number, ) of Einstein oscillators to be taken into consideration. A relatively simple version in form of a three-oscillator-based hybrid model [21] had been used in [50] for a large variety of III–V and II–VI materials with partly different crystal structures. Within the present (more detailed) study, for ZnSe and ZnTe, we shall consider a structure-specific constellation of four discrete oscillators (in analogy to [14]), which offers a physically adequate association with the contributions of the 6 phonon branches in zinc blende structure materials (consisting of 2 TA branches and 1 LA branch, for the lower parts, and 2 TO branches and 1 LO branch, for the upper parts, of the respective PDOS spectra [7, 32, 36, 45, 46, 66, 67]). Let us briefly sketch the basic arguments in favour of the corresponding four-oscillator model [14].

4.1. Four-Oscillator Constellation Chosen for ZnSe and ZnTe

In the course of numerous analyses of experimental data sets by different versions of the multioscillator hybrid model [14, 21, 50] we have repeatedly made the experience that a necessary condition for the contributions of two oscillator peaks to be actually discernible, in the course of least-mean-square fitting processes, consists in a sufficiently large distance (of more than 10%, at least) between their respective ( versus ) energy positions. Considering under this aspect the known positions of the dominating optical ( versus ) phonon peaks [7, 32, 36, 45, 46, 66, 67] in ZnSe and ZnTe we see that their relative distances are markedly smaller than 10%, so that their individual contributions to the measured dependencies are practically indiscernible. Consequently, the contributions of the three optical phonon branches can be lumped together, being thus representable within the four-oscillator hybrid model [14] by a single oscillator with the highest energy, , and a fixed weight of

Qualitatively quite different is the state of affairs for the three acoustical phonon branches. From the theoretically calculated PDOS spectra [7, 32, 36, 45, 46, 66, 67] one can see that the distances of the dominating short-wave LA peaks from the prominent TA peaks (as well as from LO/TO peaks) are larger than 30%. The contributions of the dominating LA peaks, (= centres of gravity of the LA branches), are thus well discernible from the contributions of TA phonons as well as from those of the LO/TO phonons. This corresponds in reasonable approximation to a fixed weight offor the -oscillator. Finally, with respect to the 2 TA phonon branches, we observe that the energy differences, for example, between the special points and , are larger than 20% [7, 32, 36, 45, 46, 66, 67], so that their individual contributions have good chances to be still discernible within the frame of duly careful numerical analysis procedures. This observation suggests representing the whole variety of short-wave TA phonons by two discrete model oscillators, ( = 1 and 2). Concerning their actual weights, and , we observe that, owing to the preceding fixations of the weights for the LA and LO/TO oscillators (i.e., (21) and (22)), the global normalization condition (8) for the weights of all contributions associated with the two TA phonon branches (i.e., and , in combination with and ) reduces, consequently, to a normalization condition of the TA phonon-related weighting factors of the special form [14]

Accordingly, only 3 of these 4 weighting factors are actually playing the role of independent fitting parameters; that is, the magnitude of one of them is automatically fixed by the set of the three other ones (e.g., ). In summary, within this special 4-oscillator hybrid model [14], which turns out to be obviously well suited for binary materials with zinc blende structure, the normalized heat capacity shape function, (9), reduces to a parameterized function of just 7 independent parameters (consisting of the 4 discrete phonon temperatures ,   = 1 to 4, in combination with the three weighting factors , , and ).

Furthermore we would still like to observe that, in many cases, one knows at least approximate values for the 0 limiting values of Debye temperatures, (either as results of former analyses of cryogenic heat capacity data [13, 34] or from independent estimations on the basis of elastic constants [35, 7377]). This preliminary knowledge allows us, if necessary, to substitute the weight, , of the Debye-like component, , by (according to (20)).

4.2. Additional Empirical Parameters Involved by the Isobaric Heat Capacities

It has already been pointed out above (in Section 2) that, for numerical simulations of measured isobaric heat capacity data, , by means of (5), it is generally necessary to envisage the inclusion of further empirical parameters (like , , and/or ) into the fitting process. Such additional parameters have been introduced above for a quantification of the deviations, , of measured values from the theoretical (harmonic) lattice heat capacities, (see (4)), due to lattice anharmonicities and possible electronic contributions. The actual inclusion (or exclusion) of one or the other additional parameter depends, of course, on the concrete constellation of the whole set of experimental data in consideration.

The crucial experimental basis of the present heat capacity analyses, for ZnSe and ZnTe, was provided by Kremer et al. [46] in form of unusually comprehensive and informative sets of data points comprising the regions from 2 K up to 277 K, for natural ZnSe, and up to about 150 K, for natural ZnTe (see the sets of and data points [46] represented by empty circles in Figures 2 and 3, resp., and cf. Figure in [46]). Yet, particularly from the inset of Figure 2 we see that the sequence of the respective data points () for ZnSe are indicative of some nonmonotonic (local minimum) behaviour in the liquid-helium region. This is at variance in particular to earlier data (+) given for the interval 1.7 K to 24.8 K by Birch (see Table of [34]), which showed an exclusively monotonic behaviour of the sequence of the respective data points (+) in the same region.

The more or less pronounced nonmonotonic behaviours of Kremer’s data points for natural ZnSe (Figure 2) and ZnTe (Figure 3), in the vicinities of 3.7 K or 2.6 K, respectively, can be ascribed to relatively weak (nevertheless not completely negligible) contributions of degenerate electronic systems to the measured dependencies (in analogy to similar observations reported for semiconductor materials like Si [61], Ge [62], several III–V materials [23], and FeGa3 [63]). Such electronic contributions are analytically represented within the present context by the respective linear (electronic) term, , in (5).

The subsequent least-mean-square fitting processes are thus to be performed on the basis of (5) in combination with (9), (14) to (17), and (20). For present purposes we shall limit the sets of high-temperature data, which are to be actually included into the fitting processes, to temperatures up to 600 K (for both materials under study). Owing to these limitations we found that couples of just the first two anharmonicity-related expansion coefficients, and (in (5)), are sufficient for good simulations of the respective high-temperature curve sections. The present least-mean-square fitting processes are thus involving simultaneous adjustments of the (more or less strongly correlated) magnitudes of altogether 10 model-specific parameters.

4.3. Combined Fittings of Low- and High-Temperature Data and Results

Concerning the case of ZnSe (Figure 2) we have performed a simultaneous fit of the comprehensive set of data () given by Kremer et al. [46] in combination with a series of directly measured values (△) given for the region 400 K 600 K by Pashinkin and Malkova (cf. Table   in [78]). The resulting total set of adjusted parameters is listed in Table 1. For comparisons of the presently analyzed data with some former ones given by various authors we have still shown in Figure 2 the fine  K data points due to Birch [34] (+), the more or less strongly deviating data points due to Adachi et al. [10, 33] (◊) and Sirota et al. [79] (□), two data points given (for = 300 and 600 K) by Gadzhiev et al. [80] (●), and several smoothed values estimated by Pashinkin and Malkova [78] ().

For the case of ZnTe (Figure 3), the present fit comprises the set of data () given by Kremer et al. [46] in combination with the upper section (102 K < < 327 K) of the set of data points (□) given by Gavrichev et al. (in Table of [38]), including a set of approximate data points () due to Gadzhiev et al. (which we have redigitalized from Figure  1 shown in [80]), and a series of equidistant (smoothed) data points () available from the SGTE data review [81]. The resulting total set of adjusted parameters is listed in Table 1. For comparisons of the presently analyzed data with some former ones given by various authors we have still shown in Figure 3 the data points due to Adachi et al. [10, 33] (◊), Demidenko and Maltsev [31] (+), and the compatible sections (up to 600 K) of a series of smoothed values (□) quoted by Gavrichev et al. [82].

The fitted (continuous) and dependencies, which are implied by the material-specific sets of fitted parameter values (listed for ZnSe and ZnTe in the upper part of Table 1), are represented by solid and dashed curves, respectively, in Figures 2 and 3. In the insets of these figures we have represented the respective magnitudes of ratios. The series of empty circles [46] as well as the corresponding solid curves, , show pronounced non-Debye behaviours (maxima) in the vicinities of = 16.9 K (for ZnSe) and = 13.6 K (for ZnTe). The levels indicate the estimated limiting values, , of the curves, which are due exclusively to the lattice contributions to the measured heat capacities; that is, (within the liquid-helium-hydrogen region).

In the lower part of Table 1 we have still quoted a series of characteristic, material-specific quantities, the actual values of which are resulting in unambiguous way from the fitted parameter sets. This concerns in particular the associated magnitudes (due to (19)) of the coefficients and (referring to the conventional cryogenic odd-order series expansion (6)), the first- and second-order moments of the PDOS model function, and (due to (13)), the respective average phonon temperatures [14, 21, 50], , the dispersion coefficients [14, 21, 50], , and the high-temperature () limiting values of the true Debye temperatures [14, 21, 50], (cf. (42) below). Furthermore, with respect to the frequently considered (thermochemical) reference temperature, 298.15 K, we have still quoted the fitted (smoothed) isobaric heat capacity values, [14, 50], and the respective entropy values and enthalpy differences defined by [14, 15, 50]

We have given in Table 2, for a series of selected -values, a list of smoothed isobaric lattice heat capacity values, , including the corresponding (effective) lattice Debye temperatures, , which have been calculated on the basis of precision formulas to be presented below (in Section 5).

5. Transformation of Heat Capacities into Debye Temperatures

Expressive visualizations of typical non-Debye features of the temperature dependencies of heat capacities of solids are based above all on the results of point-by-point transformations of isobaric heat capacity data, , into the respective (effective) Debye temperature values, . The latter are well known from Debye’s classical paper [1] to be defined in implicit way [14, 1315, 56] by an integral of the form (1), where the magnitudes of the characteristic (dimensionless) variable, , are representing the upper limits of integration. The original task of calculating temperature dependencies of heat capacities within the frame of Debye’s theory [1] involved thus, primarily, the necessity of preparing good approximations for the dependence of Debye’s model-specific heat capacity function, , on . A detailed analytical study of the corresponding function (see (A.1)) is performed in Appendix A. Crucial results of this study are given, above all, in terms of approximate algebraic expressions for the respective low- and high-capacity dependencies, and , which can be, fortunately, transformed in exact way into corresponding algebraic expressions for the respective inverse dependencies, (see (B.1)) and (see (B.3)). Starting from these approximate expressions for low- and high-capacity curve sections, we succeeded to construct in Appendix B an unprecedented, integral-free analytical apparatus for approximate (more or less rough) estimations up to high-precision calculations of the complete dependence. Let us display in explicit form the integral-free formulas which are actually to be used here for the transformations of the heat capacity data under study (for ZnSe and ZnTe) into the respective Debye temperatures.

5.1. Asymptotic Low- and High-Temperature Formulas

The asymptote of Debye’s low-temperature expansion [1, 3], (cf. (A.4)), was already found in [13] to be of primary importance for a proper interpretation of the behaviour in the liquid-helium-hydrogen region. From the latter followed readily for the inverse function, , an asymptotic dependence of the form (see (B.1)). The respective limiting () dependence of the Debye temperature, , was thus explicitly given by the asymptotic expression [13, 15] (cf. Eq. in [13]), where the last representation of the dependence in (25) is in accordance with the general relation (cf. Eq. (6) in [13]). However, the inherent deviations (overestimations) of approximate values (due to (25)) from exact values are found to be lower than 1% only for very low magnitudes of heat capacities, 0 < 0.03, which are corresponding to values higher than 14 (cf. Table 9 and Figure 9).

A good basis for establishing an appropriate algebraic formula for the qualitatively different behaviour at considerably smaller values is provided by the unprecedented Taylor series expansion (see (B.3)), which we have shown here to follow in unambiguous way for the reciprocal version, (see (A.11) and (B.2)), of the function (see (A.10)). Accordingly, the temperature dependence of the Debye temperature, , for moderate-to-high heat capacities, is clearly confirmed to be given by the approximate algebraic expression [15]

Note that the latter formula corresponds just to Eq. () in [15], where it had already been taken into consideration for data discussions of various III–V materials. This formula is seen to provide a reasonable approximation (within possible deviations up to ±0.7%) to the exact curve throughout the range 0.05 < 1 (as shown in the upper inset of Figure 9 by the relative deviations, , of the approximate curve from the exact dependence).

5.2. Couple of High-Precision Formulas for Effective Debye Temperatures

In view of the sample character of the numerical analyses performed within the present study for the unusually fine and comprehensive data provided by Kremer et al. [46], we have constructed in Appendix B (Section B.3) a high-precision framework for the dependence. This consists in a combination of a highly elaborate (12-parameter) low-capacity formula, (see (B.8)), with a duly modified (11-parameter) high-capacity formula, (see (B.9)) (where = 0.26 represents the position of the point of crossing between the two complementary branches). The respective low- and high-temperature formulas for the Debye temperature, , are thus given explicitly by the couple of analytical expressions of the formfor the interval (with expansion coefficients , , and listed in the second column of Table 8), andfor the interval (with expansion coefficients listed in the third column of Table 8). Note that we have still endowed the functions in (27) and (28) by an additional subscript “” (included in brackets) in order to indicate that precisely the same formulas can also be used, alternatively, for high-accuracy calculations of the “true” Debye temperatures, , which are visualizing the characteristic non-Debye features of the harmonic parts of lattice heat capacity curves, (see (4); Figures 4 and 5).

By means of this couple of the high-accuracy formulas (27) and (28) we have performed, first of all, the point-by-point transformations, , for all those isobaric heat capacity data points, , with respect to which the concept of effective Debye temperatures is actually applicable, that is, for (cf. Figures 2 and 3). The respective Debye temperature values, > 0, are represented for ZnSe (in Figure 4) and for ZnTe (in Figure 5) by the same symbols as the original data points (shown in Figures 2 and 3, resp.).

Furthermore, we have transformed in the same way (i.e., via (27) or (28)) the continuous curves (solid curves shown in Figures 2 and 3), which resulted from the preceding numerical fittings in Section 4. These continuous curves (solid curves, in Figures 4 and 5) are naturally ending at those temperature points, , where the fitted curves (in Figure 2 or 3) are crossing just the classical Delong-Petit heat capacity value, , that is, at about 401 K for ZnSe (cf. Figures 2 and 4) and about 323 K for ZnTe (cf. Figures 3 and 5).

5.3. Calculation of “True” (Harmonic Lattice) Debye Temperatures

Of considerable interest, particularly from theoretical points of view [3, 14, 21, 28, 29, 32, 36, 56, 72], is the determination of the “true” (-related) Debye temperature curves, . The latter are defined with respect to the harmonic (isochoric) lattice heat capacities, (see (4)), by [14, 21]where the dimensionless variables are representing the respective upper boundaries of integration. Comparing (29) with (A.1) we see that the analytical connection of the upper boundaries of integration, , with the respective normalized (harmonic) heat capacity ratios, (in (29)), is just the same as the preceding connection between and (in (A.1)). Consequently, all the approximation formulas derived in Appendix B, for low- and high-capacity branches of the function (cf. Figure 9), can be readily used for calculations of the dependence implied by (29). This means, in particular, that (27) and (28) can be simultaneously used (as already indicated by the inclusion of the subscript “h”) for high-precision calculations of the dependencies of the “true” Debye temperatures, , on the respective harmonic (isochoric) heat capacities, . The corresponding material-specific curves are represented by dashed curves (for ZnSe and ZnTe in Figures 4 and 5, resp.).

With respect to this alternative use of (27) and (28), for calculations of curves, we would still like to point out that, apart from the liquid-helium-hydrogen region (see below), appreciable differences between and values, similarly to the differences between the underlying and values (cf. Section 4), are encountered as a rule only in ranges of moderately low to higher temperatures, , where the heat capacities in question are higher than about 50% of the classical Dulong-Petit limiting value, . This region is automatically comprised by the range of validity of (28). Thus it is mostly not necessary to involve (27) into separate calculations of curve sections.

On the other hand, significant differences between and values may also occur at very low (liquid-helium-hydrogen) temperatures, when we are (accidentally) concerned with nonnegligible contributions, , of a degenerate electronic system to the measured heat capacities, (cf. (5)). Consequently, in order to determine the respective temperature dependence of the curve within the liquid-helium-hydrogen region, we have to insert , instead of , into the respective asymptotic (low-) expression (see (25)). The monotonically decreasing sections, , of the “true” Debye temperature curves at sufficiently low temperatures, (see the insets of Figures 4 and 5), are thus adequately described by a corresponding asymptotic expression of the form [13]where and (see also [13] and cf. the insets of Figures 2 and 3). These asymptotic dependencies are represented by dotted curve sections in the insets of Figures 4 and 5.

6. Discussion

The present investigations were devoted above all to the finding of practicable solutions for longstanding computational problems which are frequently emerging within interpretations of the results of experimental heat capacity studies utilizing the somewhat troublesome concept of Debye temperatures [1, 3, 4, 6, 1015]. The notorious computational complications are due to circumstance that the respective values are involved into the definition of the upper limits of integration, , due to the integral expression for isobaric heat capacities, (see (1)). Consequently, in view of the previous lack of integral-free formulas for dependencies on values, it was thus necessary either to calculate the individual values via numerical integration procedures (using (1)) for all the individual points in consideration or to determine the values approximately, via point-by-point interpolations, using some of the available tables [6, 8, 4749].

Within our detailed study of the analytical properties of the conventional dependence (in Appendix A) we succeeded to derive, first of all, an unprecedented, rapidly converging high-temperature Taylor series expansion, (see (A.11)). Combining the latter with Debye’s conventional low-temperature expansion [1], (see (A.4)), it became possible to perform high-accuracy calculations of values (up to at least 10 significant figures; cf. Table 7) without involving any numerical integration procedure. The subsequent adoption of respective and formulas to the Debye-like component, (see (10)), of the long-wave TA phonon contribution to the cryogenic heat capacity, including the derivation of analogous formulas for the associated non-Debye component, (see (10)), enabled us to display (in Section 3.2) an integral-free version of the multioscillator hybrid model [14, 21, 50]. This has been used here (in Section 4) for numerical fittings of the data sets under study.

The high-accuracy results obtained in Appendix A for the dependence provided us, among other things, with the numerical basis for the solution of another and practically even more important mathematical problem, namely, the construction of appropriate integral-free formulas for the inverse Debye function dependence, that is, the function (see Figure 9 and Table 9), which describes the functional dependence of the characteristic Debye versus lattice temperature ratios, , on the ratios, , of given heat capacities versus Dulong-Petit classical heat capacity limit.

We have derived (in Appendix B) a series of more or less elaborate formulas for low- and high-capacity curve sections of the curve (solid curve in Figure 9). The simplest one of these approximate expressions is given in form of a junction formula, (see (B.4a)/(B.4b)), the actual use of which, however, is still connected with deviations from the exact dependence up to 0.6% (cf. the upper inset of Figure 9). Reductions of residual deviations to ±0.13%, or even to ±0.0153%, could be achieved (in Sections B.1 and B.2) via least-mean-square fittings of the exact dependence, adopting either the still relatively simple (2-parameter) formula, (see (B.5a) and (B.5b)), or the refined (9-parameter) formula, (see (B.6)) (cf. the upper and middle insets of Figure 9, resp.). Furthermore, for the sake of high-precision calculations of Debye temperatures, we have still constructed (in Section B.3) an elaborate couple of formulas consisting of a significantly refined low-capacity (12-parameter) formula, (see (B.8)), in combination with a similarly accurate high-capacity (11-parameter) formula, (see (B.9)). From the lower inset of Figure 9 we see that this couple of high-precision formulas effectuates an enormous reduction of residual deviations from the exact dependence by several orders of magnitude (down to ±0.000016%). A set of corresponding high-precision values (consisting of 7 significant figures) is given, for a series of equidistant -points, in Table 9.

Another important aim of the present study was an exhibition of the significant progress (in comparison with [14]) in the analytical and numerical description of the heat capacity properties of ZnSe and ZnTe, which we could achieve here owing to the availability of the unusually comprehensive and informative data sets given by Kremer et al. [46]. A remarkable improvement for the case of ZnSe is due, above all, to the dense set of novel data points () for the region 2 K 277 K, which can be seen from Figure 2 to differ significantly, for  K, from the relatively scarce (and obviously less accurate) set of former data points (□) due to [79], which had been considered previously in [14]. The qualitative difference between the novel and the former data sets for ZnSe is especially apparent from the visualization of the respective Debye temperature points, , in Figure 4. As a consequence of this qualitative difference we note, above all, that the present numerical analyses (in Section 4) led to a limiting value for the “true” (harmonic lattice) Debye temperature of = 317.1 K (cf. Table 1 and Figure 4), which is by 1.5% lower than the former estimate (of about 322 K; cf. [14])

Another significant progress is brought about by the unusually dense, informative sets of data points available from [46] for the liquid-helium-hydrogen region . These novel cryogenic data were of great importance particularly for the case of ZnTe, because in experimental papers published up to 2011 no data at all had been given for temperatures lower than 13 K (and those quoted in Table of [38] for the interval were obviously not sufficiently accurate; cf. the insets of Figures 3 and 5). This was the reason why it was practically impossible, prior the year 2012, to perform a direct determination of the caloric value for the limiting Debye temperature value, (cf. [13]).

Just with respect to the latter, however, we are now concerned, for ZnTe, with a significant difference (of about 4.5%) between the present estimate of 220 K (see the inset of Figure 5) and a former estimate of 230 K due to Kremer et al. [46]. Even a larger difference (of about 7%) is an analogous discrepancy, for ZnSe, between the present value of 270 K (see the inset of Figure 4) and a former estimate of 289 K due to Kremer et al. [46]. Let us examine in due detail the obvious causes of the discrepancies between values determined for ZnSe and ZnTe within the present study versus the former estimates quoted in [46].

6.1. Limiting Debye Temperatures and Electronic Heat Capacity Contributions

A provisional estimation of values can be made in simple way by viewing the sets of cryogenic Debye temperature points, (empty circles visualized in the insets of Figures 4, 5, and 6), which are exactly corresponding to the experimental data points provided for ZnSe and ZnTe by Kremer et al. [46] (as shown above by empty circles in the insets of Figures 2 and 3). Concerning the lower part of the liquid-helium-hydrogen region (2 K < < 10 K), we see readily that the respective values are throughout ranging below the following upper boundaries:

It is remarkable that these strict limitations refer even to the lowest experimental points (which are located in the vicinity of 2 K, in both cases).

In accordance with these mere visual estimations, it followed also from our detailed analyses (in Section 4) that the limiting ( 0) values of Debye temperatures, (cf. Table 1), are nearly coinciding (within possible uncertainties of order ±2 K) with these upper boundaries (31); that is,

These results are, of course, in clear contrast to the former estimates by Kremer et al. [46] (of about 289 K for ZnSe and 230 K for ZnTe). Let us reveal the basic cause of these discrepancies.

It is well known from various earlier studies of low-temperature heat capacities of nonmetals that, within relatively small temperature regions (from 0 up to about ), the isobaric heat capacities can as a rule be well approximated by odd-order power series expansions of the conventional form [13, 16, 18, 23, 35, 39, 65]where the term [4, 6, 13, 23, 6063] represents a possible electronic contribution and the subsequent three odd-order terms are due to harmonic lattice contributions. Of special interest in this connection is the lowest-order lattice contribution, , which corresponds to Debye’s classical ( limiting) cubic law [1, 3, 6, 8, 21], . The corresponding proportionality factor, , is known to be connected in a simple way with the limiting Debye temperature, , by the relation [4, 6, 13, 15, 46, 59, 62]

Consequently, within the frame of the frequently adopted versus representation [4, 6, 13, 15, 46, 62, 63, 83]the coefficient corresponds just to the slopes of the versus curves in the limit (which are represented in Figure 6 by dash-dot lines). However, these asymptotic lines used to be nearly coincident with the actual versus curves only within very narrow intervals, from 0 to about [16, 18, 35, 83, 84] (i.e., here up to only about 3 K for ZnSe and 2 K for ZnTe). Consequently, in view of the complete lack of data for the interval from 0 to 2 K, there remains only the possibility of determining in combination with the other expansion coefficients (i.e., , , and ) by careful fittings via (35) of the somewhat larger versus curve sections extending from 2 K up to about (see the solid curves in Figure 6). The parameter values obtained in this way are listed in Table 3. Viewing the respective values (listed in the last line of Table 3) we satisfy ourselves that, within the existing uncertainty of 2 K (as indicated above, in (32)), they are in accordance with the results of the analyses performed in Section 4 (see Table 1 and cf. Figures 4 and 5).

On the background of the curves shown in Figure 6 one can now easily see why the rougher simulations of data (in [46]) via a shortened series expansion of typefor up to 10 K, were a priori incapable of giving adequate results. We see that a shortened expansion of type (36) (due to taking = 0, in (35)) corresponds to the dash-double-point curves (shown in Figure 6), which begin to deviate significantly from the exact (solid) curves already in the vicinity of 6 K, for ZnSe, or 4 K, for ZnTe (i.e., at temperatures of order ). Thus one comes unavoidably to inadequate parameter values when a shortened (two-parameter) series expansion like (36) is carelessly applied to a markedly larger interval (up to 10 K [46]). In accordance with this observation we have still verified that, if necessary, even by means of a shortened series expansion like (36) it is, in principle, possible to detect physically reasonable values (closely approaching to those quoted in Table 3), provided the fitting ranges had been strictly limited to  K, for ZnSe, and  K, for ZnTe.

In summary one can assess that the obviously erroneous values quoted in [46] are the inevitable consequence of a careless application of an unnecessarily shortened series expansion (of type (36)) to disproportionately extended temperature intervals (up to 10 K [46]). At the same time we have found that, even within the frame of careful applications of the full-blown low-temperature series expansion (of type (35)), one is still concerned with residual uncertainties (of about 1%) for estimated values. These moderate uncertainties are obviously due to certain contributions of degenerate electronic systems to the measured heat capacities, , which found their numerical expression in nonvanishing magnitudes for the coefficients (cf. Tables 1 and 3). Concerning the actual magnitudes of the latter we see, among other things, that the respective value for ZnSe is by a factor of about 3 higher than its counterpart for ZnTe. Accordingly, the electronic contribution for ZnSe finds its graphical representations in form of a sufficiently well discernible nonmonotonic behaviour in the vicinity of 3.7 K, whereas the indications for a similar behaviour in the case of ZnTe, in the vicinity of 2.6 K, are markedly weaker (see the more or less pronounced local minimum behaviours of the data points [46] in the insets of Figures 2 and 3 and the respective local maximum behaviours of the corresponding data points in the insets of Figures 4, 5, and 6). Our analyses have shown that the presently available cryogenic data sets are not yet comprehensive enough for a more incisive quantification of the electronic contributions to the heat capacities of ZnSe and ZnTe. A significant refinement would require rather detailed heat capacity data to be available especially for the interval from 1 K to 2 K. The complete lack of such data is clearly the reason of the still existing uncertainties (of order 2 K.)

We satisfy ourselves, at least, that the present estimation of 270 K (2 K), including the particular value of = 271.4 K (Table 3), for ZnSe, is in good accordance with earlier caloric values of 271 K (1 K or 2 K) due to Birch [34] and Collins et al. [35]. Furthermore we can make the satisfactory statement that the present estimation of 220 K (2 K), including the particular value of = 219.3 K (Table 3), for ZnTe, is in good agreement with several values derived from elastic constants. This concerns in particular the values of 223.2 K (5 K) due to [35], 223.3 K due to [75], and 219.3 K due to [76]. On the other hand, one can also find in literature somewhat higher as well as markedly lower values for ZnTe, such as 225.3 K due to [74] and 228 K due to [77] or 190 K (±20 K) [38], 210 K due to [73], and 204.5 K due to [85]. Yet, in striking contrast to this multitude of earlier (mutually more or less compatible) estimates of values, for ZnTe, we have still found in contemporary literature a drastically deviating quotation of 758 K by Guo et al. [86], which appears to be classifiable as a physically obviously absurd one.

Concerning the and curves shown in the insets of Figure 6 we would still like to note that, by inserting the conventional odd-order power series expansions (see (33)) for the low-temperature dependence into the low-capacity formula (see (25)) for effective Debye temperatures, one comes readily to a relatively simple (algebraic) expression of type [13] in accordance with (34). This formula has been used here, with the parameter values quoted in Table 3, for calculations of the low-temperature sections of the dependencies (solid curves in the insets of Figure 6). The associated dependencies (due to (30)) for the “true” Debye temperatures, which are represented by dashed curves in the insets of Figure 6, have equally been calculated via a formula of type (37) where, merely, the electronic coefficient has been removed (). With respect to this characteristic behaviour it is necessary, of course, to point out that simple algebraic formulas [13] of type (37), for and , apply only to very narrow cryogenic intervals, 0 < < /30, where the respective Debye temperature curves show an exclusively convex shape (cf. the insets of Figure 6). Equivalently one can state that algebraic formula of the type (37) is not applicable to temperatures higher than /30 (i.e., to > 9 K for ZnSe or > 7 K ZnTe), where the (nearly coinciding) and dependencies are showing concave shapes (cf. the insets of Figures 4 and 5).

6.2. Minima of Debye Temperature Curves and Their Interpretation

The monotonic fall of the true Debye temperatures, , from their limiting ( 0) levels, , to the respective minima, , is described in good approximation by the relatively simple (algebraic) low-temperature formula, (see (30)). We see from the graphical representation of the latter by dotted curves (in the insets of Figures 4 and 5) that this expression provides, indeed, close approaches to the exact curves from 0 to . With respect to these limited cryogenic regions, , one can thus describe the monotonic decrease of the ratios in good approximation by the relationship [13, 14] in accordance with (30). The latter allows us to establish a useful approximate relationship between the minima, , of the curves, on the one hand, and the maxima, of the curves, on the other hand. Observing that the relative shifts of positions (indicated in the insets of Figures 4 and 5) versus the positions (indicated in the insets of Figures 2 and 3) are smaller than 2.5%, in both cases (cf. also Table 4), we can hence approximate the minima of the curves by . Accordingly, by specializing (38) for , we come to an approximate relationship between the extreme values of the versus curves of the form [14]

Comparing the actual and values (quoted in Table 4) we see that approximate relation (39) is indeed well fulfilled (within deviations smaller than 0.25%) by the empirical and values for ZnSe and ZnTe.

It is instructive to compare the presently obtained (empirical) values with their theoretical counterparts that had been estimated in various earlier papers. For the case of ZnSe we see that the previously estimated (theoretical) values of about 207 K (cf. Figure in [34]) and 218 K (cf. Figure  14 in [36]) are in reasonable agreement with the presently determined (empirical) value of 206.3 K (cf. Figure 4 and Table 4). In contrast to the latter, however, the theoretical curve presented in Figure of [87] shows an exclusively monotonic behaviour for any > 10 K, without giving any indication of a possible local minimum behaviour in the vicinity of 17 K (cf. the inset of Figure 4). It is suggested by Figure of [87] that the magnitude of the Debye temperature at 10 K should be of order 120 K, which is about 40% lower than the actual value. A detailed numerical comparison with Figure 4 in the present paper (including Table 2) reveals, consequently, that the curve shown for ZnSe in Figure   of [87] is quite inadequate within a rather large cryogenic interval, from 10 K up to about 80 K (note that an analogous criticism applies also to the obviously erroneous, monotonic curve shown for ZnS in Figure of the same paper).

For the case of ZnTe, we see that the previously estimated values of about 157 K (cf. Figure  5 in [32]), 162 K (cf. Figure  1 in [35]), 159 K (cf. Figure  15 in [36]), and 168 K (cf. Figure  3 in [38]) are throughout in reasonable agreement with the presently determined (empirical)