# Unprecedented Integral-Free Debye Temperature Formulas: Sample Applications to Heat Capacities of ZnSe and ZnTe

Authors: R. Pässler
Journal: Advances in Condensed Matter Physics
 (2017)

---

## Abstract

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,CVh(T), 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, Cp(T). 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 Cp(T) 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.


---

## Body

## 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 ΘD values (so-called “Debye temperatures”) within numerous thermophysical research papers, including various representative review articles [2–5] and books [6–10]. A brief inspection of the enormous ΘD data material published hitherto, however, showed readily that it is obviously not possible to find unique ΘD 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 (T-dependent) Debye temperatures, ΘD(T). Accordingly, Debye’s heat capacity model function [1] used to be employed in practice in the general form [3, 4, 6, 10–15](1)CpTCVh∞≡κT=3xDT3∫0xDTdxx4exex-12≤1,where the upper limits of integration are given by the ratios, xD(T)=ΘD(T)/T, of adjustable Debye temperatures versus the respective lattice temperatures and CVh(∞) represents the familiar Dulong-Petit limiting value for the isochoric lattice heat capacity in harmonic approximation (i.e., in particular CVh(∞)=6R, for binary materials).In the course of numerical fittings of measured (isobaric) heat capacities,Cp(T), 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 T-dependencies of the material-specific Debye temperatures, ΘD(T). This statement applies above all to the liquid-helium-hydrogen region, where the adjusted ΘD(T) values are as a rule rapidly falling [3–5, 12–17] from their T→0 limiting levels, ΘD(0), to certain minimum values, ΘDmin, which used to be located in the vicinities of T≈ΘD(0)/15. Furthermore, towards higher temperatures, many ΘD(T) curves show a more or less pronounced increase up to certain material-specific maxima, ΘDmax, the actual magnitudes of which are in many cases significantly higher than the respective ΘDmin levels. Such typical nonmonotonic behaviours of ΘD(T) curves in the cryogenic region have been found in particular for Si and Ge [18–21] as well as for a large variety of III–V [11, 12, 15, 22–30], II–VI [14, 31–38], 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 ΘD 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 ΘD 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 [41–44] for integrals of Debye’s type and imputing, at the same time, only some fixed (T-independent) values for Debye temperatures (like ΘD = 632 K for GaN [43] or ΘD = 920 K for ZnO [44] and ΘD = 946 K for MgO [44]).The basic cause of the occurrence of relatively deep minima,ΘDmin, 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 ofshort-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 gP(ε) 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], gD(ε)∝ε2, are continually confirmed to be indeed the main cause of the obviously typical, nonmonotonic (non-Debye) behaviours of ΘD(T) curves.Of course, when a representative set of properly determinedΘD(T) 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 x^D values quoted in some of the more or less comprehensive C(x^D) or C(x^D)/CVh(∞) tables presented by various authors [6, 8, 47–49] (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 ignoredinverse problem, namely, the reliable determination of effective Debye temperatures, ΘD(T), on the basis of measured heat capacities, Cp(T). 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 xD-ratios via interpolations between discrete x^D-ratios quoted for the couples of the two neighbouring Cp(T^) (or Cp(T^)/CVh(∞)) 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 Cp(T) data points. However, in contemporary literature, one is repeatedly concerned with rather large sets (of hundreds) of Cp(T) data points (like the presently analyzed Cp(T) 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 theinverse relationship, which should hence enable performing straightforward evaluations of the actual xD(T)=ΘD(T)/T-ratios directly on the basis of given heat capacity values, Cp(T). 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 (T<ΘD(T)/12),(2)ΘDT⟶ΘDl0T=T·4π4CVh∞5CpT1/3,and in [15] for the complementary range of higher temperatures (T>ΘD(T)/12),(3)ΘDT≈ΘDh0T=T·-35+5·49+56CVh∞CpT-122.However, due to a lack of space, the somewhat unusual derivation procedure for this highly useful high-temperature formula (i.e., Eq.(26) 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 Section2 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 Cp(T) and CVh(T) 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,CVh(T), are well known to be given by corresponding heat capacity shape functions, κh(T), of the general form [4, 6, 45, 50–55](4)CVhTCVh∞≡κhT=∫dεgPεε/2kBTsinhε/2kBT2≤1,where we have denoted by ε=ħω the energies of the individual phonons and gP(ε) 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, ∫dεgP(ε)=1.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 theoreticalCVh(T) functions, on the one hand, and those of measured (isobaric) heat capacities, Cp(T), on the other hand. The latter are well known to be throughout somewhathigher, Cp(T)>CVh(T), than their theoretical (isochoric) counterparts pertaining to theharmonic lattice regime [4, 50, 54]. The respective differences, Cp(T)-CVh(T)>0, are usually found to be very small from 0 up to temperatures of order Th, where the heat capacity amounts to about 50% of the Delong-Petit limiting value, Cp(Th)≅CVh(Th)=CVh(∞)/2. On the other hand, one is usually concerned with a relatively strong monotonic increase of these differences, d(Cp(T)-CVh(T))/dT>0, at higher temperatures (T>Th). 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 Th up to temperatures of order 5Th (at least), the rapid increase of these differences can be simulated in good approximation by a proportionality of type Cp(T)-CVh(T)∝T×(CVh(T))2∝T×(κh(T))2 (note that the structure of the latter is analogous to the known Nernst-Lindemann formula [2, 3, 58], Cp(T)-CV(T)∝T×(Cp(T))2, 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 T may come into play, it appeared reasonable to simulate the T-dependencies of isobaric heat capacities by a duly general algebraic expression of the form [14, 15, 21, 50, 54, 59](5)CpT=CVh∞κhT+κhT2∑n=1∞AnTn+c1T,where the empirical Taylor series expansion coefficients, An(n=1,2,…), are quantifying the material-specific weights of contributions of the individual, anharmonicity-related Tn power terms [15, 21, 54, 57]. Furthermore we have still included into (5) an additional linear term, C(E)(T)=c1T, which has to represent possible contributions of a degenerate electronic system [4, 6, 13, 23, 60–63] 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,gP(ε) (like those shown, e.g., in Figure  3 of [45] or in Figure  4 of [46], for ZnSe and ZnTe), which are determining the T-dependencies of the respective harmonic heat capacity shape functions, κh(T) (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 Cp(T) 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 variousprominent 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 ofodd-order terms [13, 16, 18, 23, 35, 39, 65](6)CpLT⟶CVhT=c3T3+c5T5+⋯.This limitingCVh(T) 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 exclusivelyeven-order power terms, gP(cont.)(ε)=γ2ε2+γ4ε4+⋯, of the phonon energy. An incorporation of a correspondingcontinuous low-energy section into theoretical gP(ε) 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., gP(cont.)(ε≥ε1)=0) at the positions of the first critical points [14, 21, 50, 64], ε1, 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., ε1=εTA(L) [7, 32, 36, 45, 46, 66, 67]).


### 3.1. General Analytical Framework of Multioscillator Hybrid Models

Less critical than the general ansatz for the0<ε<ε1 region is the modelling of gP(ε) spectra for energies higher than the first critical points, ε>ε1. Numerous analyses have shown that excellent simulations of CVh(T)∝κh(T) dependencies (4), from cryogenic to room temperatures, can be realized on the basis of combinations of several (macroscopic) Einstein oscillators [68], ∝wnδ(ε-εn)(n = 1 to nE), the ordered energy positions of which, ε1<ε2<⋯<εnE, including their actual weights, wn, 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](7)gPε=ϑε1-ε∑k=12wCk2k+1ε2kε12k+1+∑n=1nEwnδε-εn,where we have denoted by ϑ(ε1-ε) (=1 for ε<ε1 and 0 for ε>ε1) the step function, which realizes the cut of the continuous part of the gP(ε) 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 gP(ε) function to unity is equivalent to a normalization condition of the form [14, 21, 50](8)1=∫0→∞dεgPε=∑k=12wCk+∑n=1nEwn,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 gP(ε) model function (7).Inserting, finally, thegP(ε) function (see (7)) into the integral (4) and denoting by Θn≡εn/kB (n = 1 to nE) 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](9)κhT=∑k=12wCkκCkx1T+∑n=1nEwnΘn/2Tsinh⁡Θn/2T2,where the κCk(x1(T)) functions (being due to the two continuous components, ∝ε2k, in (7)) are given in form of integrals of Debye and non-Debye type [14, 21](10)κCkx1T=2k+1x1T2k+1∫0x1Tdxx2k+2exex-12.Here we have denoted by(11)x1T≡Θ1Tthe ratio of the characteristic phonon temperature associated with the first (lowest) special point, Θ1, versus lattice temperature, T (in analogy with the notation used within Debye’s theory, xD(T)≡ΘD/T; cf. (1) and Appendix A).Informative byproducts of numerical fittings of measuredCp(T) data sets via (5), in combinations with (9) for the harmonic heat capacity shape function, xh(T) (including (10), for the κCk(x1(T)) components), are themoments, μP(m), of material-specific PDOS spectral functions [3, 4, 16, 18, 19, 21, 50, 52, 56, 69–72], which are generally defined by(12)μPm≡∫dεgPεεm,for-3<m≠0.Inserting thegP(ε) hybrid model function (see (7)) into the integral for μP(m) (see (12)), we obtain for these moments the model-specific expression [14, 21](13)μPm=3wC13+m+5wC25+mε1m+∑n=1nEwnεnm.


### 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, κC1/C2(x1(T)) (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, κCk=1(x1(T)) (see Figure 1). A detailed analyses (in Section A.2 of Appendix A) has shown that the corresponding κC1(x1) 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, κC1lnD=1(x1≥xc1(D)) (see (A.14)) and κC1hmR=7(x1≤xc1(D)) (see (A.15)), assuming a switch between the two complementary branches just at their crossing point, xc1(D) = 6.5135 (cf. the upper inset of Figure 1). This means for the low-temperature branch, x1(T)≡Θ1/T≥xc1(D)=6.5135,(14)κC1x1⟶κC1lnD=1x1=4π45x13-3x1ex1-1-12e-x11+3x1+6x12+6x13,and for the high-temperature branch, 0≤x1(T)≡Θ1/T≤xc1(D)=6.5135,(15)κC1x1⟶κC1hmR=7x1=1+∑m=1mR=7c2mDRx12m-1 (the respective expansion coefficients of which, c2mDR, are listed in Table 6).Figure 1
Graphical representation of the exactx1-dependencies (—) of the normalized low-phonon-energy contributions to lattice heat capacities of Debye and non-Debye type, κC1(x1) and κC2(x1) (see (10)), in comparison with the respective low-temperature approximations (⋅–⋅–⋅–⋅), due to the functions κC1lnD=1(x1≥xc1(D)) (see (14)) or κC2lnN=1(x1≥xc2(D)) (see (16)), the associated high-temperature approximations (– – – – –), due to the functions κC1hmR=7(x1≤xc1(D)) (see (15)) or κC2hmR=6(x1≤xc2(D)) (see (17)), and the normalized heat capacity contribution of the first Einstein oscillator (– ⋅ ⋅ –), κE1(x1)=(x1/2)/sinh(x1/2)2.Concerning the integral (10) for the second low-energy component, κCk=2(x1(T)) (see Figure 1), we have found within the frame of a detailed study (in Section A.2, of Appendix A) that the corresponding κC2(x1) 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, κC2lnN=1(x1≥xc2(N)) (see (A.18)) and κC2hmR=6(x1≤xc2(N)) (see (A.23)), assuming again a switch between the two complementary branches at the respective crossing point, xc2(N)=7.0962 (cf. the lower inset of Figure 1). This means for the low-temperature branch, x1T≡Θ1/T≥xc2(N)=7.0962, explicitly(16)κC2x1⟶κC2lnN=1=80π621x15-5x1ex1-1-30e-x11+5x1+20x12+60x13+120x14+120x15,and for the high-temperature branch, 0≤x1(T)≡Θ1/T≤xc2(N)=7.0962,(17)κC2x1⟶κC2hmR=6x1=1+∑m=1mR=6c2mNRx12m-1 (the respective expansion coefficients of which, c2mNR, are listed in Table 6).The large scale representations of theκC1(x1) and κC2(x1) 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 κCk(x1) functions are really excellent only in the regions x1 < 5 and x1 > 9, whereas one is concerned in the vicinities of the crossing points, that is, at x1≈xc1D = 6.5135, for the Debye-like component, and at x1≈xc2(N) = 7.0962, for the non-Debye component, with certain overestimations of the exact κCk(x1) values by the approximate formulas (up to a maximum order of 0.008%, for k = 1, or of 0.015%, for k = 2). Yet, such computational inaccuracies are by 1 to 2 orders of magnitude smaller than the typical uncertainties of measured heat capacities, Cp(T). 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, wherex1(T)≡Θ1/T>12, the first terms (∝x1-3 or ∝x1-5) in the respective low-temperature expressions (see (14) and (16)) are strongly dominating. Thus the Debye-like and non-Debye components, κCk(x1(T)) ((see (10) for k = 1 and 2, resp.), tend here to the known ∝T3 and ∝T5 asymptotes [14, 21](18)κC1x1T⟶4π45·TΘ13,κC2x1T⟶80π621·TΘ15.The asymptote of thek=1 component is of Debye type [1], whereas the asymptote of the k=2 component is responsible for the onset of the typical non-Debye behaviour. The latter is frequently emerging already at very low temperatures (of order T≈ΘD(0)/50), where the Einstein oscillator terms in (9) for the κh(T) dependence are practically negligible (owing to their plateau behaviour, Θn/2Tsinh(Θn/2T)2→Θn/T2exp-Θn/T→0; cf. also the κE1(x1) curve in Figure 1). In this way we satisfy ourselves that the low-temperature dependence of the harmonic lattice heat capacity, CVh(T)=CVh(∞)κh(T) (see (4)), reduces, as expected, to a couple of two odd-order terms [13, 16, 18, 23, 35, 39, 65], CVh(T)→c3T3+c5T5 (see (6)). The respective expansion coefficients, c3 and c5, are readily seen to be given in terms of the weighting factors wC1 and wC2, the characteristic phonon temperature associated with the first Einstein peak, Θ1≡ε1/kB, and the classical Dulong-Petit heat capacity value, CVh(∞), by the relations(19)c3=wC14π45Θ13CVh∞,c5=wC280π627Θ15CVh∞.Comparing, finally, the analytical expression for the limiting (cubic) asymptote implied by the present formalism,CVh(T)→c3T3=wC14π4/5T/Θ13CVh(∞) (due to (19)), with its conventional predecessor, CVh(T)→4π4/5T/ΘD(0)3CVh(∞) (due to Debye’s theory [1–3, 8, 13]), we see that these two CVh(T→0) expressions are coincident when the weighting factor for the Debye-like component, wC1, is just given by the cube [14, 21] of the ratio between the characteristic phonon temperature associated with the first Einstein peak, Θ1, and the T→0 limit value, ΘD(0) of the conventional Debye temperature; that is,(20)wC1=Θ1ΘD03.


## 3.1. General Analytical Framework of Multioscillator Hybrid Models

Less critical than the general ansatz for the0<ε<ε1 region is the modelling of gP(ε) spectra for energies higher than the first critical points, ε>ε1. Numerous analyses have shown that excellent simulations of CVh(T)∝κh(T) dependencies (4), from cryogenic to room temperatures, can be realized on the basis of combinations of several (macroscopic) Einstein oscillators [68], ∝wnδ(ε-εn)(n = 1 to nE), the ordered energy positions of which, ε1<ε2<⋯<εnE, including their actual weights, wn, 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](7)gPε=ϑε1-ε∑k=12wCk2k+1ε2kε12k+1+∑n=1nEwnδε-εn,where we have denoted by ϑ(ε1-ε) (=1 for ε<ε1 and 0 for ε>ε1) the step function, which realizes the cut of the continuous part of the gP(ε) 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 gP(ε) function to unity is equivalent to a normalization condition of the form [14, 21, 50](8)1=∫0→∞dεgPε=∑k=12wCk+∑n=1nEwn,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 gP(ε) model function (7).Inserting, finally, thegP(ε) function (see (7)) into the integral (4) and denoting by Θn≡εn/kB (n = 1 to nE) 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](9)κhT=∑k=12wCkκCkx1T+∑n=1nEwnΘn/2Tsinh⁡Θn/2T2,where the κCk(x1(T)) functions (being due to the two continuous components, ∝ε2k, in (7)) are given in form of integrals of Debye and non-Debye type [14, 21](10)κCkx1T=2k+1x1T2k+1∫0x1Tdxx2k+2exex-12.Here we have denoted by(11)x1T≡Θ1Tthe ratio of the characteristic phonon temperature associated with the first (lowest) special point, Θ1, versus lattice temperature, T (in analogy with the notation used within Debye’s theory, xD(T)≡ΘD/T; cf. (1) and Appendix A).Informative byproducts of numerical fittings of measuredCp(T) data sets via (5), in combinations with (9) for the harmonic heat capacity shape function, xh(T) (including (10), for the κCk(x1(T)) components), are themoments, μP(m), of material-specific PDOS spectral functions [3, 4, 16, 18, 19, 21, 50, 52, 56, 69–72], which are generally defined by(12)μPm≡∫dεgPεεm,for-3<m≠0.Inserting thegP(ε) hybrid model function (see (7)) into the integral for μP(m) (see (12)), we obtain for these moments the model-specific expression [14, 21](13)μPm=3wC13+m+5wC25+mε1m+∑n=1nEwnεnm.


## 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, κC1/C2(x1(T)) (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, κCk=1(x1(T)) (see Figure 1). A detailed analyses (in Section A.2 of Appendix A) has shown that the corresponding κC1(x1) 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, κC1lnD=1(x1≥xc1(D)) (see (A.14)) and κC1hmR=7(x1≤xc1(D)) (see (A.15)), assuming a switch between the two complementary branches just at their crossing point, xc1(D) = 6.5135 (cf. the upper inset of Figure 1). This means for the low-temperature branch, x1(T)≡Θ1/T≥xc1(D)=6.5135,(14)κC1x1⟶κC1lnD=1x1=4π45x13-3x1ex1-1-12e-x11+3x1+6x12+6x13,and for the high-temperature branch, 0≤x1(T)≡Θ1/T≤xc1(D)=6.5135,(15)κC1x1⟶κC1hmR=7x1=1+∑m=1mR=7c2mDRx12m-1 (the respective expansion coefficients of which, c2mDR, are listed in Table 6).Figure 1
Graphical representation of the exactx1-dependencies (—) of the normalized low-phonon-energy contributions to lattice heat capacities of Debye and non-Debye type, κC1(x1) and κC2(x1) (see (10)), in comparison with the respective low-temperature approximations (⋅–⋅–⋅–⋅), due to the functions κC1lnD=1(x1≥xc1(D)) (see (14)) or κC2lnN=1(x1≥xc2(D)) (see (16)), the associated high-temperature approximations (– – – – –), due to the functions κC1hmR=7(x1≤xc1(D)) (see (15)) or κC2hmR=6(x1≤xc2(D)) (see (17)), and the normalized heat capacity contribution of the first Einstein oscillator (– ⋅ ⋅ –), κE1(x1)=(x1/2)/sinh(x1/2)2.Concerning the integral (10) for the second low-energy component, κCk=2(x1(T)) (see Figure 1), we have found within the frame of a detailed study (in Section A.2, of Appendix A) that the corresponding κC2(x1) 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, κC2lnN=1(x1≥xc2(N)) (see (A.18)) and κC2hmR=6(x1≤xc2(N)) (see (A.23)), assuming again a switch between the two complementary branches at the respective crossing point, xc2(N)=7.0962 (cf. the lower inset of Figure 1). This means for the low-temperature branch, x1T≡Θ1/T≥xc2(N)=7.0962, explicitly(16)κC2x1⟶κC2lnN=1=80π621x15-5x1ex1-1-30e-x11+5x1+20x12+60x13+120x14+120x15,and for the high-temperature branch, 0≤x1(T)≡Θ1/T≤xc2(N)=7.0962,(17)κC2x1⟶κC2hmR=6x1=1+∑m=1mR=6c2mNRx12m-1 (the respective expansion coefficients of which, c2mNR, are listed in Table 6).The large scale representations of theκC1(x1) and κC2(x1) 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 κCk(x1) functions are really excellent only in the regions x1 < 5 and x1 > 9, whereas one is concerned in the vicinities of the crossing points, that is, at x1≈xc1D = 6.5135, for the Debye-like component, and at x1≈xc2(N) = 7.0962, for the non-Debye component, with certain overestimations of the exact κCk(x1) values by the approximate formulas (up to a maximum order of 0.008%, for k = 1, or of 0.015%, for k = 2). Yet, such computational inaccuracies are by 1 to 2 orders of magnitude smaller than the typical uncertainties of measured heat capacities, Cp(T). 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, wherex1(T)≡Θ1/T>12, the first terms (∝x1-3 or ∝x1-5) in the respective low-temperature expressions (see (14) and (16)) are strongly dominating. Thus the Debye-like and non-Debye components, κCk(x1(T)) ((see (10) for k = 1 and 2, resp.), tend here to the known ∝T3 and ∝T5 asymptotes [14, 21](18)κC1x1T⟶4π45·TΘ13,κC2x1T⟶80π621·TΘ15.The asymptote of thek=1 component is of Debye type [1], whereas the asymptote of the k=2 component is responsible for the onset of the typical non-Debye behaviour. The latter is frequently emerging already at very low temperatures (of order T≈ΘD(0)/50), where the Einstein oscillator terms in (9) for the κh(T) dependence are practically negligible (owing to their plateau behaviour, Θn/2Tsinh(Θn/2T)2→Θn/T2exp-Θn/T→0; cf. also the κE1(x1) curve in Figure 1). In this way we satisfy ourselves that the low-temperature dependence of the harmonic lattice heat capacity, CVh(T)=CVh(∞)κh(T) (see (4)), reduces, as expected, to a couple of two odd-order terms [13, 16, 18, 23, 35, 39, 65], CVh(T)→c3T3+c5T5 (see (6)). The respective expansion coefficients, c3 and c5, are readily seen to be given in terms of the weighting factors wC1 and wC2, the characteristic phonon temperature associated with the first Einstein peak, Θ1≡ε1/kB, and the classical Dulong-Petit heat capacity value, CVh(∞), by the relations(19)c3=wC14π45Θ13CVh∞,c5=wC280π627Θ15CVh∞.Comparing, finally, the analytical expression for the limiting (cubic) asymptote implied by the present formalism,CVh(T)→c3T3=wC14π4/5T/Θ13CVh(∞) (due to (19)), with its conventional predecessor, CVh(T)→4π4/5T/ΘD(0)3CVh(∞) (due to Debye’s theory [1–3, 8, 13]), we see that these two CVh(T→0) expressions are coincident when the weighting factor for the Debye-like component, wC1, is just given by the cube [14, 21] of the ratio between the characteristic phonon temperature associated with the first Einstein peak, Θ1, and the T→0 limit value, ΘD(0) of the conventional Debye temperature; that is,(20)wC1=Θ1ΘD03.


## 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,nE) of Einstein oscillators to be taken into consideration. A relatively simple version in form of athree-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 offour 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 experimentalCp(T) 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 (εn+1 versus εn) energy positions. Considering under this aspect the known positions of the dominating optical (εLO versus εTO) 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 Cp(T) 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 asingle oscillator with the highest energy, ε4≡kBΘ4=kBΘTO&LO, and a fixed weight of(21)w4=wTO&LO=12.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, ε3≡kBΘ3=kBΘLA (= 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 of(22)w3=wLA=16,for the ε3-oscillator. Finally, with respect to the 2 TA phonon branches, we observe that the energy differences, for example, between the special points εTA(L) and εTA(X), 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, εn≡kBΘn=kBΘTAn (n = 1 and 2). Concerning their actual weights, w1 and w2, we observe that, owing to the preceding fixations of the weights for the LA and LO/TO oscillators (i.e., w4=1/2 (21) and w3=1/6 (22)), the global normalization condition (8) for the weights of all contributions associated with the two TA phonon branches (i.e., w1 and w2, in combination with wC1 and wC2) reduces, consequently, to a normalization condition of the TA phonon-related weighting factors of the special form [14](23)wC1+wC2+w1+w2=wTA=13.Accordingly, only 3 of these 4 weighting factors are actually playing the role ofindependent fitting parameters; that is, the magnitude of one of them is automatically fixed by the set of the three other ones (e.g., wC2=1/3-wC1-w1-w2). 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, κh(T) (9), reduces to a parameterized function of just 7 independent parameters (consisting of the 4 discrete phonon temperatures Θn,   n = 1 to 4, in combination with the three weighting factors w1, w2, and wC1).Furthermore we would still like to observe that, in many cases, one knows at least approximate values for theT→ 0 limiting values of Debye temperatures, ΘD(0) (either as results of former analyses of cryogenic heat capacity data [13, 34] or from independent estimations on the basis of elastic constants [35, 73–77]). This preliminary knowledge allows us, if necessary, to substitute the weight, wC1, of the Debye-like component, κC1(x1(T)), by wC1→Θ1/ΘD(0)3 (according to (20)).


### 4.2. Additional Empirical Parameters Involved by the Isobaric Heat Capacities

It has already been pointed out above (in Section2) that, for numerical simulations of measured isobaric heat capacity data, Cp(T), by means of (5), it is generally necessary to envisage the inclusion of further empirical parameters (like An, n=1,2,…, and/or c1) into the fitting process. Such additional parameters have been introduced above for a quantification of the deviations, Cp(T)-CVh(T)>0, of measured Cp(T) values from the theoretical (harmonic) lattice heat capacities, CVh(T)∝κh(T) (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 Cp(T) 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 Cplow(T) 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 Cplow(T) and Cplow(T)/T3 data points [46] represented by empty circles in Figures 2 and 3, resp., and cf. Figure 8 in [46]). Yet, particularly from the inset of Figure 2 we see that the sequence of the respective ρp(T)≡Cp(T)/T3 data points (○) for ZnSe are indicative of somenonmonotonic (local minimum) behaviour in the liquid-helium region. This is at variance in particular to earlier Cplow(T) data (+) given for the interval 1.7 K to 24.8 K by Birch (see Table 1 of [34]), which showed an exclusively monotonic behaviour of the sequence of the respective ρp(T)≡Cp(T)/T3 data points (+) in the same region.Figure 2
Fitting of the comprehensive set ofCplowT data (○) provided for ZnSe by Kremer et al. [46] in combination with a series of directly measured CphighT values (△) given for 400K≤T≤600K by Pashinkin and Malkova (Table 1 in [78]). For comparisons with other data available for ZnSe, we have also visualized the low-temperature data points given by Irwin and Lacombe [33] (◊), Birch [34] (+), and Sirota et al. [79] (□), including the smoothed high-temperature data points quoted by Gadzhiev et al. [80] (●) and Pashinkin and Malkova [78] (▲) [Cp(T) curve (— ; see (5)) and CVh(T) curve (– – – – –; see (4), in combination with (9) and (14) to (17)].Figure 3
Fitting of the set ofCplow(T) data (○) provided for ZnTe by Kremer et al. [46] in combination with the upper section (102 K < T < 327 K) of the Cplow(T) data set given by Gavrichev et al. [38] (□), a set of approximate Cphigh(T) data points due to Gadzhiev et al. [80] (■) (which we have redigitalized from their Figure 1), and a series of equidistant (smoothed) Cphigh(T) data points (∗) available from the SGTE data review [81]. For comparisons with other data available for ZnTe we have also visualized the low-temperature data points given by Demidenko and Maltsev [31] (+) and Irwin and Lacombe [33] (◊), including the smoothed high-temperature data points quoted by Gavrichev et al. [82] (□) [Cp(T) curve (— ; see (5)) and CVh(T) curve (– – – – –; see (4), in combination with (9) and (14) to (17)].The more or less pronounced nonmonotonic behaviours of Kremer’sρp(T) 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 Cplow(T) 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, C(E)(T)=c1T, 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 Cp(T) 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 firsttwo anharmonicity-related expansion coefficients, A1 and A2 (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 (Figure2) we have performed a simultaneous fit of the comprehensive set of Cplow(T) data (○) given by Kremer et al. [46] in combination with a series of directly measured Cphigh(T) values (△) given for the region 400 K ≤T≤ 600 K by Pashinkin and Malkova (cf. Table  1 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 Cplow(T<25 K) data points due to Birch [34] (+), the more or less strongly deviating Cplow(T) data points due to Adachi et al. [10, 33] (◊) and Sirota et al. [79] (□), two Cphigh(T) data points given (for T = 300 and 600 K) by Gadzhiev et al. [80] (●), and several smoothed Cphigh(T) values estimated by Pashinkin and Malkova [78] (▲).Table 1
Sets of the adjusted model-specific parameter values, which have been obtained (for properly fixedΘD(0) values) via combined least-mean-square fittings of the presently considered sets of low- and high-temperature isobaric heat capacity data, Cp(T), for ZnSe (Figure 2) and for ZnTe (Figure 3), using (5) in combination with (9) and (14) to (17). In the lower part of Table 1 we have still quoted a series of characteristic, material-specific quantities (cf. the text), the actual values of which are resulting in unambiguous way from the fitted parameter values listed previously.
Material
ZnSe
ZnTe
T-ranges (K)
2 to 600
2 to 600
Θ
D
(
f
i
x
e
d
)
(
0
) (K)
270
220
Θ
1
(
T
A
) (K)
78.32
61.93
Θ
2
(
T
A
) (K)
109.4
88.01
Θ
3
(
L
A
) (K)
213.0
167.1
Θ
4
(
T
O
&
L
O
) (K)
314.4
263.4
w
1
0.0632
0.0698
w
2
0.2349
0.2331
A
1 (K−1)
0.539 × 10−4
0.751 × 10−4
A
2 (K−2)
0.600 × 10−7
1.011 × 10−7
c
1 (J⋅K−2mol−1)
0.120 × 10−3
0.045 × 10−3
wC1=(Θ1(TA)/ΘD(0))3
0.0244
0.0223
w
C
2
=
1
/
3
-
w
C
1
-
w
1
-
w
2
0.0108
0.0081
c
3 (J⋅K−4mol−1)
0.1975 × 10−3
0.3651 × 10−3
c
5 (J⋅K−6mol−1)
0.667 × 10−6
1.626 × 10−6
μ
P
(
1
) (meV)
19.43
16.05
μ
P
(
2
) ((meV)2)
448.2
308.0
Θ
P (K)
225.5
185.8
Δ
P
0.43
0.45
Θ
D
h
(
∞
) (K)
317.1
262.9
C
p
(
T
r
) (J/K)
48.19
49.48
S
p
(
T
r
) (J/K)
72.12
82.12
H
p
(
T
r
)
-
H
p
(
0
) (kJ)
10.188
10.980For the case of ZnTe (Figure3), the present fit comprises the set of Cplow(T) data (○) given by Kremer et al. [46] in combination with the upper section (102 K < T < 327 K) of the set of Cplow(T) data points (□) given by Gavrichev et al. (in Table 1 of [38]), including a set of approximate Cphigh(T) data points (■) due to Gadzhiev et al. (which we have redigitalized from Figure  1 shown in [80]), and a series of equidistant (smoothed) Cphigh(T) 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 Cplow(T) 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 Cphigh(T) values (□) quoted by Gavrichev et al. [82].The fitted (continuous)Cp(T) and CVh(T) 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 ρp/h(T)≡Cp/Vh(T)/T3 ratios. The series of empty circles [46] as well as the corresponding solid curves, ρp(T), show pronounced non-Debye behaviours (maxima) in the vicinities of Tmax = 16.9 K (for ZnSe) and Tmax = 13.6 K (for ZnTe). The c3 levels indicate the estimated T→0 limiting values, c3=ρh(T→0), of the ρh(T)≡CVh(T)/T3 curves, which are due exclusively to thelattice contributions to the measured heat capacities; that is, CVh(T)≅Cp(T)-c1T (within the liquid-helium-hydrogen region).In the lower part of Table1 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 c3 and c5 (referring to the conventional cryogenic odd-order series expansion (6)), the first- and second-order moments of the PDOS model function, μP(1) and μP(2) (due to (13)), the respective average phonon temperatures [14, 21, 50], ΘP≡μP(1)/kB, the dispersion coefficients [14, 21, 50], ΔP≡μP(2)/(μP(1))2-12, and the high-temperature (T→∞) limiting values of the true Debye temperatures [14, 21, 50], ΘDh(∞)≡(5/3)μP(2)2/kB (cf. (42) below). Furthermore, with respect to the frequently considered (thermochemical) reference temperature, Tr≡ 298.15 K, we have still quoted the fitted (smoothed) isobaric heat capacity values, Cp(Tr) [14, 50], and the respective entropy values and enthalpy differences defined by [14, 15, 50](24)SpTr=∫0TrdTCpTT,HpTr-Hp0=∫0TrdTCpT.We have given in Table2, for a series of selected T-values, a list of smoothed isobaric lattice heat capacity values, Cp(L)(T)=Cp(T)-c1T, including the corresponding (effective) lattice Debye temperatures, ΘD(L)(T), which have been calculated on the basis of precision formulas to be presented below (in Section 5).Table 2
Smoothed values of the isobariclattice heat capacities, Cp(L)(T)=Cp(T)-c1T, and of the corresponding (effective) Debye temperatures, ΘD(L)(T) (due to (27) and (28)).
ZnSe
ZnSe
ZnTe
ZnTe
T (K)
C
p
(
L
) (J⋅K−1 mol−1)
Θ
D
(
L
) (K)
C
p
(
L
) (J⋅K−1 mol−1)
Θ
D
(
L
) (K)
2
0.0015949
269.16
0.0029434
219.43
4
0.013255
265.76
0.025043
214.98
6
0.048985
257.84
0.10158
202.19
8
0.13589
244.67
0.31540
184.80
10
0.32260
229.26
0.75959
172.32
12
0.65875
216.85
1.44623
166.75
14
1.15330
209.86
2.3390
165.39
16
1.78730
207.06
3.3764
166.53
18
2.5396
206.75
4.4832
169.36
20
3.3793
208.10
5.6015
173.32
25
5.6132
216.48
8.3098
185.43
30
7.7928
228.48
10.903
197.66
35
9.9059
240.78
13.439
208.50
40
12.002
251.94
15.931
217.65
45
14.096
261.58
18.361
225.22
50
16.176
269.74
20.704
231.40
55
18.221
276.58
22.938
236.40
60
20.210
282.27
25.045
240.41
65
22.125
287.00
27.014
243.60
70
23.952
290.91
28.841
246.12
75
25.682
294.13
30.525
248.07
80
27.309
296.78
32.071
249.57
85
28.832
298.95
33.485
250.68
90
30.251
300.72
34.775
251.47
95
31.569
302.14
35.950
251.99
100
32.790
303.28
37.021
252.28
110
34.964
304.85
38.822
252.29
120
36.819
305.68
40.429
251.69
130
38.401
305.94
41.721
250.58
140
39.755
305.73
42.807
249.01
150
40.916
305.14
43.727
247.01
160
41.917
304.19
44.512
244.57
170
42.784
302.91
45.188
241.67
180
43.540
301.31
45.774
238.29
190
44.203
299.38
46.288
234.38
200
44.787
297.12
46.740
229.88
220
45.766
291.54
47.503
218.84
240
46.522
284.38
48.122
204.41
260
47.197
275.36
48.640
185.43
280
47.735
264.13
49.083
159.77
300
48.194
250.18
49.471
122.68
350
49.100
197.64
50.283
—
400
49.792
77.81
50.961
—
450
50.363
—
51.573
—
500
50.863
—
52.153
—
550
51.321
—
52.722
—
600
51.757
—
53.293
—


## 4.1. Four-Oscillator Constellation Chosen for ZnSe and ZnTe

In the course of numerous analyses of experimentalCp(T) 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 (εn+1 versus εn) energy positions. Considering under this aspect the known positions of the dominating optical (εLO versus εTO) 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 Cp(T) 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 asingle oscillator with the highest energy, ε4≡kBΘ4=kBΘTO&LO, and a fixed weight of(21)w4=wTO&LO=12.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, ε3≡kBΘ3=kBΘLA (= 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 of(22)w3=wLA=16,for the ε3-oscillator. Finally, with respect to the 2 TA phonon branches, we observe that the energy differences, for example, between the special points εTA(L) and εTA(X), 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, εn≡kBΘn=kBΘTAn (n = 1 and 2). Concerning their actual weights, w1 and w2, we observe that, owing to the preceding fixations of the weights for the LA and LO/TO oscillators (i.e., w4=1/2 (21) and w3=1/6 (22)), the global normalization condition (8) for the weights of all contributions associated with the two TA phonon branches (i.e., w1 and w2, in combination with wC1 and wC2) reduces, consequently, to a normalization condition of the TA phonon-related weighting factors of the special form [14](23)wC1+wC2+w1+w2=wTA=13.Accordingly, only 3 of these 4 weighting factors are actually playing the role ofindependent fitting parameters; that is, the magnitude of one of them is automatically fixed by the set of the three other ones (e.g., wC2=1/3-wC1-w1-w2). 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, κh(T) (9), reduces to a parameterized function of just 7 independent parameters (consisting of the 4 discrete phonon temperatures Θn,   n = 1 to 4, in combination with the three weighting factors w1, w2, and wC1).Furthermore we would still like to observe that, in many cases, one knows at least approximate values for theT→ 0 limiting values of Debye temperatures, ΘD(0) (either as results of former analyses of cryogenic heat capacity data [13, 34] or from independent estimations on the basis of elastic constants [35, 73–77]). This preliminary knowledge allows us, if necessary, to substitute the weight, wC1, of the Debye-like component, κC1(x1(T)), by wC1→Θ1/ΘD(0)3 (according to (20)).


## 4.2. Additional Empirical Parameters Involved by the Isobaric Heat Capacities

It has already been pointed out above (in Section2) that, for numerical simulations of measured isobaric heat capacity data, Cp(T), by means of (5), it is generally necessary to envisage the inclusion of further empirical parameters (like An, n=1,2,…, and/or c1) into the fitting process. Such additional parameters have been introduced above for a quantification of the deviations, Cp(T)-CVh(T)>0, of measured Cp(T) values from the theoretical (harmonic) lattice heat capacities, CVh(T)∝κh(T) (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 Cp(T) 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 Cplow(T) 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 Cplow(T) and Cplow(T)/T3 data points [46] represented by empty circles in Figures 2 and 3, resp., and cf. Figure 8 in [46]). Yet, particularly from the inset of Figure 2 we see that the sequence of the respective ρp(T)≡Cp(T)/T3 data points (○) for ZnSe are indicative of somenonmonotonic (local minimum) behaviour in the liquid-helium region. This is at variance in particular to earlier Cplow(T) data (+) given for the interval 1.7 K to 24.8 K by Birch (see Table 1 of [34]), which showed an exclusively monotonic behaviour of the sequence of the respective ρp(T)≡Cp(T)/T3 data points (+) in the same region.Figure 2
Fitting of the comprehensive set ofCplowT data (○) provided for ZnSe by Kremer et al. [46] in combination with a series of directly measured CphighT values (△) given for 400K≤T≤600K by Pashinkin and Malkova (Table 1 in [78]). For comparisons with other data available for ZnSe, we have also visualized the low-temperature data points given by Irwin and Lacombe [33] (◊), Birch [34] (+), and Sirota et al. [79] (□), including the smoothed high-temperature data points quoted by Gadzhiev et al. [80] (●) and Pashinkin and Malkova [78] (▲) [Cp(T) curve (— ; see (5)) and CVh(T) curve (– – – – –; see (4), in combination with (9) and (14) to (17)].Figure 3
Fitting of the set ofCplow(T) data (○) provided for ZnTe by Kremer et al. [46] in combination with the upper section (102 K < T < 327 K) of the Cplow(T) data set given by Gavrichev et al. [38] (□), a set of approximate Cphigh(T) data points due to Gadzhiev et al. [80] (■) (which we have redigitalized from their Figure 1), and a series of equidistant (smoothed) Cphigh(T) data points (∗) available from the SGTE data review [81]. For comparisons with other data available for ZnTe we have also visualized the low-temperature data points given by Demidenko and Maltsev [31] (+) and Irwin and Lacombe [33] (◊), including the smoothed high-temperature data points quoted by Gavrichev et al. [82] (□) [Cp(T) curve (— ; see (5)) and CVh(T) curve (– – – – –; see (4), in combination with (9) and (14) to (17)].The more or less pronounced nonmonotonic behaviours of Kremer’sρp(T) 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 Cplow(T) 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, C(E)(T)=c1T, 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 Cp(T) 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 firsttwo anharmonicity-related expansion coefficients, A1 and A2 (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 (Figure2) we have performed a simultaneous fit of the comprehensive set of Cplow(T) data (○) given by Kremer et al. [46] in combination with a series of directly measured Cphigh(T) values (△) given for the region 400 K ≤T≤ 600 K by Pashinkin and Malkova (cf. Table  1 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 Cplow(T<25 K) data points due to Birch [34] (+), the more or less strongly deviating Cplow(T) data points due to Adachi et al. [10, 33] (◊) and Sirota et al. [79] (□), two Cphigh(T) data points given (for T = 300 and 600 K) by Gadzhiev et al. [80] (●), and several smoothed Cphigh(T) values estimated by Pashinkin and Malkova [78] (▲).Table 1
Sets of the adjusted model-specific parameter values, which have been obtained (for properly fixedΘD(0) values) via combined least-mean-square fittings of the presently considered sets of low- and high-temperature isobaric heat capacity data, Cp(T), for ZnSe (Figure 2) and for ZnTe (Figure 3), using (5) in combination with (9) and (14) to (17). In the lower part of Table 1 we have still quoted a series of characteristic, material-specific quantities (cf. the text), the actual values of which are resulting in unambiguous way from the fitted parameter values listed previously.
Material
ZnSe
ZnTe
T-ranges (K)
2 to 600
2 to 600
Θ
D
(
f
i
x
e
d
)
(
0
) (K)
270
220
Θ
1
(
T
A
) (K)
78.32
61.93
Θ
2
(
T
A
) (K)
109.4
88.01
Θ
3
(
L
A
) (K)
213.0
167.1
Θ
4
(
T
O
&
L
O
) (K)
314.4
263.4
w
1
0.0632
0.0698
w
2
0.2349
0.2331
A
1 (K−1)
0.539 × 10−4
0.751 × 10−4
A
2 (K−2)
0.600 × 10−7
1.011 × 10−7
c
1 (J⋅K−2mol−1)
0.120 × 10−3
0.045 × 10−3
wC1=(Θ1(TA)/ΘD(0))3
0.0244
0.0223
w
C
2
=
1
/
3
-
w
C
1
-
w
1
-
w
2
0.0108
0.0081
c
3 (J⋅K−4mol−1)
0.1975 × 10−3
0.3651 × 10−3
c
5 (J⋅K−6mol−1)
0.667 × 10−6
1.626 × 10−6
μ
P
(
1
) (meV)
19.43
16.05
μ
P
(
2
) ((meV)2)
448.2
308.0
Θ
P (K)
225.5
185.8
Δ
P
0.43
0.45
Θ
D
h
(
∞
) (K)
317.1
262.9
C
p
(
T
r
) (J/K)
48.19
49.48
S
p
(
T
r
) (J/K)
72.12
82.12
H
p
(
T
r
)
-
H
p
(
0
) (kJ)
10.188
10.980For the case of ZnTe (Figure3), the present fit comprises the set of Cplow(T) data (○) given by Kremer et al. [46] in combination with the upper section (102 K < T < 327 K) of the set of Cplow(T) data points (□) given by Gavrichev et al. (in Table 1 of [38]), including a set of approximate Cphigh(T) data points (■) due to Gadzhiev et al. (which we have redigitalized from Figure  1 shown in [80]), and a series of equidistant (smoothed) Cphigh(T) 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 Cplow(T) 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 Cphigh(T) values (□) quoted by Gavrichev et al. [82].The fitted (continuous)Cp(T) and CVh(T) 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 ρp/h(T)≡Cp/Vh(T)/T3 ratios. The series of empty circles [46] as well as the corresponding solid curves, ρp(T), show pronounced non-Debye behaviours (maxima) in the vicinities of Tmax = 16.9 K (for ZnSe) and Tmax = 13.6 K (for ZnTe). The c3 levels indicate the estimated T→0 limiting values, c3=ρh(T→0), of the ρh(T)≡CVh(T)/T3 curves, which are due exclusively to thelattice contributions to the measured heat capacities; that is, CVh(T)≅Cp(T)-c1T (within the liquid-helium-hydrogen region).In the lower part of Table1 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 c3 and c5 (referring to the conventional cryogenic odd-order series expansion (6)), the first- and second-order moments of the PDOS model function, μP(1) and μP(2) (due to (13)), the respective average phonon temperatures [14, 21, 50], ΘP≡μP(1)/kB, the dispersion coefficients [14, 21, 50], ΔP≡μP(2)/(μP(1))2-12, and the high-temperature (T→∞) limiting values of the true Debye temperatures [14, 21, 50], ΘDh(∞)≡(5/3)μP(2)2/kB (cf. (42) below). Furthermore, with respect to the frequently considered (thermochemical) reference temperature, Tr≡ 298.15 K, we have still quoted the fitted (smoothed) isobaric heat capacity values, Cp(Tr) [14, 50], and the respective entropy values and enthalpy differences defined by [14, 15, 50](24)SpTr=∫0TrdTCpTT,HpTr-Hp0=∫0TrdTCpT.We have given in Table2, for a series of selected T-values, a list of smoothed isobaric lattice heat capacity values, Cp(L)(T)=Cp(T)-c1T, including the corresponding (effective) lattice Debye temperatures, ΘD(L)(T), which have been calculated on the basis of precision formulas to be presented below (in Section 5).Table 2
Smoothed values of the isobariclattice heat capacities, Cp(L)(T)=Cp(T)-c1T, and of the corresponding (effective) Debye temperatures, ΘD(L)(T) (due to (27) and (28)).
ZnSe
ZnSe
ZnTe
ZnTe
T (K)
C
p
(
L
) (J⋅K−1 mol−1)
Θ
D
(
L
) (K)
C
p
(
L
) (J⋅K−1 mol−1)
Θ
D
(
L
) (K)
2
0.0015949
269.16
0.0029434
219.43
4
0.013255
265.76
0.025043
214.98
6
0.048985
257.84
0.10158
202.19
8
0.13589
244.67
0.31540
184.80
10
0.32260
229.26
0.75959
172.32
12
0.65875
216.85
1.44623
166.75
14
1.15330
209.86
2.3390
165.39
16
1.78730
207.06
3.3764
166.53
18
2.5396
206.75
4.4832
169.36
20
3.3793
208.10
5.6015
173.32
25
5.6132
216.48
8.3098
185.43
30
7.7928
228.48
10.903
197.66
35
9.9059
240.78
13.439
208.50
40
12.002
251.94
15.931
217.65
45
14.096
261.58
18.361
225.22
50
16.176
269.74
20.704
231.40
55
18.221
276.58
22.938
236.40
60
20.210
282.27
25.045
240.41
65
22.125
287.00
27.014
243.60
70
23.952
290.91
28.841
246.12
75
25.682
294.13
30.525
248.07
80
27.309
296.78
32.071
249.57
85
28.832
298.95
33.485
250.68
90
30.251
300.72
34.775
251.47
95
31.569
302.14
35.950
251.99
100
32.790
303.28
37.021
252.28
110
34.964
304.85
38.822
252.29
120
36.819
305.68
40.429
251.69
130
38.401
305.94
41.721
250.58
140
39.755
305.73
42.807
249.01
150
40.916
305.14
43.727
247.01
160
41.917
304.19
44.512
244.57
170
42.784
302.91
45.188
241.67
180
43.540
301.31
45.774
238.29
190
44.203
299.38
46.288
234.38
200
44.787
297.12
46.740
229.88
220
45.766
291.54
47.503
218.84
240
46.522
284.38
48.122
204.41
260
47.197
275.36
48.640
185.43
280
47.735
264.13
49.083
159.77
300
48.194
250.18
49.471
122.68
350
49.100
197.64
50.283
—
400
49.792
77.81
50.961
—
450
50.363
—
51.573
—
500
50.863
—
52.153
—
550
51.321
—
52.722
—
600
51.757
—
53.293
—


## 5. Transformation of Heat Capacities into Debye Temperatures

Expressive visualizations of typicalnon-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, Cp(T), into the respective (effective) Debye temperature values, ΘD(T). The latter are well known from Debye’s classical paper [1] to be defined in implicit way [1–4, 13–15, 56] by an integral of the form (1), where the magnitudes of the characteristic (dimensionless) variable, xD(T)≡ΘD(T)/T, 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, CD(T)→CVh(∞)κD(T), on xD(T). A detailed analytical study of the corresponding κD(xD) 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, κDl(xD>14) and κDh(xD<12), which can be, fortunately, transformed inexact way into corresponding algebraic expressions for the respectiveinverse dependencies, xDl0(0<κ<0.03) (see (B.1)) and xDh0(0.05<κ≤1) (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 xD(0<κ≤1) 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ΘD(T) Formulas

ThexD→∞ asymptote of Debye’s low-temperature expansion [1, 3], κ(xD)→(4π4/5xD3) (cf. (A.4)), was already found in [13] to be of primary importance for a proper interpretation of the ΘD(T) behaviour in the liquid-helium-hydrogen region. From the latter followed readily for theinverse function, xD(κ), an asymptotic dependence of the form xD(κ→0)→xDl0(κ)=4π4/5κ3 (see (B.1)). The respective limiting (T→0) dependence of the Debye temperature, ΘD(T)≡T·xD(T), was thus explicitly given by the asymptotic expression [13, 15] (25)ΘDl0T=T·4π45κT1/3=T·4π4CVh∞5CpT1/3=T·ΘD0c3CpT1/3 (cf. Eq. (9) in [13]), where the last representation of the ΘDl0(T) dependence in (25) is in accordance with the general relation c3=(4/5)π4CVh(∞)/(ΘD(0))3 (cf. Eq. (6) in [13]). However, the inherent deviations (overestimations) of approximate ΘDl0(T) values (due to (25)) from exact ΘD(T) values are found to be lower than 1% only for very low magnitudes of heat capacities, 0 ≤κ(T)=Cp(T)/CVh(∞) < 0.03, which are corresponding to xD(T)=ΘD(T)/T 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 smallerxD values is provided by the unprecedented Taylor series expansion xDh0(κD) (see (B.3)), which we have shown here to follow in unambiguous way for the reciprocal version, κDhmR=2(xD)≡(κDRhmR=2(xD))-1 (see (A.11) and (B.2)), of the κDRhmR=2(xD) function (see (A.10)). Accordingly, the temperature dependence of the Debye temperature, ΘD(T)≡T·xD(T), for moderate-to-high heat capacities, is clearly confirmed to be given by the approximate algebraic expression [15](26)ΘDT≅ΘDh0T=T·-35+5·49+56CVh∞CpT-122.Note that the latter formula corresponds just to Eq. (26) in [15], where it had already been taken into consideration for ΘD(T) data discussions of various III–V materials. This ΘDh0(T) formula is seen to provide a reasonable approximation (within possible deviations up to ±0.7%) to the exact ΘD(T) curve throughout the range 0.05 < κ(T)=Cp(T)/CVh(∞)≤ 1 (as shown in the upper inset of Figure 9 by the relative deviations, δxD/xD, of the approximate xDh0(κ) curve from the exact xD(κ) 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 comprehensiveCplow(T) data provided by Kremer et al. [46], we have constructed in Appendix B (Section B.3) a high-precision framework for the xD(κ) dependence. This consists in a combination of a highly elaborate (12-parameter) low-capacity formula, xDl12(κ≤κc) (see (B.8)), with a duly modified (11-parameter) high-capacity formula, xDh11(κ≥κc) (see (B.9)) (where κc = 0.26 represents the position of the point of crossing between the two complementary xD(κ) branches). The respective low- and high-temperature formulas for the Debye temperature, ΘDl/h(T)=T·xDl12/h11(κ(T)), are thus given explicitly by the couple of analytical expressions of the form(27)ΘDhlT=T·4π45κhT1/3-∑n=210anκhTn/2-1-exp-pκhT∑n=23gnκhTn,for the interval 0≤Cp(T)/CVh(∞)=κ(h)(T)≤κc=0.26 (with expansion coefficients an, p, and gn listed in the second column of Table 8), and(28)ΘDhhT=T·-35+5·49+561/κhT-1221+∑n=212bn1-κhTn,for the interval 0.26=κc≤κ(h)(T)=Cp(Vh)(T)/CVh(∞)≤1 (with expansion coefficients bn listed in the third column of Table 8). Note that we have still endowed the κ(T) functions in (27) and (28) by an additional subscript “h” (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, ΘDh(T), which are visualizing the characteristic non-Debye features of theharmonic parts of lattice heat capacity curves, CVh(T) (see (4); Figures 4 and 5).Figure 4
Effective Debye temperature values,ΘD(T), as resulting for ZnSe from measured (isobaric) heat capacities, Cp(T) ([46] (○), [34] (+), [10] (◊), [79] (□), [80] (●), and [78] (▲)). The continuous curve for the true Debye temperature, ΘDh(T) (– – – – –), is due to (28) in combination with (9) for the isochoric (harmonic) heat capacity shape function, κh(T)≡CVh(T)/6R (according to the parameters listed in Table 1). The continuous curve for the effective Debye temperature, ΘD(T) (solid curve), is due to (27) and (28) in combination with (5) and (9) for the associated isobaric (anharmonic) heat capacity shape function, κ(T)≡Cp(T)/6R. The dotted curve (in the inset) shows the low-temperature dependence of the true Debye temperature, ΘDh(T), according to (30).Figure 5
Effective Debye temperature values,ΘD(T), as resulting for ZnTe from measured (isobaric) heat capacities, Cp(T) ([46] (○), [38] (□), [31] (+), [10] (◊), [80] (■), and [81] (∗)). The continuous curve for the true Debye temperature, ΘDh(T) (– – – – –), is due to (28) in combination with (9) for the isochoric (harmonic) heat capacity shape function, κh(T)≡CVh(T)/6R (according to the parameters listed in Table 1). The continuous curve for effective Debye temperature, ΘD(T) (solid curve), is due to (27) and (28) in combination with (5) and (9) for the associated isobaric (anharmonic) heat capacity shape function, κ(T)≡Cp(T)/6R. The dotted curve (in the inset) shows the low-temperature dependence of the true Debye temperature, ΘDh(T), according to (30).By means of this couple of the high-accuracy formulas (27) and (28) we have performed, first of all, the point-by-point transformations, Cp(T)⇒ΘD(T), for all those isobaric heat capacity data points, Cp(T), with respect to which the concept ofeffective Debye temperatures is actually applicable, that is, for Cp(T)<CVh(∞)=6R (cf. Figures 2 and 3). The respective Debye temperature values, ΘD(T)> 0, are represented for ZnSe (in Figure 4) and for ZnTe (in Figure 5) by the same symbols as the original Cp(T) 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 Cp(T) curves (solid curves shown in Figures 2 and 3), which resulted from the preceding numerical fittings in Section 4. These continuous ΘD(T)>0 curves (solid curves, in Figures 4 and 5) are naturally ending at those temperature points, Tf, where the fitted Cp(T) curves (in Figure 2 or 3) are crossing just the classical Delong-Petit heat capacity value, Cp(Tf)=CVh(∞)=6R, 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” (CVh(T)-related) Debye temperature curves, ΘDh(T). The latter are defined with respect to the harmonic (isochoric) lattice heat capacities, CVh(T) (see (4)), by [14, 21](29)CVhTCVh∞≡κhT=3xDhT3∫0xDhTdxx4exex-12≤1,where the dimensionless variables xDh(T)≡ΘDh(T)/T 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, xDh, with the respective normalized (harmonic) heat capacity ratios, κh (in (29)), is just the same as the preceding connection between xD and κ (in (A.1)). Consequently, all the approximation formulas derived in Appendix B, for low- and high-capacity branches of the xD(κ) function (cf. Figure 9), can be readily used for calculations of the xDh(κh) 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, ΘDh(T)=T·xDh(T), on the respective harmonic (isochoric) heat capacities, CVh(T)=κh(T)·CVh(∞). The corresponding material-specific ΘDh(T) 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 ΘDh(T) curves, we would still like to point out that, apart from the liquid-helium-hydrogen region (see below), appreciable differences between ΘDh(T) and ΘD(T) values, similarly to the differences between the underlying CVh(T) and Cp(T) values (cf. Section 4), are encountered as a rule only in ranges of moderately low to higher temperatures, T>Th, where the heat capacities in question are higher than about 50% of the classical Dulong-Petit limiting value, κh(T)=CVh(T)/CVh(∞)>0.5. 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 ΘDh(T) curve sections.On the other hand, significant differences betweenΘDh(T) and ΘD(T) values may also occur at very low (liquid-helium-hydrogen) temperatures, when we are (accidentally) concerned with nonnegligible contributions, C(E)(T)=c1T, of a degenerate electronic system to the measured heat capacities, Cp(T)→CVh(T)+c1T (cf. (5)). Consequently, in order to determine the respective temperature dependence of the ΘDh(T) curve within the liquid-helium-hydrogen region, we have to insert κh(T)=(Cp(T)-c1(T))/CVh(∞), instead of κ(T)=Cp(T)/CVh(∞), into the respective asymptotic (low-T) expression (see (25)). The monotonically decreasing sections, dΘDh(T)/dT<0, of the “true” Debye temperature curves at sufficiently low temperatures, T<Tmin (see the insets of Figures 4 and 5), are thus adequately described by a corresponding asymptotic expression of the form [13](30)ΘDhl0T=T·4π45κhT1/3=T·4π4CVh∞5CVhT1/3=ΘD0ρh0ρhT1/3,where ρh(T)≡CVh(T)/T3 and ρh(0)=c3=(4/5)π4CVh(∞)/(ΘD(0))3 (see also [13] and cf. the insets of Figures 2 and 3). These asymptotic ΘDhl0(T) dependencies are represented by dotted curve sections in the insets of Figures 4 and 5.


## 5.1. Asymptotic Low- and High-TemperatureΘD(T) Formulas

ThexD→∞ asymptote of Debye’s low-temperature expansion [1, 3], κ(xD)→(4π4/5xD3) (cf. (A.4)), was already found in [13] to be of primary importance for a proper interpretation of the ΘD(T) behaviour in the liquid-helium-hydrogen region. From the latter followed readily for theinverse function, xD(κ), an asymptotic dependence of the form xD(κ→0)→xDl0(κ)=4π4/5κ3 (see (B.1)). The respective limiting (T→0) dependence of the Debye temperature, ΘD(T)≡T·xD(T), was thus explicitly given by the asymptotic expression [13, 15] (25)ΘDl0T=T·4π45κT1/3=T·4π4CVh∞5CpT1/3=T·ΘD0c3CpT1/3 (cf. Eq. (9) in [13]), where the last representation of the ΘDl0(T) dependence in (25) is in accordance with the general relation c3=(4/5)π4CVh(∞)/(ΘD(0))3 (cf. Eq. (6) in [13]). However, the inherent deviations (overestimations) of approximate ΘDl0(T) values (due to (25)) from exact ΘD(T) values are found to be lower than 1% only for very low magnitudes of heat capacities, 0 ≤κ(T)=Cp(T)/CVh(∞) < 0.03, which are corresponding to xD(T)=ΘD(T)/T 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 smallerxD values is provided by the unprecedented Taylor series expansion xDh0(κD) (see (B.3)), which we have shown here to follow in unambiguous way for the reciprocal version, κDhmR=2(xD)≡(κDRhmR=2(xD))-1 (see (A.11) and (B.2)), of the κDRhmR=2(xD) function (see (A.10)). Accordingly, the temperature dependence of the Debye temperature, ΘD(T)≡T·xD(T), for moderate-to-high heat capacities, is clearly confirmed to be given by the approximate algebraic expression [15](26)ΘDT≅ΘDh0T=T·-35+5·49+56CVh∞CpT-122.Note that the latter formula corresponds just to Eq. (26) in [15], where it had already been taken into consideration for ΘD(T) data discussions of various III–V materials. This ΘDh0(T) formula is seen to provide a reasonable approximation (within possible deviations up to ±0.7%) to the exact ΘD(T) curve throughout the range 0.05 < κ(T)=Cp(T)/CVh(∞)≤ 1 (as shown in the upper inset of Figure 9 by the relative deviations, δxD/xD, of the approximate xDh0(κ) curve from the exact xD(κ) 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 comprehensiveCplow(T) data provided by Kremer et al. [46], we have constructed in Appendix B (Section B.3) a high-precision framework for the xD(κ) dependence. This consists in a combination of a highly elaborate (12-parameter) low-capacity formula, xDl12(κ≤κc) (see (B.8)), with a duly modified (11-parameter) high-capacity formula, xDh11(κ≥κc) (see (B.9)) (where κc = 0.26 represents the position of the point of crossing between the two complementary xD(κ) branches). The respective low- and high-temperature formulas for the Debye temperature, ΘDl/h(T)=T·xDl12/h11(κ(T)), are thus given explicitly by the couple of analytical expressions of the form(27)ΘDhlT=T·4π45κhT1/3-∑n=210anκhTn/2-1-exp-pκhT∑n=23gnκhTn,for the interval 0≤Cp(T)/CVh(∞)=κ(h)(T)≤κc=0.26 (with expansion coefficients an, p, and gn listed in the second column of Table 8), and(28)ΘDhhT=T·-35+5·49+561/κhT-1221+∑n=212bn1-κhTn,for the interval 0.26=κc≤κ(h)(T)=Cp(Vh)(T)/CVh(∞)≤1 (with expansion coefficients bn listed in the third column of Table 8). Note that we have still endowed the κ(T) functions in (27) and (28) by an additional subscript “h” (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, ΘDh(T), which are visualizing the characteristic non-Debye features of theharmonic parts of lattice heat capacity curves, CVh(T) (see (4); Figures 4 and 5).Figure 4
Effective Debye temperature values,ΘD(T), as resulting for ZnSe from measured (isobaric) heat capacities, Cp(T) ([46] (○), [34] (+), [10] (◊), [79] (□), [80] (●), and [78] (▲)). The continuous curve for the true Debye temperature, ΘDh(T) (– – – – –), is due to (28) in combination with (9) for the isochoric (harmonic) heat capacity shape function, κh(T)≡CVh(T)/6R (according to the parameters listed in Table 1). The continuous curve for the effective Debye temperature, ΘD(T) (solid curve), is due to (27) and (28) in combination with (5) and (9) for the associated isobaric (anharmonic) heat capacity shape function, κ(T)≡Cp(T)/6R. The dotted curve (in the inset) shows the low-temperature dependence of the true Debye temperature, ΘDh(T), according to (30).Figure 5
Effective Debye temperature values,ΘD(T), as resulting for ZnTe from measured (isobaric) heat capacities, Cp(T) ([46] (○), [38] (□), [31] (+), [10] (◊), [80] (■), and [81] (∗)). The continuous curve for the true Debye temperature, ΘDh(T) (– – – – –), is due to (28) in combination with (9) for the isochoric (harmonic) heat capacity shape function, κh(T)≡CVh(T)/6R (according to the parameters listed in Table 1). The continuous curve for effective Debye temperature, ΘD(T) (solid curve), is due to (27) and (28) in combination with (5) and (9) for the associated isobaric (anharmonic) heat capacity shape function, κ(T)≡Cp(T)/6R. The dotted curve (in the inset) shows the low-temperature dependence of the true Debye temperature, ΘDh(T), according to (30).By means of this couple of the high-accuracy formulas (27) and (28) we have performed, first of all, the point-by-point transformations, Cp(T)⇒ΘD(T), for all those isobaric heat capacity data points, Cp(T), with respect to which the concept ofeffective Debye temperatures is actually applicable, that is, for Cp(T)<CVh(∞)=6R (cf. Figures 2 and 3). The respective Debye temperature values, ΘD(T)> 0, are represented for ZnSe (in Figure 4) and for ZnTe (in Figure 5) by the same symbols as the original Cp(T) 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 Cp(T) curves (solid curves shown in Figures 2 and 3), which resulted from the preceding numerical fittings in Section 4. These continuous ΘD(T)>0 curves (solid curves, in Figures 4 and 5) are naturally ending at those temperature points, Tf, where the fitted Cp(T) curves (in Figure 2 or 3) are crossing just the classical Delong-Petit heat capacity value, Cp(Tf)=CVh(∞)=6R, 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” (CVh(T)-related) Debye temperature curves, ΘDh(T). The latter are defined with respect to the harmonic (isochoric) lattice heat capacities, CVh(T) (see (4)), by [14, 21](29)CVhTCVh∞≡κhT=3xDhT3∫0xDhTdxx4exex-12≤1,where the dimensionless variables xDh(T)≡ΘDh(T)/T 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, xDh, with the respective normalized (harmonic) heat capacity ratios, κh (in (29)), is just the same as the preceding connection between xD and κ (in (A.1)). Consequently, all the approximation formulas derived in Appendix B, for low- and high-capacity branches of the xD(κ) function (cf. Figure 9), can be readily used for calculations of the xDh(κh) 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, ΘDh(T)=T·xDh(T), on the respective harmonic (isochoric) heat capacities, CVh(T)=κh(T)·CVh(∞). The corresponding material-specific ΘDh(T) 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 ΘDh(T) curves, we would still like to point out that, apart from the liquid-helium-hydrogen region (see below), appreciable differences between ΘDh(T) and ΘD(T) values, similarly to the differences between the underlying CVh(T) and Cp(T) values (cf. Section 4), are encountered as a rule only in ranges of moderately low to higher temperatures, T>Th, where the heat capacities in question are higher than about 50% of the classical Dulong-Petit limiting value, κh(T)=CVh(T)/CVh(∞)>0.5. 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 ΘDh(T) curve sections.On the other hand, significant differences betweenΘDh(T) and ΘD(T) values may also occur at very low (liquid-helium-hydrogen) temperatures, when we are (accidentally) concerned with nonnegligible contributions, C(E)(T)=c1T, of a degenerate electronic system to the measured heat capacities, Cp(T)→CVh(T)+c1T (cf. (5)). Consequently, in order to determine the respective temperature dependence of the ΘDh(T) curve within the liquid-helium-hydrogen region, we have to insert κh(T)=(Cp(T)-c1(T))/CVh(∞), instead of κ(T)=Cp(T)/CVh(∞), into the respective asymptotic (low-T) expression (see (25)). The monotonically decreasing sections, dΘDh(T)/dT<0, of the “true” Debye temperature curves at sufficiently low temperatures, T<Tmin (see the insets of Figures 4 and 5), are thus adequately described by a corresponding asymptotic expression of the form [13](30)ΘDhl0T=T·4π45κhT1/3=T·4π4CVh∞5CVhT1/3=ΘD0ρh0ρhT1/3,where ρh(T)≡CVh(T)/T3 and ρh(0)=c3=(4/5)π4CVh(∞)/(ΘD(0))3 (see also [13] and cf. the insets of Figures 2 and 3). These asymptotic ΘDhl0(T) 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, 10–15]. The notorious computational complications are due to circumstance that the respective ΘD(T) values are involved into the definition of theupper limits of integration, xD(T)≡ΘD(T)/T, due to the integral expression for isobaric heat capacities, Cp(T) (see (1)). Consequently, in view of the previous lack of integral-free formulas for Cp(T) dependencies on xD(T) values, it was thus necessary either to calculate the individual Cp(T) values via numerical integration procedures (using (1)) for all the individual xD(T) points in consideration or to determine the Cp(T)=CVh(∞)κ(T) values approximately, via point-by-point interpolations, using some of the available κD(xD) tables [6, 8, 47–49].Within our detailed study of the analytical properties of the conventionalκD(xD) dependence (in Appendix A) we succeeded to derive, first of all, an unprecedented, rapidly converging high-temperature Taylor series expansion, κDh(mR=7)(xD) (see (A.11)). Combining the latter with Debye’s conventional low-temperature expansion [1], κDl(nD=15)(xD) (see (A.4)), it became possible to perform high-accuracy calculations of κD(xD) values (up to at least 10 significant figures; cf. Table 7) without involving any numerical integration procedure. The subsequent adoption of respective κDh(xD) and κDl(xD) formulas to the Debye-like component, κC1(x1(T)) (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, κC2(x1(T)) (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 Cp(T) data sets under study.The high-accuracy results obtained in AppendixA for the κ(xD) 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 theinverse Debye function dependence, that is, the xD(κ) function (see Figure 9 and Table 9), which describes the functional dependence of the characteristic Debye versus lattice temperature ratios, xD(T)=ΘD(T)/T, on the ratios, κ(T)=Cp(T)/CVh(∞), of given heat capacities versus Dulong-Petit classical heat capacity limit.We have derived (in AppendixB) a series of more or less elaborate formulas for low- and high-capacity curve sections of the xD(κ) curve (solid curve in Figure 9). The simplest one of these approximate xD(κ) expressions is given in form of a junction formula, x~D(κ) (see (B.4a)/(B.4b)), the actual use of which, however, is still connected with deviations from the exact xD(κ) 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 xD(κ) dependence, adopting either the still relatively simple (2-parameter) formula, xDl2(0<κ≤κ2) (see (B.5a) and (B.5b)), or the refined (9-parameter) formula, xDl9(0<κ≤κ9) (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, xDl12(0<κ≤0.26) (see (B.8)), in combination with a similarly accurate high-capacity (11-parameter) formula, xDh11(0.26≤κ≤1) (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 xD(κ) dependence by several orders of magnitude (down to ±0.000016%). A set of corresponding high-precision xD(κ) 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 Cp(T) 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 Cp(T) data points (○) for the region 2 K ≤T≤ 277 K, which can be seen from Figure 2 to differ significantly, for T>100 K, from the relatively scarce (and obviously less accurate) set of former Cp(T) data points (□) due to [79], which had been considered previously in [14]. The qualitative difference between the novel and the former Cp(T) data sets for ZnSe is especially apparent from the visualization of the respective Debye temperature points, ΘD(T), 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 T→∞ limiting value for the “true” (harmonic lattice) Debye temperature of ΘDh(∞) = 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 ofCplow(T) data points available from [46] for the liquid-helium-hydrogen region (2K≤T≤20K). These novel cryogenic data were of great importance particularly for the case of ZnTe, because in experimental papers published up to 2011 no Cp(T) data at all had been given for temperatures lower than 13 K (and those quoted in Table 1 of [38] for the interval 13K<T<18K 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 thecaloric value for the T→0 limiting Debye temperature value, ΘD(0) (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ΘD(0)≅ 220 K (see the inset of Figure 5) and a former estimate of ΘD(0)≈ 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 ΘD(0)≅ 270 K (see the inset of Figure 4) and a former estimate of ΘD(0)≈ 289 K due to Kremer et al. [46]. Let us examine in due detail the obvious causes of the discrepancies between ΘD(0) 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ΘD(0) values can be made in simple way by viewing the sets of cryogenic Debye temperature points, ΘD(T) (empty circles visualized in the insets of Figures 4, 5, and 6), which are exactly corresponding to the experimental Cp(T)/T3 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 < T < 10 K), we see readily that the respective ΘD(T) values are throughout rangingbelow the following upper boundaries:(31)ΘDT<270Kfor  ZnSe,ΘDT<220Kfor  ZnTe.Figure 6
Fittings via (35) of the Cp(T)/T data points (○) corresponding to the cryogenic Cp(T) data provided by Kremer et al. [46] for ZnSe and ZnTe [—, least-mean square fittings from 2 K to 8.7 K or to 6.4 K, resp.; ⋅⋅–⋅⋅–⋅⋅–⋅⋅, tail behaviour, for c7→0; ⋅–⋅–⋅–⋅, T→0 asymptote, for c7=0 and c5=0]. The insets show the effective ΘD(T) values (○) following directly (via (25)) from Kremer’s cryogenic Cp(T) data [46] for ZnSe and ZnTe. For comparison we have still included the ΘD(T) values (+) following from former Cp(T) data quoted for ZnSe by Birch [34], which show no indication for an appreciable electronic contribution. The continuous curves in the insets are resulting from (25), implying the cryogenic Cp(T) series expansion (see (33)), with fitted Taylor series expansion coefficients listed in Table 3. [—, ΘD(T) curve, for c1 > 0;   – – – – –,    ΘDh(T) curve, for c1→ 0.].Table 3
Conventional Taylor series expansion coefficients determined by least-mean square fittings of cryogenicCp(T)/T data [46] for ZnSe and ZnTe (see Figure 6) via (35). The respective ΘD(0) values are resulting from fitted c3 values via (34).
ZnSe
ZnTe
Intervals (K)
2.0 to 8.7
1.9 to 6.4
c
1 (mJ⋅K−2 mol−1)
0.1339
0.0409
c
3 (mJ⋅K−4 mol−1)
0.1937
0.3686
c
5 (mJ⋅K−6 mol−1)
0.602 × 10−3
0.537 × 10−3
c
7 (mJ⋅K−8 mol−1)
0.815 × 10−5
0.638 × 10−4
Θ
D
(
0
) (K)
271.7
219.3It 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 Section4) that the limiting (T→ 0) values of Debye temperatures, ΘD(0) (cf. Table 1), are nearly coinciding (within possible uncertainties of order ±2 K) with these upper boundaries (31); that is, (32)ΘD0≅270K±2K,for  ZnSe,ΘD0≅220K±2K,for  ZnTe.These results are, of course, in clear contrast to the formerΘD(0) 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ΘD(0)/30), 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](33)CpT=∑n=03c2n+1T2n+1,where the c1T term [4, 6, 13, 23, 60–63] 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, c3T3, which corresponds to Debye’s classical (T→0 limiting) cubic law [1, 3, 6, 8, 21], CVh(T→0)∝T3. The corresponding proportionality factor, c3, is known to be connected in a simple way with the limiting Debye temperature, ΘD(0), by the relation [4, 6, 13, 15, 46, 59, 62](34)ΘD0=4π4CVh∞5c33.Consequently, within the frame of the frequently adoptedCp(T)/T versus T2 representation [4, 6, 13, 15, 46, 62, 63, 83](35)CpTT=c1+c3T2+c5T4+c7T6,the coefficient c3 corresponds just to the slopes of the Cp(T)/T versus T2 curves in the T→0 limit (which are represented in Figure 6 by dash-dot lines). However, these asymptotic lines used to be nearly coincident with the actual Cp(T)/T versus T2 curves only within very narrow intervals, from 0 to about ΘD(0)/100 [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 Cp(T) data for the interval from 0 to 2 K, there remains only the possibility of determining c3 in combination with the other expansion coefficients (i.e., c1, c5, and c7) by careful fittings via (35) of the somewhat larger Cp(T)/T versus T2 curve sections extending from 2 K up to about ΘD(0)/30 (see the solid curves in Figure 6). The parameter values obtained in this way are listed in Table 3. Viewing the respective ΘD(0) 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 Figure6 one can now easily see why the rougher simulations of Cp(T)/T data (in [46]) via a shortened series expansion of type(36)CpTT≅c3T2+c5T4,for T up to 10 K, were a prioriincapable of giving adequate results. We see that a shortened Cp(T)/T expansion of type (36) (due to taking c7 = 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 ΘD(0)/50). 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 ΘD(0) values (closely approaching to those quoted in Table 3), provided the fitting ranges had been strictly limited to T<6 K, for ZnSe, and T<4 K, for ZnTe.In summary one can assess that the obviously erroneousΘD(0) 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 ΘD(0) values. These moderate uncertainties are obviously due to certain contributions of degenerate electronic systems to the measured heat capacities, Cp(T), which found their numerical expression in nonvanishing magnitudes for the coefficients c1 (cf. Tables 1 and 3). Concerning the actual magnitudes of the latter we see, among other things, that the respective c1 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 ρp(T)≡Cp(T)/T3 data points [46] in the insets of Figures 2 and 3 and the respective local maximum behaviours of the corresponding ΘD(T) data points in the insets of Figures 4, 5, and 6). Our analyses have shown that the presently available cryogenic Cp(T) 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ΘD(0)≅ 270 K (±2 K), including the particular value of ΘD(0) = 271.4 K (Table 3), for ZnSe, is in good accordance with earlier caloric values of ΘD(0)≅ 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 ΘD(0)≅ 220 K (±2 K), including the particular value of ΘD(0) = 219.3 K (Table 3), for ZnTe, is in good agreement with several ΘDel(0) values derived from elastic constants. This concerns in particular the ΘDel(0) 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 ΘDel(0) 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 ΘDel(0) 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ΘD(T) and ΘDh(T) 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 Cp(T) dependence into the low-capacity formula (see (25)) for effective Debye temperatures, one comes readily to a relatively simple (algebraic) ΘD(T) expression of type [13](37)ΘDl0T⟶4π4/5·CVh∞c1/T2+c3+c5T2+c7T43=ΘD01+c1/T2+c5T2+c7T4/c33, 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 ΘDl0(T) dependencies (solid curves in the insets of Figure 6). The associated ΘDhl0(T) 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 (c1→0). With respect to this characteristic ΘDhl0(T→0) behaviour it is necessary, of course, to point out that simple algebraic formulas [13] of type (37), for ΘDl0(T) and ΘDhl0(T), apply only to very narrow cryogenic intervals, 0 < T < ΘD(0)/30, where the respective Debye temperature curves show an exclusivelyconvex shape (cf. the insets of Figure 6). Equivalently one can state that algebraic formula of the type (37) isnot applicable to temperatures higher than ΘD(0)/30 (i.e., to T > 9 K for ZnSe or T > 7 K ZnTe), where the (nearly coinciding) ΘDl0(T) and ΘDhl0(T) dependencies are showingconcave 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,ΘDh(T), from their limiting (T→ 0) levels, ΘD(0), to the respective minima, ΘDmin≡ΘD(Tmin), is described in good approximation by the relatively simple (algebraic) low-temperature formula, ΘDhl0(T) (see (30)). We see from the graphical representation of the latter by dotted curves (in the insets of Figures 4 and 5) that this ΘDhl0(T) expression provides, indeed, close approaches to the exact ΘDh(T) curves from 0 to Tmin≈ΘD(0)/15. With respect to these limited cryogenic regions, 0<T<ΘD(0)/15, one can thus describe the monotonic decrease of the ratios ΘD(0)/ΘDh(T) in good approximation by the relationship [13, 14](38)ΘD0ΘDhT≅ΘD0ΘDhl0T=ρhTρh03, in accordance with (30). The latter allows us to establish a useful approximate relationship between the minima, ΘDmin=ΘDh(Tmin), of the ΘDh(T)→ΘD(T) curves, on the one hand, and the maxima, ρmax=ρh(Tmax) of the ρh(T)→ρp(T) curves, on the other hand. Observing that the relative shifts of Tmin(Θ) positions (indicated in the insets of Figures 4 and 5) versus the Tmax(ρ) 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 ΘDh(T) curves by ΘDmin≅ΘDh(Tmax(ρ)). Accordingly, by specializing (38) for T→Tmax(ρ)≅Tmin(Θ), we come to an approximate relationship between the extreme values of the ΘDh(T) versus ρh(T) curves of the form [14](39)ΘD0ΘDmin≅ρmaxρh03=ρmaxc33.Table 4
Pieces of evidence of the close correlation between the characteristic cryogenic non-Debye features of theρp/h(T)≡Cp/Vh(T)/T3 curves (see the insets of Figures 2 and 3) and the respective Debye temperature curves, ΘD(T) (see Figures 4 and 5). The prominent features of the ρp/h(T) curves are characterized by the T→0 levels for the lattice heat capacity, ρh(0)=c3, the coordinates of their maxima, Tmax and ρmax≡ρ(Tmax), and the respective elevation ratios, ρmax/ρh(0) > 1. The ΘD(h)(T) curves are characterized above all by the limiting (T→0) levels, ΘD(0), the coordinates of their minima, Tmin and ΘDmin≡ΘD(Tmin), and the corresponding drop ratios, ΘD(min)/ΘD(0) < 1. The material-specific ratios ΘD(min)/Tmin are quoted for comparisons with the general condition ΘD(min)/Tmin>14 (cf. Section 5) of close approach to asymptotic behaviour due to (25). Furthermore, the maximum values of the effective Debye temperatures, ΘD(max), the theoretical high-temperature limiting value of the true Debye temperature, ΘDh(∞) (42), and several characteristic ratios, the approximate equality of which, for both materials, is indicating their high degree of similarity, are quoted. Particularly the large deviations from unity of the cryogenic drop ratios, ΘD(min)/ΘD(0)≈ 0.76, and of the ratios for the subsequent rise to maxima, ΘD(max)/ΘD(min)≈ 1.5, illustrate the considerable strength of T-related variations of Debye temperatures for both materials under study.
Material
ZnSe
ZnTe
ρ
h
(
0
)
=
c
3 (μJ⋅K−4 mol−1)
197.5
365.0
T
m
a
x (K)
16.9
13.6
ρ
m
a
x (μJ⋅K−4 mol−1)
440.0
855.7
ρ
m
a
x
/
ρ
h
(
0
)
2.228
2.344
ρ
m
a
x
/
ρ
h
(
0
)
3
1.306
1.328
Θ
D
(
0
) (K)
270
220
T
m
i
n (K)
17.3
13.9
Θ
D
(
0
)
/
T
m
i
n
15.6
15.8
Θ
D
(
m
i
n
) (K)
206.3
165.4
Θ
D
(
0
)
/
Θ
D
m
i
n
1.309
1.330
Θ
D
(
m
i
n
)
/
T
m
i
n
11.9
11.9
Θ
D
(
m
a
x
) (K)
305.7
252.3
Θ
D
(
m
a
x
)
/
Θ
D
(
m
i
n
)
1.482
1.525
Θ
D
h
(
∞
) (K)
317
262
Θ
D
(
0
)
/
Θ
D
h
(
∞
)
0.852
0.840
Θ
D
(
m
i
n
)
/
Θ
D
h
(
∞
)
0.651
0.631Comparing the actualΘD(0)/ΘDmin and ρmax/ρh(0)3=ρmax/c33 values (quoted in Table 4) we see that approximate relation (39) is indeed well fulfilled (within deviations smaller than 0.25%) by the empirical ΘDmin and ρmax values for ZnSe and ZnTe.It is instructive to compare the presently obtained (empirical)ΘDmin 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) ΘDmin values of about 207 K (cf. Figure 2 in [34]) and 218 K (cf. Figure  14 in [36]) are in reasonable agreement with the presently determined (empirical) ΘDmin value of 206.3 K (cf. Figure 4 and Table 4). In contrast to the latter, however, the theoretical ΘD(T) curve presented in Figure 7 of [87] shows an exclusively monotonic behaviour for any T > 10 K, without giving any indication of a possible local minimum behaviour in the vicinity of Tmin≈ 17 K (cf. the inset of Figure 4). It is suggested by Figure 7 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 ΘDmin value. A detailed numerical comparison with Figure 4 in the present paper (including Table 2) reveals, consequently, that the ΘD(T) curve shown for ZnSe in Figure  7 of [87] is quiteinadequate 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 ΘD(T>16K) curve shown for ZnS in Figure 7 of the same paper).For the case of ZnTe, we see that the previously estimatedΘDmin 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) ΘDmin value of 165.4 K (cf. Figure 5 and Table 4).Concerning the physical interpretation of the occurrence of the closely correlated maxima and minima of the respectiveρp/h(T) and ΘD/Dh(T) curves, we reconsider briefly the typical experimental manifestation of a heat capacity contribution, CE(T)=(ΘE/2T)2/(sinh(ΘE/2T))2, which is due to an individual Einstein oscillator located at εE=kBΘE. It can easily be shown that the maximum of the respective ρE(T)≡CE(T)/T3 curve is located at Tmax(E)=0.203ΘE (according to Figure A1 of the supporting information to [50]). Thus, when we connect the maximum of a given ρp/h(T) curve in roughest approximation with the contribution of a single Einstein oscillator, we can estimate its characteristic phonon temperature to amount to ΘE=4.93·Tmax(E). Accordingly, one can readily derive from the known Tmax positions (of 16.9 K for ZnSe and 13.6 K for ZnTe) values of ΘE = 83.3 K, for ZnSe, and ΘE = 70.0 K, for ZnTe, for the characteristic phonon temperatures of apparently dominating low-energy Einstein oscillators. Comparing these ΘE values with the set of Θn values obtained by means of the detailed analyses via the four-oscillator hybrid model (in Section 4) we see that they are ranging, for both materials, within the lower parts of the intervals between the first and second Einstein oscillators,(40)Θ1TA<ΘE<Θ2TA (cf. Table 1). Moreover, from the TA phonon dispersion relations shown for ZnSe and ZnTe in Figures  1(a) and 2(a) of [46] we infer that the characteristic TA phonon temperatures associated with the L and X points of the Brillouin zone boundaries, ΘTA(L) and ΘTA(X), are about 75 K and 93 K, for ZnSe, and about 59 K and 77 K, for ZnSe. Accordingly, we find that the effective Einstein oscillator temperatures, ΘE, are located close to the centre of gravity of the intervals between the ΘTA(L) and ΘTA(X) values; that is, ΘE≈(ΘTA(L)+ΘTA(X))/2 (for both materials). These findings can be considered as a clear confirmation of various earlier assessments [14, 20, 84], according to which the occurrences of the maximum behaviours of ρp/h(T) curves (like those shown in the insets of Figures 2 and 3) as well as the closely related minimum behaviours of the associated ΘD(h)(T) curves (shown in the insets of Figures 4 and 5) are characteristic manifestations of strongly increasing contributions of short-wave TA phonons to cryogenic heat capacities.


### 6.3. Debye Temperature Behaviour at Higher Temperatures

In view of the analogy of the analytical connections (1) and (29), between the xD(T)=ΘD(T)/T ratios with Cp(T) values, on the one hand, and the xDh(T)=ΘDh(T)/T ratios with CVh(T) values, on the other hand, it is obvious that the physical causes of the general differences, ΘDh(T)>ΘD(T) (see Figures 4 and 5), between “true” and effective Debye temperatures, must be sought in the mechanisms which are responsible for the inherent differences, Cp(T)-CVh(T)>0 (Section 2), between measured (isobaric) and theoretical (isochoric/harmonic) lattice heat capacities. For nonmagnetic materials these differences are well known to be due to(i)
possible electronic contributions,C(E)(T)=c1T (in (5)), which are responsible for qualitatively different asymptotic behaviours, that is, Cp(T)→c1T versus CVh(T)→c3T3, in the T→ 0 limit,(ii)
the cumulative effects of lattice expansion and lattice anharmonicities, which are responsible for markedly increasing differences atT>Th (cf. Figures 2 and 3).Concerning the latter, it has been found within a series of earlier studies [14–18, 21, 71, 88] that a relatively simple way for detection of temperature dependencies of “true” Debye temperatures, ΘDh(T), is offered via fittings of properly confined (i.e., exclusively concave) sections of ΘD(T) versus T-2 curves on the basis of (C.1) (see Figure 10, for the presently considered cases of ZnSe and ZnTe). This conventional method is as a rule capable of providing good results for materials showing a sufficiently high degree of phonon dispersion like the present ones, ΔP>1/3 (see Table 1; note that this frequently encountered state of affairs is in analogy to that found, for example, for silicon and germanium [21] as well as many II–VI and III–V materials [14, 50]). Thislarge dispersion regime is expressively visualized by ΘD(T) curves showing pronouncednonmonotonic behaviours in the vicinities of their minima (like those shown in Figures 4 and 5). The respective large degree of variability of ΘD(T) values, particularly with respect to the region T>Tmin, is characterized by rise ratios, ΘDmax/ΘDmin, that are significantlyhigher than unity (like ΘDmax/ΘDmin≈ 1.5, for ZnSe and ZnTe; cf. Table 4). In contrast to this, for materials with relativelysmall phonon dispersion, 1/3>ΔP>1/4 (like diamond and BN [50], or MgO [17]), for which the characteristic (material-specific) ΘDmax/ΘDmin ratios are known to be closely approaching to unity, this conventional method [14–18, 21, 71, 88] usually does not give good results for high-temperature behaviours, ΘDh(T→∞), of “true” Debye temperatures (for an illustration of significant residual uncertainties involved by a tentative application of this conventional ΘDh(T) estimation method to a low-dispersion material like MgO see [17]).In view of the limited applicability of this conventionalΘDh(T)→ΘD(T) fitting method (cf. Appendix C) we have displayed here in due detail (in Sections 2–4) a more general, alternative method, which is based on a properly adapted multioscillator hybrid model. This model involves, above all, a physically plausible decomposition of isobaric heat capacities, Cp(T) (see (5)), into prevailing harmonic lattice heat capacity parts, CVh(T)∝κh(T) (see (4)), and a rapidly increasing rest term, Cp(T)-CVh(T)∝(κh(T))2T, which represents the cumulative effects of lattice expansion and anharmonicities. Using this analytical framework for careful least-mean-square fittings of comprehensive sets of measured Cp(T) data, we have determined the respective parameter sets (Table 1). On the basis of the latter it was then possible to calculate in unambiguous way the correspondingharmonic lattice heat capacities, CVh(T)=CVh(∞)κh(T) (using (9), in combination with (14) to (17)). Consequently, we could here readily determine the T-dependencies of the “true” Debye temperatures, ΘDh(T), using the high-precision formulas (27) and (28). The respective ΘDh(T) dependencies are shown (by dashed curves) in Figures 4 and 5.Comparing the limitingΘDh(T→∞) values listed in Table 4 (indicated by arrows in Figures 4 and 5) with their counterparts resulting from the conventional estimation procedure listed in Table 10 (indicated in Figure 10 by the intercepts of the dashed curves with the ordinate), we can discern only minor differences (of less than 0.3%) between the limiting ΘDh(∞) values estimated via the two alternative analysis procedures. The slightly different results due to the alternative analyses procedures can thus be summarized in the form(41)ΘDh∞≅317K±1Kfor  ZnSe,ΘDh∞≅263K±1Kfor  ZnTe.The rather small size of residual uncertainties can be considered as clear indications of an obviously high degree of internal consistency and compatibility of the two qualitatively differentΘDh(∞) estimation procedures.An additional confirmation of the apparent adequacy of these caloricΘDh(∞) values (see (41)) can be obtained by comparison with independent theoretical estimations based on first-principle calculations of PDOS spectra. It is well known that the magnitudes of ΘDh(∞) are connected in unambiguous way with the second moments, μP(2), by [14, 16, 17, 21, 54, 56, 71, 72](42)ΘDh∞=5/3μP22kB≡5/32·εP2kB,where εP(2)≡μP(2)2 represents the associated (second order) phonon energy. We have performed within a former study [56] evaluations of moments, μP(m) (see (12)), on the basis of theoretically calculated PDOS spectral functions, gP(ε), which had been published for a variety of group-IV, III–V, and II–VI materials in a variety of earlier papers. Among the numerical results listed in Table I of [56] one can find, among other things, corresponding ΘDh(∞) values (due to (42)) for both materials here under study, namely, 324 K or 319 K for ZnSe (detected from the spectra presented in Figure 8 of [36] and Figure  2 of [66], resp.), and 264 K or 259 K for ZnTe (detected from the spectra presented in Figure  9 of [36] and Figure 3 of [32], resp.). Instructive are also comparisons with graphical representations of precalculated ΘDh(T) curves, from which one can detect limiting ΘDh(T→∞) values of about 324 K for ZnSe and 259 K for ZnTe (as indicated by Figures  14 and 15 in [36]) and 259 K for ZnTe (as indicated by Figure  5 in [32]). In this way we satisfy ourselves that the presently obtained (caloric) ΘDh(∞) values (41) show only moderate deviations (<2.5%) from former theoretical estimates [32, 36].Concerning the high-temperature dependencies of theeffective (anharmonicity-related) Debye temperatures, ΘD(T), we would still like to point out that, owing to the strict limitation of their definition (see (1)) to ratios Cp(T)/CVh(∞)≤1, their magnitudes are automatically dropping to zero at those positions, Tf, where the individual Cp(T) curves are crossing the classical Delong-Petit level; that is, (43)ΘDTf=0forCpTf=CVh∞=6R.This basic property is manifested by the rapidΘD→0 drops of the solid ΘD(T) curves (in Figures 4 and 5), which are located just at those positions, Tf = 401 K (for ZnSe) or 321 K (for ZnTe), where the respective Cp(T) curves are crossing the classical Dulong-Petit level, CVh(∞)=6R (as shown in Figures 2 and 3).


### 6.4. Characteristic PDOS-Related Phonon Energy Values

Within the frame of the presently considered four-oscillator hybrid model, the totality of quantitative information on the material-specific PDOS spectra is given (in Table1) in terms of the adjusted sets of effective Einstein oscillator temperatures, Θn (n = 1 to 4), and the respective weights, wn (n = 1 to 4), including the weights wC1 and wC2 for the associated continuous low-energy curve section components (of Debye and non-Debye type, resp.). The corresponding Einstein peakenergy positions, εn≡kBΘn (n = 1 to 4), are quoted in Table 5. On the basis of these sets of fitted parameter values one can readily calculate a variety of characteristic dispersion-related quantities, some of which are directly comparable with their theoretical counterparts being detectible from theoretically calculated PDOS spectra presented in various former papers [32, 36, 45, 46, 66, 89].Table 5
Phonon energy values for the 4 discrete (Einstein) oscillators,εn=kBΘn (n = 1 to 4, corresponding to the Θn values given in Table 1), the average phonon energies for the TA phonon sections, ε-TA (44), and the whole acoustic phonon energy sections, ε-TA&LA (45), including a series of equivalent (moment-related) phonon energies, εP(m) (due to (46) and (47)).
Material
ZnSe
ZnTe
ε
1
=
ε
T
A
1 (meV)
6.75
5.34
ε
2
=
ε
T
A
2 (meV)
9.43
7.58
ε
3
=
ε
-
L
A (meV)
18.36
14.40
ε
4
=
ε
-
T
O
&
L
O
≡
ε
-
U (meV)
27.10
22.69
ε
-
T
A (meV)
8.48
6.80
ε
-
T
A
&
L
A
=
ε
-
L (meV)
11.77
9.33
ε
P
(
-
5
/
2
) (meV)
10.13
8.18
ε
P
(
-
2
) (meV)
11.78
9.49
ε
P
(
-
3
/
2
) (meV)
13.12
10.58
ε
P
(
-
1
) (meV)
14.44
11.68
ε
P
(
-
1
/
2
) (meV)
15.79
12.83
ε
P
(
0
) (meV)
17.11
13.97
ε
P
(
1
)
=
ε
-
L
&
U (meV)
19.43
16.01
ε
P
(
2
) (meV)
21.17
17.55
ε
P
(
4
) (meV)
23.22
19.36
ε
P
(
6
) (meV)
24.27
20.29
ε
P
(
8
) (meV)
24.89
20.83
ε
P
(
10
) (meV)
25.30
21.18Table 6
The 8 lowest-order expansion coefficients,c2≤2m≤16DR, implied by the unprecedented high-capacity (high-temperature) xD2m-series expansion, κDhmR(xD) (A.11), for the conventional Debye function, κD(xD) (A.1). The same expansion coefficients, c2≤2m≤16DR, apply also the analogous x12m-series expansion, κC1hmR(x1) (A.15), for the Debye-like component, κC1(x1) (A.13), of the low-energy TA phonon section of the lattice heat capacity. With respect to the latter we give here also the set of the 8 lowest-order expansion coefficients, c2≤2m≤16NR, which are implied by the analogous x12m-series expansion, κC2hmR(x1) (A.23), for the associated non-Debye heat capacity component, κC2(x1) (A.17).
m
c
2
m
D
R
c
2
m
N
R
1
1
20
=
0.05
5
84
=
0.05952380952
2
1
1400
=
7.14285714286
×
1
0
-
4
13
10584
=
1.228269085
×
1
0
-
3
3
1
648000
=
1.54320987654
×
1
0
-
6
205
19559232
=
1.048098412
×
1
0
-
5
4
-
73
3492720000
=
-
2.09006161387
×
1
0
-
8
4559
160190110080
=
2.845993424
×
1
0
-
8
5
1
82555200000
=
1.21131073512
×
1
0
-
11
-
503
4485323082240
=
-
1.121435381
×
1
0
-
10
6
5591
5148967824000000
=
1.08584869650
×
1
0
-
12
-
478871
1585247736956083200
=
-
3.020795986
×
1
0
-
13
7
-
1879
212199886080000000
=
-
8.85485866515
×
1
0
-
15
-
48299
202404310545527029760
=
-
2.386261987
×
1
0
-
17
8
563743
11383109222284800000000
=
4.9524518213
×
1
0
-
17
887853331
16576922903367866373734400
=
5.355959827
×
1
0
-
17Table 7
High-accuracy sample values for Debye’s normalized heat capacity shape function,κD(xD) (A.1), which are directly comparable with their somewhat more restricted precedents listed in [6, 47] (up to 6 or 8 significant figures, resp.). The individual high-accuracy κD(xD) values have been calculated here for the interval 0≤xD≤1.8, via the unprecedented high-capacity (high-temperature) Taylor series expansion, κDhmR=7(xD) (A.11) and, for the interval 1.8≤xD≤40, by means of Debye’s conventional low-capacity (low-temperature) expansion, κDlnD=15(xD) (A.4) (note that a simple asymptotic expression, κD(xD)→77.92727286/xD3, applies to xD > 40). For comparisons with numerical results quoted in mathematical standard literature (see, e.g., Table 27.1 in [90]) we have listed here also the respective magnitudes of the conventional Debye function, D3(xD) (for 3-dimensional systems), D3(xD)=κD(xD)+3xD/(exp(xD)-1)/4 (according to (A.2) and (A.3)).
x
D
κ
D
(
x
D
)
D
3
(
x
D
)
x
D
κ
D
(
x
D
)
D
3
(
x
D
)
0.0
1.0000000000
1.0000000000
5.0
0.3686348236
0.1175974118
0.1
0.9995001785
0.9629999405
5.5
0.3132547732
0.09524080272
0.2
0.9980028536
0.9269990483
6.0
0.2655968891
0.07758132473
0.3
0.9955144242
0.8919951865
6.5
0.2250561623
0.06360434235
0.4
0.9920454895
0.8579848069
7.0
0.1908564046
0.05250585097
0.5
0.9876107521
0.8249629690
7.5
0.1621687997
0.04365502115
0.6
0.9822288834
0.7929233677
8.0
0.1381873779
0.03656029567
0.7
0.9759223556
0.7618583673
8.5
0.1181717647
0.03084031600
0.8
0.9687172408
0.7317590428
9.0
0.1014670944
0.02619989259
0.9
0.9606429816
0.7026152269
9.5
0.08750918665
0.02241065587
1.0
0.9517321357
0.6744155641
10.0
0.07582100303
0.01929576569
1.1
0.9420200996
0.6471475686
11.0
0.05773082402
0.01457049734
1.2
0.9315448136
0.6207976880
12.0
0.04478011684
0.01125032746
1.3
0.9203464535
0.5953513709
13.0
0.03534728217
0.008858858804
1.4
0.9084671129
0.5707931370
14.0
0.02835178383
0.007096677015
1.5
0.8959504801
0.5471066514
15.0
0.02307128730
0.005771263228
1.6
0.8828415127
0.5242747994
16.0
0.01901817415
0.004755893959
1.7
0.8691861161
0.5022797638
17.0
0.01585873694
0.003965212076
1.8
0.8550308264
0.4811031024
18.0
0.01336097602
0.003340449610
1.9
0.8404225040
0.4607258256
19.0
0.01136091855
0.002840309476
2.0
0.8254080384
0.4411284737
20.0
0.009740756601
0.002435220067
2.2
0.7943467192
0.4041938106
21.0
0.008414505229
0.002103638250
2.4
0.7622123728
0.3701368830
22.0
0.007318466926
0.001829621334
2.6
0.7293549760
0.3387926100
23.0
0.006404797367
0.001601201112
2.8
0.6961028658
0.3099951582
24.0
0.005637096943
0.001409274915
3.0
0.6627580442
0.2835798281
25.0
0.004987344231
0.001246836318
3.2
0.6295928776
0.2593847028
26.0
0.004433731487
0.001108432971
3.4
0.5968481116
0.2372520497
28.0
0.003549893920
0.0008874734947
3.6
0.5647320637
0.2170294741
30.0
0.002886195280
0.0007215488222
3.8
0.5334208309
0.1985708314
32.0
0.002378151635
0.0005945379091
4.0
0.5030593309
0.1817369138
34.0
0.001982680461
0.0004956701152
4.2
0.4737629932
0.1663959309
36.0
0.001670251904
0.0004175629760
4.4
0.4456199253
0.1524238050
38.0
0.001420164616
0.0003550411541
4.6
0.4186933946
0.1397043074
40.0
0.001217613638
0.0003044034095
4.8
0.3930244869
0.1281290579Table 8
Fitted values of the model-specific coefficients involved by the refined low-to-intermediate capacityxD(κ) formula, xDl9(0<κ≤κ9) (B.6), and for the couple of complementary high-accuracy low- and high-capacity Debye temperature formulas, xDl12(0<κ≤κc) (B.8) and xDh11(κc≤κ≤1) (B.9), which have been determined by least-mean-square fittings of a sufficiently comprehensive set of couples of high-accuracy κ versus xD curve points (in accordance with Table 7).
x
D
l
9
(
κ
) (see (B.6))
x
D
l
12
(
κ
) (see (B.8))
x
D
h
11
(
κ
) (see (B.9))
0
<
κ
≤
κ
9
=
0.8844
0
<
κ
≤
κ
c
=
0.260
0.260
=
κ
c
≤
κ
≤
1
p
=
0.071797185
a
2
=
2.8743354
×
10
-
1
b
2
=
6.0777673
×
10
-
3
g
1
=
3.042348
a
3
=
1.1377662
×
10
-
1
b
3
=
5.1466371
×
10
-
3
g
2
=
1.816742
a
4
=
-
9.5999179
×
10
-
2
b
4
=
-
6.8442924
×
10
-
3
g
3
=
-
6.892814
a
5
=
4.3062699
×
10
-
2
b
5
=
2.6051894
×
10
-
2
g
4
=
6.954053
a
6
=
-
4.1060421
×
10
-
3
b
6
=
-
3.6011454
×
10
-
2
g
5
=
5.456090
a
7
=
-
1.1955013
×
10
-
3
b
7
=
3.6169365
×
10
-
2
g
6
=
-
11.20018
a
8
=
4.9953204
×
10
-
4
b
8
=
-
3.7204523
×
10
-
2
g
7
=
2.022044
a
9
=
-
6.2681506
×
10
-
5
b
9
=
1.6804524
×
10
-
2
g
8
=
2.824376
a
10
=
3.3449835
×
10
-
6
b
10
=
-
3.0967328
×
10
-
3
—
p
=
1.2969
b
11
=
5.2219439
×
10
-
2
—
g
2
=
-
3.1541530
b
12
=
-
6.2935315
×
10
-
2
—
g
3
=
8.7184026
—Table 9
High-precision sample values of the characteristic Debye variable,xD(≡ΘD(T)/T), which have been calculated for a series of discrete (equidistant) κ-values via the couple of complementary low- and high-capacity formulas, xDl12(0<κ≤κc) (B.8) and xDh11(κc≤κ≤1) (B.9), using the model-specific coefficients listed in Table 8 (note that we have obtained precisely the same xD(κ) values also in alternative way via numerical integration procedures based directly on (A.1), for integration steps chosen to be of order 10-7).
κ
x
D
(
κ
)
κ
x
D
(
κ
)
κ
x
D
(
κ
)
κ
x
D
(
κ
)
κ
x
D
(
κ
)
0.01
19.825638
0.21
6.709531
0.41
4.666635
0.61
3.319275
0.81
2.100219
0.02
15.733335
0.22
6.568697
0.42
4.590077
0.62
3.258266
0.82
2.035428
0.03
13.737134
0.23
6.434348
0.43
4.514891
0.63
3.197532
0.83
1.969683
0.04
12.468384
0.24
6.305852
0.44
4.440994
0.64
3.137032
0.84
1.902849
0.05
11.557161
0.25
6.182662
0.45
4.368306
0.65
3.076725
0.85
1.834773
0.06
10.854261
0.26
6.064298
0.46
4.296752
0.66
3.016570
0.86
1.765274
0.07
10.285969
0.27
5.950343
0.47
4.226263
0.67
2.956525
0.87
1.694144
0.08
9.810966
0.28
5.840425
0.48
4.156772
0.68
2.896548
0.88
1.621130
0.09
9.403982
0.29
5.734216
0.49
4.088217
0.69
2.836595
0.89
1.545932
0.10
9.048539
0.30
5.631424
0.50
4.020538
0.70
2.776622
0.90
1.468177
0.11
8.733353
0.31
5.531788
0.51
3.953679
0.71
2.716581
0.91
1.387399
0.12
8.450382
0.32
5.435072
0.52
3.887587
0.72
2.656425
0.92
1.302998
0.13
8.193715
0.33
5.341064
0.53
3.822209
0.73
2.596103
0.93
1.214179
0.14
7.958888
0.34
5.249572
0.54
3.757496
0.74
2.535560
0.94
1.119852
0.15
7.742451
0.35
5.160421
0.55
3.693400
0.75
2.474740
0.95
1.018444
0.16
7.541692
0.36
5.073452
0.56
3.629874
0.76
2.413582
0.96
0.907538
0.17
7.354438
0.37
4.988520
0.57
3.566875
0.77
2.352021
0.97
0.783057
0.18
7.178928
0.38
4.905490
0.58
3.504357
0.78
2.289985
0.98
0.637031
0.19
7.013711
0.39
4.824241
0.59
3.442279
0.79
2.227398
0.99
0.448821
0.20
6.857583
0.40
4.744657
0.60
3.380599
0.80
2.164173
1.00
0.000000Table 10
Limiting values of the “true” (harmonic lattice) Debye temperatures,ΘDh(∞), including adjusted magnitudes of several even-order expansion coefficients, a2s, that are representing the weights of the individual (ΘDh(∞)/T)2s power terms in (C.1). These sets of parameter values have been determined via fittings (shown in Figure 10) of an appropriate selection of usable Debye temperature points, ΘDh(T)≅ΘD(T), implied by the experimental Cp(T) data points given by Kremer et al. [46] for natural ZnSe and ZnTe. Listed are also the magnitudes of the respective (Debye-model-related) cut-off phonon energies, εDh(∞)≡kBΘDh(∞), and several even-order phonon energy values, εP(2s)≡(μP(2s))1/2s (for 2s=2 to 10), which are resulting from the corresponding PDOS spectra moments, μP(2s), calculated on the basis of the analytical expressions (C.2a) to (C.2e) (note that the latter equations are exactly corresponding to the series of μP(2s) formulas originally derived within the frame of the supporting information to [21]; see also appendix B of [14] and/or the appendix B of [15]).
ZnSe
ZnTe
Intervals (K)
17 to 96
15 to 80
Θ
D
h
(
∞
) (K)
317.9
262.3
a
2
0.9259 × 10−2
0.9995 × 10−2
a
4
0.6895 × 10−4
0.7851 × 10−4
a
6
0.2869 × 10−6
0.3410 × 10−6
a
8
0.6533 × 10−9
0.7714 × 10−9
a
10
(0.7135 × 10−12)
(0.7038 × 10−12)
a
12
(0.2510 × 10−15)
—
ε
D
h
(
∞
) (meV)
27.40
22.60
ε
P
(
2
) (meV)
21.22
17.51
ε
P
(
4
) (meV)
23.48
19.45
ε
P
(
6
) (meV)
24.83
20.60
ε
P
(
8
) (meV)
25.75
21.38
ε
P
(
10
) (meV)
26.42
21.95Of importance for a physical understanding of the prominent role of short-wave TA phonons, which are mainly responsible for the occurrence of the characteristic non-Debye (maximum) behaviours ofCp(T)/T3 curves in the liquid hydrogen region (cf. the insets of Figures 2 and 3), is also the knowledge of the centres of gravity of TA phonon sections, ε-TA. The magnitudes of the latter are given, in accordance with the respective TA phonon normalization condition (see (23)), by [14](44)ε-TA=33wC14+5wC26+w1ε1+w2ε2 (cf. Eq. (18) in [14]). Concerning the presently estimated ε-TA values of 8.48 meV for ZnSe and 6.80 meV for ZnTe (as quoted in Table 5), we can make the satisfying statement that they are in good agreement with their counterparts (of about 8.4 meV or 8.5 meV, for ZnSe, and 6.8 meV or 6.9 meV, for ZnTe), which we have detected from the theoretically calculated PDOS spectra presented for ZnSe and ZnTe in Figure  3 of [45] or Figure  4 of [46], respectively. Accidentally, these ε-TA positions are even nearly coinciding with the positions of the highest TA peaks shown within the same figures (which are readily seen to be located at about 8.5 meV or 8.6 meV, for ZnSe, and between 6.8 meV and 6.9 meV, for ZnTe).Similarly informative are also comparisons of the average values of LA phonon energies,ε-LA, which are represented within the present formalism (due to (22)) by the third discrete oscillator, ε-LA=ε3. The presently estimated ε-LA values of 18.36 meV for ZnSe and 14.40 meV for ZnTe (cf. Table 5) are in reasonable agreement with their counterparts (of about 18.5 meV or 19.2 meV, for ZnSe, and 14.5 meV or 15.0 meV, for ZnTe), which we have detected, again, from the theoretically calculated PDOS spectra presented in Figure  3 of [45] or Figure  4 of [46], respectively.The average values for all acoustical (i.e., TA and LA) phonon energies,ε-TA&LA, are another characteristic quantity for comparisons between the results of numerical analyses of heat capacity data with theoretically calculated PDOS spectra. The analytical expression for the latter, which thus comprises the wholelower halves of PDOS spectral sections, is given within the present formalism by [14](45)ε-L=ε-TA&LA=23wC14+5wC26+w1ε1+w2ε2+ε36=2ε-TA+ε-LA3.The presently estimatedε-L=ε-TA&LA values of 11.77 meV for ZnSe and 9.33 meV for ZnTe (cf. Table 5) are found again to be in reasonable agreement with their counterparts of about 11.8 meV or 12.1 meV, for ZnSe, and 9.4 meV or 9.6 meV, for ZnTe, which followed immediately via (45) from the couples of ε-TA and ε-LA values for the average TA and LA phonon energies (which had been detected from the theoretically calculated PDOS spectra presented in Figure  3 of [45] or Figure  4 of [46], resp.). Furthermore we observe that the presently estimated (caloric) ε-L=ε-TA&LA values are in good agreement with our former estimations [56], which were based on theoretically calculated PDOS spectral functions published in various earlier papers, namely, 11.7 meV due to [66, 89], or 12.5 meV due to [36], for ZnSe, and 9.5 meV or 9.4 meV due to [32, 36], for ZnTe.Let us still consider theupper part of the PDOS spectra in consideration, which are consisting of one LO and two TO branches. The corresponding average optical phonon energy is given within the present formalism (due to (21)) just by the fourth discrete oscillator; that is, ε-U=ε-TO&LO=ε4. The presently estimated ε-U values (cf. Table 5) of 27.10 meV, for ZnSe, and 22.69 meV, for ZnTe, are found to be in more or less good agreement with their counterparts (of about 27.3 meV or 25.6 meV, for ZnSe, and 23.1 meV or 21.7 meV, for ZnTe), which we have detected again from the theoretically calculated PDOS spectra presented in Figure  4 of [46] and Figure  3 of [45], respectively. Informative is also a comparison of present (caloric) results with our former estimations [56] of ε-U values on the basis of PDOS spectral functions published in various earlier papers, namely, 27.3 meV, 27.1 meV, or 28.4 meV due to [36, 66, 89], for ZnSe, and 22.0 meV or 22.7 meV due to [32, 36], for ZnTe. Thus we are concerned in the case of the upper (optical) PDOS spectral sections with differences between approximate ε-U values of different origin up to a maximum order of 5% (for both ZnSe and ZnTe). At the same time we can see that the presently estimated (caloric) ε-U values (Table 5) are ranging near the middle between the highest and lowest ε-U estimates referring to the variety of PDOS spectra in consideration for the two materials under study. This shows that one has good chances to get reliable information on various prominent features of material-specific PDOS spectra even via dispersion-related analyses of heat capacity data.Very informative in this respect are, above all, comparisons between momentsμP(m), of material-specific PDOS spectral functions [3, 4, 16, 18, 19, 21, 50, 52, 56, 69–72], which are generally defined by (12) (for -3<m≠0). Within the present formalism of the four-oscillator hybrid PDOS model function, gP(ε) (i.e., (7), for nE = 4), one obtained for the corresponding moments, μP(m) (see (12)), the model-specific expression (13) [14, 21]. Yet, for visualizations of the m-dependencies of such PDOS spectra moments it is convenient to represent the individual mth-order moments, μP(m), in terms of their equivalent mth-order phonon energies [14, 16, 21, 50, 56], εP(m)≡μP(m)m. The latter are thus given within the frame of the present four-oscillator hybrid model (in accordance with (13)) by an expression of the form [14, 21, 50](46)εPm=ε13wC13+m+5wC25+m+w1+∑n=24wnεnε1m1/m,for-3<m<0,0<m.By means of am→0 limiting transition one obtains from (46) in particular(47)εPm⟶0=ε1exp-wC13-wC25+∑n=24wnlnεnε1.Collections of corresponding material-specificεP(m) values for a series of integers and half-integers (within the interval -3<m≤10) are listed in Table 5 and are represented by open circles in Figure 7. In view of the fact that the whole m-dependence of the εP(m) curve (due to (46)) divides into a couple ofcontinuous curve sections (for -3<m<0 and 0<m), the sequence of discrete εP(m) values (○) are perfectly connected by the respective solid curve sections. Furthermore we have shown in Figure 7 several sets of discrete εP(m) values which we have obtained via detailed numerical analyses of the theoretically calculated PDOS spectra presented for ZnSe and/or ZnTe in [32] (×), [36] (∗), [45] (◊), [46] (△), [66] (+), and [89] (☆). The comparisons between the different εP(m) values, in Figure 7, show a multitude of good agreements between the discrete εP(m) values due to the PDOS spectra published by different authors, including present (caloric) results (○). This concerns especially the region of negative orders (m<0), where the actual εP(m) values are predominantly due to the contributions of TA phonons. At the same time we are concerned with some markedly larger deviations in the region of positive orders (m > 0), where the monotonically increasing εP(m) values tend to be the stronger ones governed by optical phonons the higher their orders are. The maximum difference between estimated εP(10) values for ZnTe (due to [46] (△) versus [45] (◊)) amounts to about 6%, whereas, for ZnSe, the maximum difference (due to [89] (☆) versus [45] (◊)) amounts even to about 13%. Relatively large differences like the latter ones are obviously mainly due to significant qualitative differences between the upper (optical) sections of the theoretically calculated PDOS spectra presented by the corresponding authors (of [45] (◊) versus [89] (☆)). At the same time we can make the satisfactory statement that the presently obtained (caloric) εP(10) values (Table 5) differ from those we have derived from several other PDOS spectra in consideration for ZnSe ([36] (∗), [46] (△), and [66] (+)) and for ZnTe ([32] (×), [36] (∗), and [46] (△)) by less than 2%.Figure 7
Visualization of the effective phonon energy curves,εP(m>-3) (see (46); solid curves), for ZnSe and ZnTe. Included are the corresponding discrete εP(m) points ((○); cf. Table 5), for integer and several negative half integer orders. For comparisons with counterparts following from former theoretically calculated PDOS spectra we have shown here also the respective series of discrete εP(m) values obtained via detailed numerical analyses of the theoretically calculated PDOS spectra presented for ZnSe and/or ZnTe in [32] (×), [36] (∗), [45] (◊), [46] (△), [66] (+), and [89] (☆). The series of even-order εP(m) values (for m=2 to 10, listed in Table 10), which have been determined in Appendix C via the conventional analyses procedure for Debye temperature values representing the harmonic lattice, ΘDh(T)≅ΘD(T) (cf. Figure 10), are represented by squares (□).The latter findings underscore the above assessment, according to which one has as a rule good chances to get reliable information of various characteristic features of material-specific PDOS spectra by means of detailed dispersion-related analyses of heat capacity data. The rather fine quantification abilities of our elaborate analytical framework (displayed in Section4) concern obviously also such delicate quantities likehigher-order phonon energy values, εP(m≫1), and, consequently, even the corresponding higher-order moments, μP(m)=(εP(m))m, of PDOS spectra.Just with respect to the latter it might be instructive to consider still the alternative series of five even-orderεP(m) values (for m = 2 to 10), which have been determined in Appendix C by means of the conventional analyses procedure [14–18, 21, 71, 88] for Debye temperatures representing the harmonic lattice properties, ΘDh(T)≅ΘD(T) (the range of which is limited to T < 96 K, for ZnSe, and to T< 80 K, for ZnTe; cf. Figure 10 and Table 10). The corresponding alternative εP(m) values (listed in Table 10) are represented in Figure 7 by squares (□). Comparing the latter with their hybrid-model-related counterparts (○) (listed in the last five rows of Table 5) we find excellent agreements, at least, for the lowest two orders (m = 2 and 4). Towards higher orders (m = 6 to 10), however, we are concerned with growing deviations between the two alternative series of εP(m) values, which are seen to increase for ZnSe from 2.3% to 4.4% and for ZnTe from 1.5% to 3.6% (cf. Figure 7).Finally, let us still consider the behaviour (m-dependence) of the conventionalDebye-model-related energy curves, which are defined (in analogy to [16–19, 64, 71]) by an expression of the form(48)εDPm≡3+mμPm31/m=3+m31/mεPm. in accordance with the general definition, εP(m)≡μP(m)1/m, of the mth-order phonon energies [14, 21, 50, 56]. Observing that the m-dependencies of the two complementary εP(m) curve sections (for -3<m<0 and 0<m) are given by (46), we come thus readily (via (48)) to the material-specific εDP(m) dependencies, for ZnSe and ZnTe, we have visualized by corresponding dash-dot curves in Figure 7 (note that, for the special point m→0, the particular εDP(0) value follows from a corresponding limiting transition to adopt just the value εDP(m→0)=e1/3εP(0) (cf. [16, 17]), where εP(0) is given by (47)).A prominent feature of these (conventional) Debye-model-related energy curves,εDP(m) (see (48)), is due to the well-known property [16, 17, 19], according to which the m→-3 limit of this function is automatically (exactly) equal to the T→0 limit of the cut-off energy for Debye’s fictive phonon energy spectrum; that is,(49)εDPm⟶-3=wC1-1/3ε1=kBΘD0=εD0, in accordance with (20). The corresponding magnitudes of εD(0)=23.3 meV, for ZnSe, and εD(0) = 19.0 meV, for ZnTe, are indicated by solid triangles (▸) in Figure 7.Another characteristicεDP(m) value is that for m = 2, which follows readily from (48), in combination with (42), to be just equal to the T→∞ limit of the cut-off energy for Debye’s fictive (harmonic) phonon energy spectrum; that is,(50)εDPm=2=532εP2=kBΘDh∞≡εDh∞.The corresponding magnitudes ofεDh∞=27.3 meV, for ZnSe, and εDh(∞) = 22.6 meV, for ZnTe, are indicated by solid squares (■) in Figure 7. Accordingly, the ratios between these two limiting Debye energy (or temperature) values, εD(0)/εDh(∞)=ΘD(0)/ΘDh(∞), amount to about 0.85 or 0.84 for ZnSe and ZnTe, respectively (cf. Table 4). The significant differences (of about 15%) between these two limiting Debye temperature values are already, per se, clear indications for a pronouncednon-Debye character of both PDOS spectra in consideration.Actually, for the hypothetical case of a properly normalized phonon energy spectrum of Debye’s type, that is, for the fictive (exclusively quadratic) PDOS functiongD(ε)=3ε2/εDc3 [14, 21] (within the interval 0≤ε≤εDc, and gD(ε)=0, for εDc<ε<∞), the corresponding moments are readily seen to be given just by μD(m)=3εDcm/(3+m) [14]. Substituting hence the actual (material-specific) moments μP(m) in (48) by the moments μD(m) due to Debye’s fictive model, we see readily that expression (48) for εDP(m) reduces just to a constant(51)εDPm⟶εDDm≡3+mμDm31/m=εDc,for  any  order-3<m<∞.Such a constancy of anεDP(m) function (of type (48)), however, has never been observed in physical reality (not even in form of a relatively close approach of one or the other εDP(m) curve to a horizontal line). The rather strong variations of εDP(m) values shown in Figure 7, particularly the pronouncednonmonotonic (minimum) behaviours of both εDP(m) curves in the vicinity of m≈-2, can be looked upon as clear visual manifestations (and a global confirmation) of the obviously very pronouncednon-Debye character of the PDOS spectra [32, 36, 45, 46, 66, 89] for both materials under study.


## 6.1. Limiting Debye Temperatures and Electronic Heat Capacity Contributions

A provisional estimation ofΘD(0) values can be made in simple way by viewing the sets of cryogenic Debye temperature points, ΘD(T) (empty circles visualized in the insets of Figures 4, 5, and 6), which are exactly corresponding to the experimental Cp(T)/T3 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 < T < 10 K), we see readily that the respective ΘD(T) values are throughout rangingbelow the following upper boundaries:(31)ΘDT<270Kfor  ZnSe,ΘDT<220Kfor  ZnTe.Figure 6
Fittings via (35) of the Cp(T)/T data points (○) corresponding to the cryogenic Cp(T) data provided by Kremer et al. [46] for ZnSe and ZnTe [—, least-mean square fittings from 2 K to 8.7 K or to 6.4 K, resp.; ⋅⋅–⋅⋅–⋅⋅–⋅⋅, tail behaviour, for c7→0; ⋅–⋅–⋅–⋅, T→0 asymptote, for c7=0 and c5=0]. The insets show the effective ΘD(T) values (○) following directly (via (25)) from Kremer’s cryogenic Cp(T) data [46] for ZnSe and ZnTe. For comparison we have still included the ΘD(T) values (+) following from former Cp(T) data quoted for ZnSe by Birch [34], which show no indication for an appreciable electronic contribution. The continuous curves in the insets are resulting from (25), implying the cryogenic Cp(T) series expansion (see (33)), with fitted Taylor series expansion coefficients listed in Table 3. [—, ΘD(T) curve, for c1 > 0;   – – – – –,    ΘDh(T) curve, for c1→ 0.].Table 3
Conventional Taylor series expansion coefficients determined by least-mean square fittings of cryogenicCp(T)/T data [46] for ZnSe and ZnTe (see Figure 6) via (35). The respective ΘD(0) values are resulting from fitted c3 values via (34).
ZnSe
ZnTe
Intervals (K)
2.0 to 8.7
1.9 to 6.4
c
1 (mJ⋅K−2 mol−1)
0.1339
0.0409
c
3 (mJ⋅K−4 mol−1)
0.1937
0.3686
c
5 (mJ⋅K−6 mol−1)
0.602 × 10−3
0.537 × 10−3
c
7 (mJ⋅K−8 mol−1)
0.815 × 10−5
0.638 × 10−4
Θ
D
(
0
) (K)
271.7
219.3It 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 Section4) that the limiting (T→ 0) values of Debye temperatures, ΘD(0) (cf. Table 1), are nearly coinciding (within possible uncertainties of order ±2 K) with these upper boundaries (31); that is, (32)ΘD0≅270K±2K,for  ZnSe,ΘD0≅220K±2K,for  ZnTe.These results are, of course, in clear contrast to the formerΘD(0) 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ΘD(0)/30), 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](33)CpT=∑n=03c2n+1T2n+1,where the c1T term [4, 6, 13, 23, 60–63] 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, c3T3, which corresponds to Debye’s classical (T→0 limiting) cubic law [1, 3, 6, 8, 21], CVh(T→0)∝T3. The corresponding proportionality factor, c3, is known to be connected in a simple way with the limiting Debye temperature, ΘD(0), by the relation [4, 6, 13, 15, 46, 59, 62](34)ΘD0=4π4CVh∞5c33.Consequently, within the frame of the frequently adoptedCp(T)/T versus T2 representation [4, 6, 13, 15, 46, 62, 63, 83](35)CpTT=c1+c3T2+c5T4+c7T6,the coefficient c3 corresponds just to the slopes of the Cp(T)/T versus T2 curves in the T→0 limit (which are represented in Figure 6 by dash-dot lines). However, these asymptotic lines used to be nearly coincident with the actual Cp(T)/T versus T2 curves only within very narrow intervals, from 0 to about ΘD(0)/100 [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 Cp(T) data for the interval from 0 to 2 K, there remains only the possibility of determining c3 in combination with the other expansion coefficients (i.e., c1, c5, and c7) by careful fittings via (35) of the somewhat larger Cp(T)/T versus T2 curve sections extending from 2 K up to about ΘD(0)/30 (see the solid curves in Figure 6). The parameter values obtained in this way are listed in Table 3. Viewing the respective ΘD(0) 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 Figure6 one can now easily see why the rougher simulations of Cp(T)/T data (in [46]) via a shortened series expansion of type(36)CpTT≅c3T2+c5T4,for T up to 10 K, were a prioriincapable of giving adequate results. We see that a shortened Cp(T)/T expansion of type (36) (due to taking c7 = 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 ΘD(0)/50). 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 ΘD(0) values (closely approaching to those quoted in Table 3), provided the fitting ranges had been strictly limited to T<6 K, for ZnSe, and T<4 K, for ZnTe.In summary one can assess that the obviously erroneousΘD(0) 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 ΘD(0) values. These moderate uncertainties are obviously due to certain contributions of degenerate electronic systems to the measured heat capacities, Cp(T), which found their numerical expression in nonvanishing magnitudes for the coefficients c1 (cf. Tables 1 and 3). Concerning the actual magnitudes of the latter we see, among other things, that the respective c1 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 ρp(T)≡Cp(T)/T3 data points [46] in the insets of Figures 2 and 3 and the respective local maximum behaviours of the corresponding ΘD(T) data points in the insets of Figures 4, 5, and 6). Our analyses have shown that the presently available cryogenic Cp(T) 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ΘD(0)≅ 270 K (±2 K), including the particular value of ΘD(0) = 271.4 K (Table 3), for ZnSe, is in good accordance with earlier caloric values of ΘD(0)≅ 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 ΘD(0)≅ 220 K (±2 K), including the particular value of ΘD(0) = 219.3 K (Table 3), for ZnTe, is in good agreement with several ΘDel(0) values derived from elastic constants. This concerns in particular the ΘDel(0) 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 ΘDel(0) 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 ΘDel(0) 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ΘD(T) and ΘDh(T) 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 Cp(T) dependence into the low-capacity formula (see (25)) for effective Debye temperatures, one comes readily to a relatively simple (algebraic) ΘD(T) expression of type [13](37)ΘDl0T⟶4π4/5·CVh∞c1/T2+c3+c5T2+c7T43=ΘD01+c1/T2+c5T2+c7T4/c33, 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 ΘDl0(T) dependencies (solid curves in the insets of Figure 6). The associated ΘDhl0(T) 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 (c1→0). With respect to this characteristic ΘDhl0(T→0) behaviour it is necessary, of course, to point out that simple algebraic formulas [13] of type (37), for ΘDl0(T) and ΘDhl0(T), apply only to very narrow cryogenic intervals, 0 < T < ΘD(0)/30, where the respective Debye temperature curves show an exclusivelyconvex shape (cf. the insets of Figure 6). Equivalently one can state that algebraic formula of the type (37) isnot applicable to temperatures higher than ΘD(0)/30 (i.e., to T > 9 K for ZnSe or T > 7 K ZnTe), where the (nearly coinciding) ΘDl0(T) and ΘDhl0(T) dependencies are showingconcave 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,ΘDh(T), from their limiting (T→ 0) levels, ΘD(0), to the respective minima, ΘDmin≡ΘD(Tmin), is described in good approximation by the relatively simple (algebraic) low-temperature formula, ΘDhl0(T) (see (30)). We see from the graphical representation of the latter by dotted curves (in the insets of Figures 4 and 5) that this ΘDhl0(T) expression provides, indeed, close approaches to the exact ΘDh(T) curves from 0 to Tmin≈ΘD(0)/15. With respect to these limited cryogenic regions, 0<T<ΘD(0)/15, one can thus describe the monotonic decrease of the ratios ΘD(0)/ΘDh(T) in good approximation by the relationship [13, 14](38)ΘD0ΘDhT≅ΘD0ΘDhl0T=ρhTρh03, in accordance with (30). The latter allows us to establish a useful approximate relationship between the minima, ΘDmin=ΘDh(Tmin), of the ΘDh(T)→ΘD(T) curves, on the one hand, and the maxima, ρmax=ρh(Tmax) of the ρh(T)→ρp(T) curves, on the other hand. Observing that the relative shifts of Tmin(Θ) positions (indicated in the insets of Figures 4 and 5) versus the Tmax(ρ) 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 ΘDh(T) curves by ΘDmin≅ΘDh(Tmax(ρ)). Accordingly, by specializing (38) for T→Tmax(ρ)≅Tmin(Θ), we come to an approximate relationship between the extreme values of the ΘDh(T) versus ρh(T) curves of the form [14](39)ΘD0ΘDmin≅ρmaxρh03=ρmaxc33.Table 4
Pieces of evidence of the close correlation between the characteristic cryogenic non-Debye features of theρp/h(T)≡Cp/Vh(T)/T3 curves (see the insets of Figures 2 and 3) and the respective Debye temperature curves, ΘD(T) (see Figures 4 and 5). The prominent features of the ρp/h(T) curves are characterized by the T→0 levels for the lattice heat capacity, ρh(0)=c3, the coordinates of their maxima, Tmax and ρmax≡ρ(Tmax), and the respective elevation ratios, ρmax/ρh(0) > 1. The ΘD(h)(T) curves are characterized above all by the limiting (T→0) levels, ΘD(0), the coordinates of their minima, Tmin and ΘDmin≡ΘD(Tmin), and the corresponding drop ratios, ΘD(min)/ΘD(0) < 1. The material-specific ratios ΘD(min)/Tmin are quoted for comparisons with the general condition ΘD(min)/Tmin>14 (cf. Section 5) of close approach to asymptotic behaviour due to (25). Furthermore, the maximum values of the effective Debye temperatures, ΘD(max), the theoretical high-temperature limiting value of the true Debye temperature, ΘDh(∞) (42), and several characteristic ratios, the approximate equality of which, for both materials, is indicating their high degree of similarity, are quoted. Particularly the large deviations from unity of the cryogenic drop ratios, ΘD(min)/ΘD(0)≈ 0.76, and of the ratios for the subsequent rise to maxima, ΘD(max)/ΘD(min)≈ 1.5, illustrate the considerable strength of T-related variations of Debye temperatures for both materials under study.
Material
ZnSe
ZnTe
ρ
h
(
0
)
=
c
3 (μJ⋅K−4 mol−1)
197.5
365.0
T
m
a
x (K)
16.9
13.6
ρ
m
a
x (μJ⋅K−4 mol−1)
440.0
855.7
ρ
m
a
x
/
ρ
h
(
0
)
2.228
2.344
ρ
m
a
x
/
ρ
h
(
0
)
3
1.306
1.328
Θ
D
(
0
) (K)
270
220
T
m
i
n (K)
17.3
13.9
Θ
D
(
0
)
/
T
m
i
n
15.6
15.8
Θ
D
(
m
i
n
) (K)
206.3
165.4
Θ
D
(
0
)
/
Θ
D
m
i
n
1.309
1.330
Θ
D
(
m
i
n
)
/
T
m
i
n
11.9
11.9
Θ
D
(
m
a
x
) (K)
305.7
252.3
Θ
D
(
m
a
x
)
/
Θ
D
(
m
i
n
)
1.482
1.525
Θ
D
h
(
∞
) (K)
317
262
Θ
D
(
0
)
/
Θ
D
h
(
∞
)
0.852
0.840
Θ
D
(
m
i
n
)
/
Θ
D
h
(
∞
)
0.651
0.631Comparing the actualΘD(0)/ΘDmin and ρmax/ρh(0)3=ρmax/c33 values (quoted in Table 4) we see that approximate relation (39) is indeed well fulfilled (within deviations smaller than 0.25%) by the empirical ΘDmin and ρmax values for ZnSe and ZnTe.It is instructive to compare the presently obtained (empirical)ΘDmin 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) ΘDmin values of about 207 K (cf. Figure 2 in [34]) and 218 K (cf. Figure  14 in [36]) are in reasonable agreement with the presently determined (empirical) ΘDmin value of 206.3 K (cf. Figure 4 and Table 4). In contrast to the latter, however, the theoretical ΘD(T) curve presented in Figure 7 of [87] shows an exclusively monotonic behaviour for any T > 10 K, without giving any indication of a possible local minimum behaviour in the vicinity of Tmin≈ 17 K (cf. the inset of Figure 4). It is suggested by Figure 7 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 ΘDmin value. A detailed numerical comparison with Figure 4 in the present paper (including Table 2) reveals, consequently, that the ΘD(T) curve shown for ZnSe in Figure  7 of [87] is quiteinadequate 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 ΘD(T>16K) curve shown for ZnS in Figure 7 of the same paper).For the case of ZnTe, we see that the previously estimatedΘDmin 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) ΘDmin value of 165.4 K (cf. Figure 5 and Table 4).Concerning the physical interpretation of the occurrence of the closely correlated maxima and minima of the respectiveρp/h(T) and ΘD/Dh(T) curves, we reconsider briefly the typical experimental manifestation of a heat capacity contribution, CE(T)=(ΘE/2T)2/(sinh(ΘE/2T))2, which is due to an individual Einstein oscillator located at εE=kBΘE. It can easily be shown that the maximum of the respective ρE(T)≡CE(T)/T3 curve is located at Tmax(E)=0.203ΘE (according to Figure A1 of the supporting information to [50]). Thus, when we connect the maximum of a given ρp/h(T) curve in roughest approximation with the contribution of a single Einstein oscillator, we can estimate its characteristic phonon temperature to amount to ΘE=4.93·Tmax(E). Accordingly, one can readily derive from the known Tmax positions (of 16.9 K for ZnSe and 13.6 K for ZnTe) values of ΘE = 83.3 K, for ZnSe, and ΘE = 70.0 K, for ZnTe, for the characteristic phonon temperatures of apparently dominating low-energy Einstein oscillators. Comparing these ΘE values with the set of Θn values obtained by means of the detailed analyses via the four-oscillator hybrid model (in Section 4) we see that they are ranging, for both materials, within the lower parts of the intervals between the first and second Einstein oscillators,(40)Θ1TA<ΘE<Θ2TA (cf. Table 1). Moreover, from the TA phonon dispersion relations shown for ZnSe and ZnTe in Figures  1(a) and 2(a) of [46] we infer that the characteristic TA phonon temperatures associated with the L and X points of the Brillouin zone boundaries, ΘTA(L) and ΘTA(X), are about 75 K and 93 K, for ZnSe, and about 59 K and 77 K, for ZnSe. Accordingly, we find that the effective Einstein oscillator temperatures, ΘE, are located close to the centre of gravity of the intervals between the ΘTA(L) and ΘTA(X) values; that is, ΘE≈(ΘTA(L)+ΘTA(X))/2 (for both materials). These findings can be considered as a clear confirmation of various earlier assessments [14, 20, 84], according to which the occurrences of the maximum behaviours of ρp/h(T) curves (like those shown in the insets of Figures 2 and 3) as well as the closely related minimum behaviours of the associated ΘD(h)(T) curves (shown in the insets of Figures 4 and 5) are characteristic manifestations of strongly increasing contributions of short-wave TA phonons to cryogenic heat capacities.


## 6.3. Debye Temperature Behaviour at Higher Temperatures

In view of the analogy of the analytical connections (1) and (29), between the xD(T)=ΘD(T)/T ratios with Cp(T) values, on the one hand, and the xDh(T)=ΘDh(T)/T ratios with CVh(T) values, on the other hand, it is obvious that the physical causes of the general differences, ΘDh(T)>ΘD(T) (see Figures 4 and 5), between “true” and effective Debye temperatures, must be sought in the mechanisms which are responsible for the inherent differences, Cp(T)-CVh(T)>0 (Section 2), between measured (isobaric) and theoretical (isochoric/harmonic) lattice heat capacities. For nonmagnetic materials these differences are well known to be due to(i)
possible electronic contributions,C(E)(T)=c1T (in (5)), which are responsible for qualitatively different asymptotic behaviours, that is, Cp(T)→c1T versus CVh(T)→c3T3, in the T→ 0 limit,(ii)
the cumulative effects of lattice expansion and lattice anharmonicities, which are responsible for markedly increasing differences atT>Th (cf. Figures 2 and 3).Concerning the latter, it has been found within a series of earlier studies [14–18, 21, 71, 88] that a relatively simple way for detection of temperature dependencies of “true” Debye temperatures, ΘDh(T), is offered via fittings of properly confined (i.e., exclusively concave) sections of ΘD(T) versus T-2 curves on the basis of (C.1) (see Figure 10, for the presently considered cases of ZnSe and ZnTe). This conventional method is as a rule capable of providing good results for materials showing a sufficiently high degree of phonon dispersion like the present ones, ΔP>1/3 (see Table 1; note that this frequently encountered state of affairs is in analogy to that found, for example, for silicon and germanium [21] as well as many II–VI and III–V materials [14, 50]). Thislarge dispersion regime is expressively visualized by ΘD(T) curves showing pronouncednonmonotonic behaviours in the vicinities of their minima (like those shown in Figures 4 and 5). The respective large degree of variability of ΘD(T) values, particularly with respect to the region T>Tmin, is characterized by rise ratios, ΘDmax/ΘDmin, that are significantlyhigher than unity (like ΘDmax/ΘDmin≈ 1.5, for ZnSe and ZnTe; cf. Table 4). In contrast to this, for materials with relativelysmall phonon dispersion, 1/3>ΔP>1/4 (like diamond and BN [50], or MgO [17]), for which the characteristic (material-specific) ΘDmax/ΘDmin ratios are known to be closely approaching to unity, this conventional method [14–18, 21, 71, 88] usually does not give good results for high-temperature behaviours, ΘDh(T→∞), of “true” Debye temperatures (for an illustration of significant residual uncertainties involved by a tentative application of this conventional ΘDh(T) estimation method to a low-dispersion material like MgO see [17]).In view of the limited applicability of this conventionalΘDh(T)→ΘD(T) fitting method (cf. Appendix C) we have displayed here in due detail (in Sections 2–4) a more general, alternative method, which is based on a properly adapted multioscillator hybrid model. This model involves, above all, a physically plausible decomposition of isobaric heat capacities, Cp(T) (see (5)), into prevailing harmonic lattice heat capacity parts, CVh(T)∝κh(T) (see (4)), and a rapidly increasing rest term, Cp(T)-CVh(T)∝(κh(T))2T, which represents the cumulative effects of lattice expansion and anharmonicities. Using this analytical framework for careful least-mean-square fittings of comprehensive sets of measured Cp(T) data, we have determined the respective parameter sets (Table 1). On the basis of the latter it was then possible to calculate in unambiguous way the correspondingharmonic lattice heat capacities, CVh(T)=CVh(∞)κh(T) (using (9), in combination with (14) to (17)). Consequently, we could here readily determine the T-dependencies of the “true” Debye temperatures, ΘDh(T), using the high-precision formulas (27) and (28). The respective ΘDh(T) dependencies are shown (by dashed curves) in Figures 4 and 5.Comparing the limitingΘDh(T→∞) values listed in Table 4 (indicated by arrows in Figures 4 and 5) with their counterparts resulting from the conventional estimation procedure listed in Table 10 (indicated in Figure 10 by the intercepts of the dashed curves with the ordinate), we can discern only minor differences (of less than 0.3%) between the limiting ΘDh(∞) values estimated via the two alternative analysis procedures. The slightly different results due to the alternative analyses procedures can thus be summarized in the form(41)ΘDh∞≅317K±1Kfor  ZnSe,ΘDh∞≅263K±1Kfor  ZnTe.The rather small size of residual uncertainties can be considered as clear indications of an obviously high degree of internal consistency and compatibility of the two qualitatively differentΘDh(∞) estimation procedures.An additional confirmation of the apparent adequacy of these caloricΘDh(∞) values (see (41)) can be obtained by comparison with independent theoretical estimations based on first-principle calculations of PDOS spectra. It is well known that the magnitudes of ΘDh(∞) are connected in unambiguous way with the second moments, μP(2), by [14, 16, 17, 21, 54, 56, 71, 72](42)ΘDh∞=5/3μP22kB≡5/32·εP2kB,where εP(2)≡μP(2)2 represents the associated (second order) phonon energy. We have performed within a former study [56] evaluations of moments, μP(m) (see (12)), on the basis of theoretically calculated PDOS spectral functions, gP(ε), which had been published for a variety of group-IV, III–V, and II–VI materials in a variety of earlier papers. Among the numerical results listed in Table I of [56] one can find, among other things, corresponding ΘDh(∞) values (due to (42)) for both materials here under study, namely, 324 K or 319 K for ZnSe (detected from the spectra presented in Figure 8 of [36] and Figure  2 of [66], resp.), and 264 K or 259 K for ZnTe (detected from the spectra presented in Figure  9 of [36] and Figure 3 of [32], resp.). Instructive are also comparisons with graphical representations of precalculated ΘDh(T) curves, from which one can detect limiting ΘDh(T→∞) values of about 324 K for ZnSe and 259 K for ZnTe (as indicated by Figures  14 and 15 in [36]) and 259 K for ZnTe (as indicated by Figure  5 in [32]). In this way we satisfy ourselves that the presently obtained (caloric) ΘDh(∞) values (41) show only moderate deviations (<2.5%) from former theoretical estimates [32, 36].Concerning the high-temperature dependencies of theeffective (anharmonicity-related) Debye temperatures, ΘD(T), we would still like to point out that, owing to the strict limitation of their definition (see (1)) to ratios Cp(T)/CVh(∞)≤1, their magnitudes are automatically dropping to zero at those positions, Tf, where the individual Cp(T) curves are crossing the classical Delong-Petit level; that is, (43)ΘDTf=0forCpTf=CVh∞=6R.This basic property is manifested by the rapidΘD→0 drops of the solid ΘD(T) curves (in Figures 4 and 5), which are located just at those positions, Tf = 401 K (for ZnSe) or 321 K (for ZnTe), where the respective Cp(T) curves are crossing the classical Dulong-Petit level, CVh(∞)=6R (as shown in Figures 2 and 3).


## 6.4. Characteristic PDOS-Related Phonon Energy Values

Within the frame of the presently considered four-oscillator hybrid model, the totality of quantitative information on the material-specific PDOS spectra is given (in Table1) in terms of the adjusted sets of effective Einstein oscillator temperatures, Θn (n = 1 to 4), and the respective weights, wn (n = 1 to 4), including the weights wC1 and wC2 for the associated continuous low-energy curve section components (of Debye and non-Debye type, resp.). The corresponding Einstein peakenergy positions, εn≡kBΘn (n = 1 to 4), are quoted in Table 5. On the basis of these sets of fitted parameter values one can readily calculate a variety of characteristic dispersion-related quantities, some of which are directly comparable with their theoretical counterparts being detectible from theoretically calculated PDOS spectra presented in various former papers [32, 36, 45, 46, 66, 89].Table 5
Phonon energy values for the 4 discrete (Einstein) oscillators,εn=kBΘn (n = 1 to 4, corresponding to the Θn values given in Table 1), the average phonon energies for the TA phonon sections, ε-TA (44), and the whole acoustic phonon energy sections, ε-TA&LA (45), including a series of equivalent (moment-related) phonon energies, εP(m) (due to (46) and (47)).
Material
ZnSe
ZnTe
ε
1
=
ε
T
A
1 (meV)
6.75
5.34
ε
2
=
ε
T
A
2 (meV)
9.43
7.58
ε
3
=
ε
-
L
A (meV)
18.36
14.40
ε
4
=
ε
-
T
O
&
L
O
≡
ε
-
U (meV)
27.10
22.69
ε
-
T
A (meV)
8.48
6.80
ε
-
T
A
&
L
A
=
ε
-
L (meV)
11.77
9.33
ε
P
(
-
5
/
2
) (meV)
10.13
8.18
ε
P
(
-
2
) (meV)
11.78
9.49
ε
P
(
-
3
/
2
) (meV)
13.12
10.58
ε
P
(
-
1
) (meV)
14.44
11.68
ε
P
(
-
1
/
2
) (meV)
15.79
12.83
ε
P
(
0
) (meV)
17.11
13.97
ε
P
(
1
)
=
ε
-
L
&
U (meV)
19.43
16.01
ε
P
(
2
) (meV)
21.17
17.55
ε
P
(
4
) (meV)
23.22
19.36
ε
P
(
6
) (meV)
24.27
20.29
ε
P
(
8
) (meV)
24.89
20.83
ε
P
(
10
) (meV)
25.30
21.18Table 6
The 8 lowest-order expansion coefficients,c2≤2m≤16DR, implied by the unprecedented high-capacity (high-temperature) xD2m-series expansion, κDhmR(xD) (A.11), for the conventional Debye function, κD(xD) (A.1). The same expansion coefficients, c2≤2m≤16DR, apply also the analogous x12m-series expansion, κC1hmR(x1) (A.15), for the Debye-like component, κC1(x1) (A.13), of the low-energy TA phonon section of the lattice heat capacity. With respect to the latter we give here also the set of the 8 lowest-order expansion coefficients, c2≤2m≤16NR, which are implied by the analogous x12m-series expansion, κC2hmR(x1) (A.23), for the associated non-Debye heat capacity component, κC2(x1) (A.17).
m
c
2
m
D
R
c
2
m
N
R
1
1
20
=
0.05
5
84
=
0.05952380952
2
1
1400
=
7.14285714286
×
1
0
-
4
13
10584
=
1.228269085
×
1
0
-
3
3
1
648000
=
1.54320987654
×
1
0
-
6
205
19559232
=
1.048098412
×
1
0
-
5
4
-
73
3492720000
=
-
2.09006161387
×
1
0
-
8
4559
160190110080
=
2.845993424
×
1
0
-
8
5
1
82555200000
=
1.21131073512
×
1
0
-
11
-
503
4485323082240
=
-
1.121435381
×
1
0
-
10
6
5591
5148967824000000
=
1.08584869650
×
1
0
-
12
-
478871
1585247736956083200
=
-
3.020795986
×
1
0
-
13
7
-
1879
212199886080000000
=
-
8.85485866515
×
1
0
-
15
-
48299
202404310545527029760
=
-
2.386261987
×
1
0
-
17
8
563743
11383109222284800000000
=
4.9524518213
×
1
0
-
17
887853331
16576922903367866373734400
=
5.355959827
×
1
0
-
17Table 7
High-accuracy sample values for Debye’s normalized heat capacity shape function,κD(xD) (A.1), which are directly comparable with their somewhat more restricted precedents listed in [6, 47] (up to 6 or 8 significant figures, resp.). The individual high-accuracy κD(xD) values have been calculated here for the interval 0≤xD≤1.8, via the unprecedented high-capacity (high-temperature) Taylor series expansion, κDhmR=7(xD) (A.11) and, for the interval 1.8≤xD≤40, by means of Debye’s conventional low-capacity (low-temperature) expansion, κDlnD=15(xD) (A.4) (note that a simple asymptotic expression, κD(xD)→77.92727286/xD3, applies to xD > 40). For comparisons with numerical results quoted in mathematical standard literature (see, e.g., Table 27.1 in [90]) we have listed here also the respective magnitudes of the conventional Debye function, D3(xD) (for 3-dimensional systems), D3(xD)=κD(xD)+3xD/(exp(xD)-1)/4 (according to (A.2) and (A.3)).
x
D
κ
D
(
x
D
)
D
3
(
x
D
)
x
D
κ
D
(
x
D
)
D
3
(
x
D
)
0.0
1.0000000000
1.0000000000
5.0
0.3686348236
0.1175974118
0.1
0.9995001785
0.9629999405
5.5
0.3132547732
0.09524080272
0.2
0.9980028536
0.9269990483
6.0
0.2655968891
0.07758132473
0.3
0.9955144242
0.8919951865
6.5
0.2250561623
0.06360434235
0.4
0.9920454895
0.8579848069
7.0
0.1908564046
0.05250585097
0.5
0.9876107521
0.8249629690
7.5
0.1621687997
0.04365502115
0.6
0.9822288834
0.7929233677
8.0
0.1381873779
0.03656029567
0.7
0.9759223556
0.7618583673
8.5
0.1181717647
0.03084031600
0.8
0.9687172408
0.7317590428
9.0
0.1014670944
0.02619989259
0.9
0.9606429816
0.7026152269
9.5
0.08750918665
0.02241065587
1.0
0.9517321357
0.6744155641
10.0
0.07582100303
0.01929576569
1.1
0.9420200996
0.6471475686
11.0
0.05773082402
0.01457049734
1.2
0.9315448136
0.6207976880
12.0
0.04478011684
0.01125032746
1.3
0.9203464535
0.5953513709
13.0
0.03534728217
0.008858858804
1.4
0.9084671129
0.5707931370
14.0
0.02835178383
0.007096677015
1.5
0.8959504801
0.5471066514
15.0
0.02307128730
0.005771263228
1.6
0.8828415127
0.5242747994
16.0
0.01901817415
0.004755893959
1.7
0.8691861161
0.5022797638
17.0
0.01585873694
0.003965212076
1.8
0.8550308264
0.4811031024
18.0
0.01336097602
0.003340449610
1.9
0.8404225040
0.4607258256
19.0
0.01136091855
0.002840309476
2.0
0.8254080384
0.4411284737
20.0
0.009740756601
0.002435220067
2.2
0.7943467192
0.4041938106
21.0
0.008414505229
0.002103638250
2.4
0.7622123728
0.3701368830
22.0
0.007318466926
0.001829621334
2.6
0.7293549760
0.3387926100
23.0
0.006404797367
0.001601201112
2.8
0.6961028658
0.3099951582
24.0
0.005637096943
0.001409274915
3.0
0.6627580442
0.2835798281
25.0
0.004987344231
0.001246836318
3.2
0.6295928776
0.2593847028
26.0
0.004433731487
0.001108432971
3.4
0.5968481116
0.2372520497
28.0
0.003549893920
0.0008874734947
3.6
0.5647320637
0.2170294741
30.0
0.002886195280
0.0007215488222
3.8
0.5334208309
0.1985708314
32.0
0.002378151635
0.0005945379091
4.0
0.5030593309
0.1817369138
34.0
0.001982680461
0.0004956701152
4.2
0.4737629932
0.1663959309
36.0
0.001670251904
0.0004175629760
4.4
0.4456199253
0.1524238050
38.0
0.001420164616
0.0003550411541
4.6
0.4186933946
0.1397043074
40.0
0.001217613638
0.0003044034095
4.8
0.3930244869
0.1281290579Table 8
Fitted values of the model-specific coefficients involved by the refined low-to-intermediate capacityxD(κ) formula, xDl9(0<κ≤κ9) (B.6), and for the couple of complementary high-accuracy low- and high-capacity Debye temperature formulas, xDl12(0<κ≤κc) (B.8) and xDh11(κc≤κ≤1) (B.9), which have been determined by least-mean-square fittings of a sufficiently comprehensive set of couples of high-accuracy κ versus xD curve points (in accordance with Table 7).
x
D
l
9
(
κ
) (see (B.6))
x
D
l
12
(
κ
) (see (B.8))
x
D
h
11
(
κ
) (see (B.9))
0
<
κ
≤
κ
9
=
0.8844
0
<
κ
≤
κ
c
=
0.260
0.260
=
κ
c
≤
κ
≤
1
p
=
0.071797185
a
2
=
2.8743354
×
10
-
1
b
2
=
6.0777673
×
10
-
3
g
1
=
3.042348
a
3
=
1.1377662
×
10
-
1
b
3
=
5.1466371
×
10
-
3
g
2
=
1.816742
a
4
=
-
9.5999179
×
10
-
2
b
4
=
-
6.8442924
×
10
-
3
g
3
=
-
6.892814
a
5
=
4.3062699
×
10
-
2
b
5
=
2.6051894
×
10
-
2
g
4
=
6.954053
a
6
=
-
4.1060421
×
10
-
3
b
6
=
-
3.6011454
×
10
-
2
g
5
=
5.456090
a
7
=
-
1.1955013
×
10
-
3
b
7
=
3.6169365
×
10
-
2
g
6
=
-
11.20018
a
8
=
4.9953204
×
10
-
4
b
8
=
-
3.7204523
×
10
-
2
g
7
=
2.022044
a
9
=
-
6.2681506
×
10
-
5
b
9
=
1.6804524
×
10
-
2
g
8
=
2.824376
a
10
=
3.3449835
×
10
-
6
b
10
=
-
3.0967328
×
10
-
3
—
p
=
1.2969
b
11
=
5.2219439
×
10
-
2
—
g
2
=
-
3.1541530
b
12
=
-
6.2935315
×
10
-
2
—
g
3
=
8.7184026
—Table 9
High-precision sample values of the characteristic Debye variable,xD(≡ΘD(T)/T), which have been calculated for a series of discrete (equidistant) κ-values via the couple of complementary low- and high-capacity formulas, xDl12(0<κ≤κc) (B.8) and xDh11(κc≤κ≤1) (B.9), using the model-specific coefficients listed in Table 8 (note that we have obtained precisely the same xD(κ) values also in alternative way via numerical integration procedures based directly on (A.1), for integration steps chosen to be of order 10-7).
κ
x
D
(
κ
)
κ
x
D
(
κ
)
κ
x
D
(
κ
)
κ
x
D
(
κ
)
κ
x
D
(
κ
)
0.01
19.825638
0.21
6.709531
0.41
4.666635
0.61
3.319275
0.81
2.100219
0.02
15.733335
0.22
6.568697
0.42
4.590077
0.62
3.258266
0.82
2.035428
0.03
13.737134
0.23
6.434348
0.43
4.514891
0.63
3.197532
0.83
1.969683
0.04
12.468384
0.24
6.305852
0.44
4.440994
0.64
3.137032
0.84
1.902849
0.05
11.557161
0.25
6.182662
0.45
4.368306
0.65
3.076725
0.85
1.834773
0.06
10.854261
0.26
6.064298
0.46
4.296752
0.66
3.016570
0.86
1.765274
0.07
10.285969
0.27
5.950343
0.47
4.226263
0.67
2.956525
0.87
1.694144
0.08
9.810966
0.28
5.840425
0.48
4.156772
0.68
2.896548
0.88
1.621130
0.09
9.403982
0.29
5.734216
0.49
4.088217
0.69
2.836595
0.89
1.545932
0.10
9.048539
0.30
5.631424
0.50
4.020538
0.70
2.776622
0.90
1.468177
0.11
8.733353
0.31
5.531788
0.51
3.953679
0.71
2.716581
0.91
1.387399
0.12
8.450382
0.32
5.435072
0.52
3.887587
0.72
2.656425
0.92
1.302998
0.13
8.193715
0.33
5.341064
0.53
3.822209
0.73
2.596103
0.93
1.214179
0.14
7.958888
0.34
5.249572
0.54
3.757496
0.74
2.535560
0.94
1.119852
0.15
7.742451
0.35
5.160421
0.55
3.693400
0.75
2.474740
0.95
1.018444
0.16
7.541692
0.36
5.073452
0.56
3.629874
0.76
2.413582
0.96
0.907538
0.17
7.354438
0.37
4.988520
0.57
3.566875
0.77
2.352021
0.97
0.783057
0.18
7.178928
0.38
4.905490
0.58
3.504357
0.78
2.289985
0.98
0.637031
0.19
7.013711
0.39
4.824241
0.59
3.442279
0.79
2.227398
0.99
0.448821
0.20
6.857583
0.40
4.744657
0.60
3.380599
0.80
2.164173
1.00
0.000000Table 10
Limiting values of the “true” (harmonic lattice) Debye temperatures,ΘDh(∞), including adjusted magnitudes of several even-order expansion coefficients, a2s, that are representing the weights of the individual (ΘDh(∞)/T)2s power terms in (C.1). These sets of parameter values have been determined via fittings (shown in Figure 10) of an appropriate selection of usable Debye temperature points, ΘDh(T)≅ΘD(T), implied by the experimental Cp(T) data points given by Kremer et al. [46] for natural ZnSe and ZnTe. Listed are also the magnitudes of the respective (Debye-model-related) cut-off phonon energies, εDh(∞)≡kBΘDh(∞), and several even-order phonon energy values, εP(2s)≡(μP(2s))1/2s (for 2s=2 to 10), which are resulting from the corresponding PDOS spectra moments, μP(2s), calculated on the basis of the analytical expressions (C.2a) to (C.2e) (note that the latter equations are exactly corresponding to the series of μP(2s) formulas originally derived within the frame of the supporting information to [21]; see also appendix B of [14] and/or the appendix B of [15]).
ZnSe
ZnTe
Intervals (K)
17 to 96
15 to 80
Θ
D
h
(
∞
) (K)
317.9
262.3
a
2
0.9259 × 10−2
0.9995 × 10−2
a
4
0.6895 × 10−4
0.7851 × 10−4
a
6
0.2869 × 10−6
0.3410 × 10−6
a
8
0.6533 × 10−9
0.7714 × 10−9
a
10
(0.7135 × 10−12)
(0.7038 × 10−12)
a
12
(0.2510 × 10−15)
—
ε
D
h
(
∞
) (meV)
27.40
22.60
ε
P
(
2
) (meV)
21.22
17.51
ε
P
(
4
) (meV)
23.48
19.45
ε
P
(
6
) (meV)
24.83
20.60
ε
P
(
8
) (meV)
25.75
21.38
ε
P
(
10
) (meV)
26.42
21.95Of importance for a physical understanding of the prominent role of short-wave TA phonons, which are mainly responsible for the occurrence of the characteristic non-Debye (maximum) behaviours ofCp(T)/T3 curves in the liquid hydrogen region (cf. the insets of Figures 2 and 3), is also the knowledge of the centres of gravity of TA phonon sections, ε-TA. The magnitudes of the latter are given, in accordance with the respective TA phonon normalization condition (see (23)), by [14](44)ε-TA=33wC14+5wC26+w1ε1+w2ε2 (cf. Eq. (18) in [14]). Concerning the presently estimated ε-TA values of 8.48 meV for ZnSe and 6.80 meV for ZnTe (as quoted in Table 5), we can make the satisfying statement that they are in good agreement with their counterparts (of about 8.4 meV or 8.5 meV, for ZnSe, and 6.8 meV or 6.9 meV, for ZnTe), which we have detected from the theoretically calculated PDOS spectra presented for ZnSe and ZnTe in Figure  3 of [45] or Figure  4 of [46], respectively. Accidentally, these ε-TA positions are even nearly coinciding with the positions of the highest TA peaks shown within the same figures (which are readily seen to be located at about 8.5 meV or 8.6 meV, for ZnSe, and between 6.8 meV and 6.9 meV, for ZnTe).Similarly informative are also comparisons of the average values of LA phonon energies,ε-LA, which are represented within the present formalism (due to (22)) by the third discrete oscillator, ε-LA=ε3. The presently estimated ε-LA values of 18.36 meV for ZnSe and 14.40 meV for ZnTe (cf. Table 5) are in reasonable agreement with their counterparts (of about 18.5 meV or 19.2 meV, for ZnSe, and 14.5 meV or 15.0 meV, for ZnTe), which we have detected, again, from the theoretically calculated PDOS spectra presented in Figure  3 of [45] or Figure  4 of [46], respectively.The average values for all acoustical (i.e., TA and LA) phonon energies,ε-TA&LA, are another characteristic quantity for comparisons between the results of numerical analyses of heat capacity data with theoretically calculated PDOS spectra. The analytical expression for the latter, which thus comprises the wholelower halves of PDOS spectral sections, is given within the present formalism by [14](45)ε-L=ε-TA&LA=23wC14+5wC26+w1ε1+w2ε2+ε36=2ε-TA+ε-LA3.The presently estimatedε-L=ε-TA&LA values of 11.77 meV for ZnSe and 9.33 meV for ZnTe (cf. Table 5) are found again to be in reasonable agreement with their counterparts of about 11.8 meV or 12.1 meV, for ZnSe, and 9.4 meV or 9.6 meV, for ZnTe, which followed immediately via (45) from the couples of ε-TA and ε-LA values for the average TA and LA phonon energies (which had been detected from the theoretically calculated PDOS spectra presented in Figure  3 of [45] or Figure  4 of [46], resp.). Furthermore we observe that the presently estimated (caloric) ε-L=ε-TA&LA values are in good agreement with our former estimations [56], which were based on theoretically calculated PDOS spectral functions published in various earlier papers, namely, 11.7 meV due to [66, 89], or 12.5 meV due to [36], for ZnSe, and 9.5 meV or 9.4 meV due to [32, 36], for ZnTe.Let us still consider theupper part of the PDOS spectra in consideration, which are consisting of one LO and two TO branches. The corresponding average optical phonon energy is given within the present formalism (due to (21)) just by the fourth discrete oscillator; that is, ε-U=ε-TO&LO=ε4. The presently estimated ε-U values (cf. Table 5) of 27.10 meV, for ZnSe, and 22.69 meV, for ZnTe, are found to be in more or less good agreement with their counterparts (of about 27.3 meV or 25.6 meV, for ZnSe, and 23.1 meV or 21.7 meV, for ZnTe), which we have detected again from the theoretically calculated PDOS spectra presented in Figure  4 of [46] and Figure  3 of [45], respectively. Informative is also a comparison of present (caloric) results with our former estimations [56] of ε-U values on the basis of PDOS spectral functions published in various earlier papers, namely, 27.3 meV, 27.1 meV, or 28.4 meV due to [36, 66, 89], for ZnSe, and 22.0 meV or 22.7 meV due to [32, 36], for ZnTe. Thus we are concerned in the case of the upper (optical) PDOS spectral sections with differences between approximate ε-U values of different origin up to a maximum order of 5% (for both ZnSe and ZnTe). At the same time we can see that the presently estimated (caloric) ε-U values (Table 5) are ranging near the middle between the highest and lowest ε-U estimates referring to the variety of PDOS spectra in consideration for the two materials under study. This shows that one has good chances to get reliable information on various prominent features of material-specific PDOS spectra even via dispersion-related analyses of heat capacity data.Very informative in this respect are, above all, comparisons between momentsμP(m), of material-specific PDOS spectral functions [3, 4, 16, 18, 19, 21, 50, 52, 56, 69–72], which are generally defined by (12) (for -3<m≠0). Within the present formalism of the four-oscillator hybrid PDOS model function, gP(ε) (i.e., (7), for nE = 4), one obtained for the corresponding moments, μP(m) (see (12)), the model-specific expression (13) [14, 21]. Yet, for visualizations of the m-dependencies of such PDOS spectra moments it is convenient to represent the individual mth-order moments, μP(m), in terms of their equivalent mth-order phonon energies [14, 16, 21, 50, 56], εP(m)≡μP(m)m. The latter are thus given within the frame of the present four-oscillator hybrid model (in accordance with (13)) by an expression of the form [14, 21, 50](46)εPm=ε13wC13+m+5wC25+m+w1+∑n=24wnεnε1m1/m,for-3<m<0,0<m.By means of am→0 limiting transition one obtains from (46) in particular(47)εPm⟶0=ε1exp-wC13-wC25+∑n=24wnlnεnε1.Collections of corresponding material-specificεP(m) values for a series of integers and half-integers (within the interval -3<m≤10) are listed in Table 5 and are represented by open circles in Figure 7. In view of the fact that the whole m-dependence of the εP(m) curve (due to (46)) divides into a couple ofcontinuous curve sections (for -3<m<0 and 0<m), the sequence of discrete εP(m) values (○) are perfectly connected by the respective solid curve sections. Furthermore we have shown in Figure 7 several sets of discrete εP(m) values which we have obtained via detailed numerical analyses of the theoretically calculated PDOS spectra presented for ZnSe and/or ZnTe in [32] (×), [36] (∗), [45] (◊), [46] (△), [66] (+), and [89] (☆). The comparisons between the different εP(m) values, in Figure 7, show a multitude of good agreements between the discrete εP(m) values due to the PDOS spectra published by different authors, including present (caloric) results (○). This concerns especially the region of negative orders (m<0), where the actual εP(m) values are predominantly due to the contributions of TA phonons. At the same time we are concerned with some markedly larger deviations in the region of positive orders (m > 0), where the monotonically increasing εP(m) values tend to be the stronger ones governed by optical phonons the higher their orders are. The maximum difference between estimated εP(10) values for ZnTe (due to [46] (△) versus [45] (◊)) amounts to about 6%, whereas, for ZnSe, the maximum difference (due to [89] (☆) versus [45] (◊)) amounts even to about 13%. Relatively large differences like the latter ones are obviously mainly due to significant qualitative differences between the upper (optical) sections of the theoretically calculated PDOS spectra presented by the corresponding authors (of [45] (◊) versus [89] (☆)). At the same time we can make the satisfactory statement that the presently obtained (caloric) εP(10) values (Table 5) differ from those we have derived from several other PDOS spectra in consideration for ZnSe ([36] (∗), [46] (△), and [66] (+)) and for ZnTe ([32] (×), [36] (∗), and [46] (△)) by less than 2%.Figure 7
Visualization of the effective phonon energy curves,εP(m>-3) (see (46); solid curves), for ZnSe and ZnTe. Included are the corresponding discrete εP(m) points ((○); cf. Table 5), for integer and several negative half integer orders. For comparisons with counterparts following from former theoretically calculated PDOS spectra we have shown here also the respective series of discrete εP(m) values obtained via detailed numerical analyses of the theoretically calculated PDOS spectra presented for ZnSe and/or ZnTe in [32] (×), [36] (∗), [45] (◊), [46] (△), [66] (+), and [89] (☆). The series of even-order εP(m) values (for m=2 to 10, listed in Table 10), which have been determined in Appendix C via the conventional analyses procedure for Debye temperature values representing the harmonic lattice, ΘDh(T)≅ΘD(T) (cf. Figure 10), are represented by squares (□).The latter findings underscore the above assessment, according to which one has as a rule good chances to get reliable information of various characteristic features of material-specific PDOS spectra by means of detailed dispersion-related analyses of heat capacity data. The rather fine quantification abilities of our elaborate analytical framework (displayed in Section4) concern obviously also such delicate quantities likehigher-order phonon energy values, εP(m≫1), and, consequently, even the corresponding higher-order moments, μP(m)=(εP(m))m, of PDOS spectra.Just with respect to the latter it might be instructive to consider still the alternative series of five even-orderεP(m) values (for m = 2 to 10), which have been determined in Appendix C by means of the conventional analyses procedure [14–18, 21, 71, 88] for Debye temperatures representing the harmonic lattice properties, ΘDh(T)≅ΘD(T) (the range of which is limited to T < 96 K, for ZnSe, and to T< 80 K, for ZnTe; cf. Figure 10 and Table 10). The corresponding alternative εP(m) values (listed in Table 10) are represented in Figure 7 by squares (□). Comparing the latter with their hybrid-model-related counterparts (○) (listed in the last five rows of Table 5) we find excellent agreements, at least, for the lowest two orders (m = 2 and 4). Towards higher orders (m = 6 to 10), however, we are concerned with growing deviations between the two alternative series of εP(m) values, which are seen to increase for ZnSe from 2.3% to 4.4% and for ZnTe from 1.5% to 3.6% (cf. Figure 7).Finally, let us still consider the behaviour (m-dependence) of the conventionalDebye-model-related energy curves, which are defined (in analogy to [16–19, 64, 71]) by an expression of the form(48)εDPm≡3+mμPm31/m=3+m31/mεPm. in accordance with the general definition, εP(m)≡μP(m)1/m, of the mth-order phonon energies [14, 21, 50, 56]. Observing that the m-dependencies of the two complementary εP(m) curve sections (for -3<m<0 and 0<m) are given by (46), we come thus readily (via (48)) to the material-specific εDP(m) dependencies, for ZnSe and ZnTe, we have visualized by corresponding dash-dot curves in Figure 7 (note that, for the special point m→0, the particular εDP(0) value follows from a corresponding limiting transition to adopt just the value εDP(m→0)=e1/3εP(0) (cf. [16, 17]), where εP(0) is given by (47)).A prominent feature of these (conventional) Debye-model-related energy curves,εDP(m) (see (48)), is due to the well-known property [16, 17, 19], according to which the m→-3 limit of this function is automatically (exactly) equal to the T→0 limit of the cut-off energy for Debye’s fictive phonon energy spectrum; that is,(49)εDPm⟶-3=wC1-1/3ε1=kBΘD0=εD0, in accordance with (20). The corresponding magnitudes of εD(0)=23.3 meV, for ZnSe, and εD(0) = 19.0 meV, for ZnTe, are indicated by solid triangles (▸) in Figure 7.Another characteristicεDP(m) value is that for m = 2, which follows readily from (48), in combination with (42), to be just equal to the T→∞ limit of the cut-off energy for Debye’s fictive (harmonic) phonon energy spectrum; that is,(50)εDPm=2=532εP2=kBΘDh∞≡εDh∞.The corresponding magnitudes ofεDh∞=27.3 meV, for ZnSe, and εDh(∞) = 22.6 meV, for ZnTe, are indicated by solid squares (■) in Figure 7. Accordingly, the ratios between these two limiting Debye energy (or temperature) values, εD(0)/εDh(∞)=ΘD(0)/ΘDh(∞), amount to about 0.85 or 0.84 for ZnSe and ZnTe, respectively (cf. Table 4). The significant differences (of about 15%) between these two limiting Debye temperature values are already, per se, clear indications for a pronouncednon-Debye character of both PDOS spectra in consideration.Actually, for the hypothetical case of a properly normalized phonon energy spectrum of Debye’s type, that is, for the fictive (exclusively quadratic) PDOS functiongD(ε)=3ε2/εDc3 [14, 21] (within the interval 0≤ε≤εDc, and gD(ε)=0, for εDc<ε<∞), the corresponding moments are readily seen to be given just by μD(m)=3εDcm/(3+m) [14]. Substituting hence the actual (material-specific) moments μP(m) in (48) by the moments μD(m) due to Debye’s fictive model, we see readily that expression (48) for εDP(m) reduces just to a constant(51)εDPm⟶εDDm≡3+mμDm31/m=εDc,for  any  order-3<m<∞.Such a constancy of anεDP(m) function (of type (48)), however, has never been observed in physical reality (not even in form of a relatively close approach of one or the other εDP(m) curve to a horizontal line). The rather strong variations of εDP(m) values shown in Figure 7, particularly the pronouncednonmonotonic (minimum) behaviours of both εDP(m) curves in the vicinity of m≈-2, can be looked upon as clear visual manifestations (and a global confirmation) of the obviously very pronouncednon-Debye character of the PDOS spectra [32, 36, 45, 46, 66, 89] for both materials under study.


## 7. Conclusions

The most important progress of the present study in comparison with all earlier papers, which involved in one or the other way the somewhat delicate concept of material-specific Debye temperatures,ΘD(T), consists obviously in our elaboration (in Appendix B) of appropriate integral-free formulas for direct calculations of ΘD(T) values associated with given heat capacities, Cp(T). In view of the sample character of the present applications of this novel analytical apparatus to given isobaric and/or isochoric (harmonic) heat capacity data sets, we have performed the corresponding Cp/Vh(T)⇒ΘD/Dh(T) transformations here (in Section 5) on the basis of a couple of high-accuracy formulas (27) and (28), the residual deviations from exact Debye temperature values can be seen (from the lower inset of Figure 9) to be extremely small (limited to an order of ±0.000016%). Such a high degree of accuracy, however, did not use to be actually required in practice, owing to the well-known fact that typical experimental uncertainties of measured heat capacities used to be of order ±1%. Couples of somewhat simpler low- and high-capacity formulas might thus also be largely adequate for future practical applications. In this respect we would like to draw attention above all to the somewhat simpler low-capacity formulas of type (B.5a) and (B.5b) or (B.6) which, in combination with the complementary high-capacity formula (B.3), are providing already good results (within maximum deviations of orders ±0.13% or ±0.015%, resp.; cf. Figure 9). In particular the latter choice appears to be quite sufficient for forthcoming numerical applications.A crucial prerequisite for our analytical constructions and tests of appropriate analytical formulas suitable forCp/Vh(T)⇒ΘD/Dh(T) transformations (within Appendix B) consisted in a preliminary, high-accuracy determination of the general dependence of heat capacities, Cp/Vh(T), on given Debye temperatures, ΘD/Dh(T). An excellent solution of this conventional problem has been found in Appendix A in form of a couple of an unprecedented high-capacity series expansion (see (A.11)) and Debye’s conventional low-capacity (low-temperature) representation [1–3] (see (A.4)). We have shown (in Appendix A) that it is rather easy to perform high-accuracy ΘD/Dh(T)⇒Cp/Vh(T) transformations by means of this couple of integral-free expressions. A collection of corresponding high-accuracy sample values (consisting of 10 significant figures) for Debye’s model-specific κD(xD) function (see (A.5)) is given in Table 7. In view of this extraordinarily comfortable way offered here (in Section A.1 of Appendix A) for high-accuracy calculations of heat capacities, Cp(T)≡κD(T)CVh(∞), at given magnitudes of Debye temperature versus lattice temperature ratios, xD(T)=ΘD(T)/T, it will hence no longer be necessary either to perform point-by-point calculations of κD(xD) values (via numerical integrations according to (A.5)) or to resort to cumbersome point-by-point interpolations between neighbouring κD(xD) values utilizing an available table with discrete κD(xD) values (like Table 7 in the present paper or one of those published in [6, 8, 47–49]).Another significant progress with respect to the intended elaboration of a more easy-to-use analytical apparatus for heat capacity analyses within the frame of the hybrid model (Section3) has been achieved here (in Section A.2 of Appendix A) by the derivation of couples of appropriate high- and low-capacity formulas for the normalized low-energy TA phonon contributions of Debye and non-Debye type, κC1(Θ1/T) and κC2(Θ1/T) (see (10)), to the measured cryogenic heat capacities. Incorporating the corresponding couples of high- and low-temperature expressions (i.e., (14) to (17)) into the general formalism of the multioscillator hybrid model we have established a completely integral-free version of our analytical framework. Its actual use effectuated a significant acceleration of the multiparameter least-mean-square fitting processes performed (in Section 4) for fine simulations of the temperature dependencies of the comprehensive heat capacity data sets under study.An additional important progress in quantitative respect (in comparison with our previous study [14]) could be made here owing to the availability of the unusually comprehensive and informative low-temperature heat capacity data sets given by Kremer et al. [46] for natural ZnSe and ZnTe. Of great importance just for ZnTe have been in this connection above all the unprecedented cryogenic Cp(T) data available for the liquid-helium-hydrogen region (2 K <T < 20 K), which gave us here for the first time the chance to perform a realistic estimation of thecaloricT→0 limiting value of the Debye temperature, ΘD(0)≅ 220 K (±2 K). Moreover, for the case of ZnSe, the novel Cp(T) data implied a noticeable refinement (in comparison with [14]) of the ΘD(T) behaviour in the region between 100 K and 300 K (cf. Figure 4 in the present paper with Figure  12 in [14]).Finally we would still like to point out that, in comparison with a possible use of the global (apparently universal) non-Debye formula displayed in [15], a decisive advantage of the multioscillator hybrid model used within the present study consists in the delivery of physically important byproducts in form of a multitude of adequate quantitative information (like those presented in Section 6.4) on the PDOS functions of the materials under study. At the same time one must make the reservation that the hybrid model in the specializedfour-oscillator version (as displayed in Section 4) used to function quite well just for systems, the PDOS spectra of which are consisting of 3 acoustical in combination with only 3 optical phonon branches. This favourable state of affairs is encountered above all within a large variety of III–V and II–VI materials with cubic structure, including silicon and germanium. However, for systems with PDOS spectra comprisingmore than 3 optical branches (e.g., for binary systems with hexagonal structure, or ternary systems), it will as a rule be necessary to envisage more or less different constellations of discrete oscillators in order to come again to fine numerical simulations of corresponding heat capacity data sets.



---

*Source: 9321439-2017-09-18.xml*