A systematic study of the negative thermal expansion in zinc-blende and diamond-like semiconductors

Upon heating, almost all zinc-blende (ZB) and diamond-like semiconductors undergo volume contraction at low temperature, i.e. negative thermal expansion (NTE), instead of commonly expected expansion. Specifically, CuCl has the largest NTE among these semiconductors with a coefficient comparable with the record value of ZrW2O8. So far, underlying physical mechanism remains ambiguous. Here, we present a systematic and quantitative study of the NTE in ZB and diamond-like semiconductors using first-principles calculations. We clarified that the material ionicity, which renders the softening of the bond-angle-bending and thus, the enhancement of excitation of the transverse acoustic (TA) phonon, is responsible for the NTE of ZB and diamond-like semiconductors. With the increase in the ionicity from the groups IV, III-V, IIB-VI to IB-VII ZB semiconductors, the coefficient of the maximum NTE increases due to the weakness in bond-rotation effect, which makes the relative motion between cation and anion transverse to the direction of the bond more feasible and the mode Grüneisen parameters of the TA modes more negative. Since CuCl has the highest ionicity among all ZB and diamond-like semiconductors, it is expected to have the largest NTE, in good agreement with the experimental observation. This understanding would be beneficial for tetrahedral materials with specific applications.


Introduction
Material undergoing expansion upon heating and contraction upon cooling is a common expectation. However, some materials are discovered to demonstrate a contraction upon heating, known as negative thermal expansion (NTE) [1][2][3][4][5][6][7]. Such NTE materials have attracted considerable attention for use as thermal-expansion compensators by mixing with a 'normal' material which expands on heating to fabricate composites with tailored expansivity, which have a range of potential engineering, mechanical, ferroelectric and chemical applications [8][9][10][11][12][13]. Interestingly, almost all common semiconductors, which are dominant materials in the microelectronic and optoelectronic industries, were found experimentally to have NTE at low temperature (see figure 1). The linear coefficients α L of NTE for conventional zinc-blende (ZB) and diamond-like semiconductors, such as GaAs, GaSb, ZnS, ZnSe, Si, and Ge, reach α L =−1×10 −6 K −1 [14]. Among all ZB semiconductors, cuprous chloride (CuCl) possesses the largest NTE of α L =−8.33×10 −6 K −1 , [15,16] which is comparable with the historically record value discovered in ZrW O 2 8 [17-19], and only one order smaller than that of the recently reported giant NTE in Ca RuO 2 4 [5]. Because the rapid advances in modern semiconductor technology require precise control of thermal expansion for fabrication of heterostructure semiconductor devices, understanding the underlying mechanisms of NTE in ZB and diamond-like semiconductors is fundamentally important.
So far, several dominant categories of physical mechanisms have been proposed to explain the NTE: [10,20]  Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. transition, and (v) bond and associated vibrations. Mechanism (i) has been proposed to explain the NTE of ZrW O 2 8 , where the rotation of the two-linked tetrahedral WO 4 causes significant distortions within the octahedral ZrO 6 rigid units, which consumes the open spaces in the crystal [21][22][23]. Mechanism (ii) considers high-spin state usually being relatively larger in radius than the low-spin state for an atom with given valence and coordination number [24][25][26]. Mechanism (iii) is the volume change involved with the variation of magnetic moment such as in Fe-36Ni [27]. Mechanism (iv) involves some change in symmetry or space group of a crystal [28][29][30]. As to the mechanism (v), earlier, Barron et al [1,15,20,22,[31][32][33] argued qualitatively that the tension effect is responsible for the NTE in ZB and diamond semiconductors. Upon heating, transverse motion causes the bond to shorten and decreases volume in solids. Nevertheless, a quantitative and systematic consideration of GaZ, (f 1 )-(f 3 ) InZ, and (g 1 )-(g 3 ) Si, Ge, Sn, respectively, where Xä(Cl, Br, I), Yä(S, Se, Te) and Zä(P, As, Sb). the NTE occurred in ZB-and diamond-structure semiconductors remains lacking and why the simple tetrahedral CuCl exhibits extremely large NTE, which are observed usually in complex structure materials [10,11,34].
In this paper, we present a systematic, comparative, and quantitative study of the NTE in ZB and diamondlike semiconductors at low temperature. We clarify how the material ionicity which renders the softening of bond-angle-bending determines the NTE of those semiconductors. We demonstrate that the maximum NTE coefficients increase significantly as the ionicity getting stronger from groups IV, III-V, IIB-VI to IB-VII semiconductors because the reduced bond-angle-bending force constant causes the enhancement of the vibrational magnitudes in the transverse acoustic (TA) phonon modes. Since CuCl has the highest ionicity and smallest bond-angle-bending rigidity, it is thus expected that CuCl has the largest NTE among all ZB and diamond-like semiconductors, in agreement with experimental data [15].

Method of calculations
Our calculations are performed within the first-principles density-functional theory (DFT) as implemented in quantum-espresso package [35,36]. The generalized-gradient approximation which was parametrized in the form of Perdew-Burke-Ernzerhof for exchange-correlation potentials is adopted [37]. The norm-conserving pseudopotentials for treating the valence electrons are used [38,39]. To determine the static electronic structural properties, we minimize the total energy of the system, where the kinetic energy cut-off for the plane wave basis set is chosen as 120Ry. The energy threshold for convergence is 1×10 −12 Ry, and the force threshold for each atom is 1×10 −10 Ry/Bohr. As to the self-consistent field calculations, the k-points in the Brillouin-zone (BZ) are sampled by a 24×24×24 Monkhorst-Packmesh, [40] while for the phonon calculations, we set the q-points as in a 6×6×6 mesh and then a very dense sampling is obtained using a Wannier interpolation.
To study the linear thermal expansion coefficient α L (T), we employ the thermodynamic theory as implemented in thermo-pw package [41]. We consider the anharmonic effect, in which the volume dependence of the phonon frequencies in a solid at finite temperature is included. We first take derivatives of the Helmholtz free energy with respect to strain that gives the stress acting on the solid. This yields an equation of state which correlates strain, stress and temperature. Through searching the strain for which the stress is zero, we determine the equilibrium geometry of the solid at each temperature as the minimum of the Helmholtz free energy. Using this information we calculate phonon dispersions and finally deduce the thermal expansion of a solid.

Results and discussion
3.1. First-principles determination of thermal expansion coefficient figure 1 shows the first-principles calculated thermal expansion coefficients α L (T) as a function of temperature for all investigated ZB and diamond-like semiconductors. We find that all investigated materials undergo volume contraction and expansion in the temperature range from 0 K to room temperature. As increasing temperature, α L (T) decreases initially from zero to negative giving rise to NTE, then starts to increase through archiving the minimum (which is denoted by a L M ), and finally becomes positive giving rise to positive thermal expansion (PTE) for all the investigated semiconductors. Figure 1 also shows that our predicted α L (T) is in reasonably good agreement with the available experimental data [14], especially at low temperature. The corresponding first-principles predicted a L M is summarized in table 1 in comparison with experimental values, if the latter are available. For the cuprous halides, we find that the calculated a = -´--3.96 10 K L M 6 1 (occurring at T=34 K) is indeed the largest in magnitude among all ZB and diamond-like semiconductors as observed experimentally. However, we should note that its magnitude is still smaller than the experimental measured value of a = -´--8.33 10 K L M 6 1 . This discrepancy is attributed to the underestimation of the occupied d orbitals of e.g. Cu, in DFT calculations, which artificially leads the wave functions of d electrons to be more delocalized, and hence stiffens the TA phonon modes as it reduces the materials ionicity (see detailed discussions below) [42][43][44].

Mode-dependent Grüneisen parameters
To understand why the NTE commonly happens in ZB and diamond-like semiconductors, we first examine their lattice vibrational properties. Generally speaking, the longitudinal vibration for longitudinal acoustic (LA) modes is along the bond direction which favors the crystal volume expansion, whereas the transverse vibration in the TA modes is associated with the bond-angle-bending (perpendicular to the bond) which tends to reduce the crystal volume [15,20,22,45,46]. This effect can be identified from the dimensionless mode Grüneisen parameter: [41,47,48] (ν is phonon mode, q is wavevector, and V is volume), which measures how sensitive the vibrational eigenfrequencies ω q,ν are to the variation of the volume. Hereby, the NTE coefficient α L (T) can also be expressed as a , , [21,41,49] where B is the bulk modulus and n C q, the heat capacity. It is obvious that a negative α L (T) is exclusively arising from negative γ q,ν . Because expansion of a compound upon heating depends on the relative contributions of all phonons, a feature necessary for the occurrence of NTE is the presence of low energy phonons with negative Grüneisen parameters, and a phonon gap that separates these modes from the high energy phonons in the structure. If there are enough modes with negative γ q,ν excited to outweigh those with positive γ q,ν , the thermal expansion is negative.
The mode Grüneisen parameters γ q,ν as a function of q for both CuCl and GaAs are shown in figures 2(a) and (b), respectively. One can see that, in CuCl, γ q,ν of all TA modes are negative in the whole BZ except in the vicinity of K, in sharp contrast to the LA modes as well as the optical modes in which γ q,ν are positive through out the whole BZ. To have large NTE, one should excite the TA modes as much as possible but suppress the LA modes as well as optical modes. Indeed, as increasing temperature, the low-laying TA modes are first excited before the LA and optical modes. Therefore, the NTE of semiconductors occurs at low temperature. Whereas, in GaAs, γ q,ν of TA modes are negative only in the vicinity of X and are much smaller in magnitude than that of CuCl. Consequently, GaAs has a much small NTE coefficient compared to CuCl as appears in table 1.

Vibrational density of states and phonon excitation
To examine the phonon gap that separates the TA modes from the remaining modes with positive γ q,ν , we compare in figure 2(c) the vibrational density of states (VDOS), g(ω) between GaAs and CuCl as a function of frequency ω for TA and LA modes, respectively. In the case of CuCl, g TA (ω) exhibits a high peak at w = 4.45 TA P meV (corresponding to a temperature T = 52 K), and g LA (ω) has a peak at w = 13.65 with n(ω, T) being the Bose-Einstein distribution function. It is readily noted that when T<100 K, N TA is one order of magnitude larger than N LA in both CuCl and GaAs, confirming the lattice vibration is excited dominantly by the TA phonons in the low temperature regime. Furthermore, N TA of CuCl is much larger than that of GaAs at low temperature, contributing more to NTE in CuCl compared with GaAs. Therefore, to have a large NTE coefficient for a material, w TA P should be as small as possible while w LA P has a large energy separation from w TA P , so that more TA modes with negative γ q,ν and less LA modes with positiveγ q,ν are excited at a specific temperature. Table 1. We present the maximum negative value of a L M ( --10 K 6 1 ) from both theoretical calculations and experimental measurements [14] for all the investigated ZB and diamond-like semiconductors; T c is the corresponding transition temperature, where experimental data of a L for CuBr, CuI, CdS and CdSe are not available yet in the literature; w TA P and w LA P are the phonon frequencies (meV) corresponding to the peak positions of TA and LA modes in VDOS (see figure 2(c)), respectively; w w D = -; LA P TA P k r and k θ are the DFT-calculated bondstretching and bond-angle-bending force constants (N m −1 ), respectively; a is the lattice constant (Å) at zero temperature. For a comparison, we list w TA P and w LA P of all the investigated groups IV, III-V, IIB-VI and IB-VII ZB semiconductors in table 1. Among all these compounds, we find that (i) CuCl has the lowest w TA P . (ii) For IB-VII semiconductors, as moving along CuCl→CuBr→CuI series, w TA P increases and w w D D = -( ) LA P TA P decreases, and CuI has the smallest Δ and thus lowest w LA P . This result indicates that the LA modes of CuI is relatively easy to be excited compared with other ZB materials at low temperature. Therefore, in CuI the volume contraction induced by the TA modes will be compensated by volume expansion due to the excitation of some LA modes. Furthermore, we find from figure 3 that the magnitude of negative Grüneisen parameter of TA modes in CuI is much smaller than that of CuCl. (iii) As to the IIB-VI, III-V and IV semiconductors, the magnitude of negative a L M increases from ZnS to ZnTe, from GaP to GaSb and from Si to Ge to Sn with decreasing w TA P .

A general trend of NTE against bond-angle-bending force constant
Thus far, we have illustrated the excitations of the TA modes is responsible for the NTE of semiconductors. However, phonon vibrations are not the most fundamental quantities for determining NTE regarding their dependence on the atomic factors. To further reveal the fundamental atomic origin of the NTE, we employ an intuitive but sophisticated method, i.e. the valence force field (VFF) theory to study the volume changes, which are described by a bond-stretching force constant k r and a bond-angle-bending force constant k θ , respectively. The VFF potential is written as [50,51] å å å å å å = -+ +  Here, r 0 is the equilibrium interatomic distance, and ( ) t s r , m is the position vector connecting atom s in unit t to its mth nearest neighbor. Force constant k r characterizes the bond stiffness and k θ is strongly related to the ionicity of a material [50,52,53]. Generally, if a solid has strong ionicity, it has a small k θ because the Coulomb interaction between cation and anion is less sensitive to bond-angle bending. Therefore, from groups IV, III-V, IIB-VI to IB-VII ZB semiconductors, the ionicity increases and thus the bond-angle-bending force constant k θ decreases [53,54]. Figure 4 shows VFF calculated α L as a function of k θ (solid line) for artificial solids with same k r =22.90 N m −1 and bond-length r 0 =2.38 Å but different k θ , in comparison with first-principles DFT calculated a L M (given in table 1) of all investigated ZB and diamond-like semiconductors. One can find that the first-principles calculated a L M are in good agreement with that from the VFF modeling. Since CuCl has the highest ionicity and smallest k θ among all of the ZB and diamond-like semiconductors, it possesses the largest NTE among the examined materials. As moving from CuCl through CuBr to CuI, the increase in k θ combined with the decrease in the phonon gap Δ remarkably amplifies their difference of a L M compared with remaining semiconductors. While, the covalent group IV semiconductors including Si and Ge without ionic characters are among solids having the smallest NTE. The relation that the negative expansion is more marked in more ionic crystals had been noted previously based on empirical observations [22] without clear explanation.

Understanding of NTE from charge density distributions
To elucidate the strong correlation between the ionicity and NTE for a solid, we compare in figure 5 the charge density distributions of groups IB-VII, IIB-VI, III-V, and IV semiconductors, taking CuCl, ZnSe, GaAs, and Si as examples, respectively. Clearly, in IB-VII semiconductors, the charge is mostly located on cation and anion with only few components on the bond. Specifically, Cu loses an electron to become a positively charged cation and Cl accepts this electron to become a negatively charged anion, forming a strong ionic bond. The relative motion between cation and anion perpendicular to the direction of the bond is rather easy in such strong ionic bond,  rendering more negative to the mode Grüneisen parameters of the TA modes. From group IB-VII through IIB-VI, to III-V semiconductors, the ionic character becomes weaker and covalency enhances. Figure 5 shows that the component of change density on the bond significantly gets larger. Specifically, group IV element semiconductors are covalent solids without ionic character and their charge is located mostly on the bonds. The component of the charge on the bond will bound two end atoms of a bond, producing a net torque tending to rotate the bond away from the direction of relative motion. The bond-rotation makes the relative motion transverse to the direction of the bond more difficult, and therefore, the mode Grüneisen parameters of the TA modes are less negative. Furthermore, when the relative motion of two atoms has components both along the bond and transverse to it, the bond-stretching and tension effects contribute simultaneously to the thermal expansion. Hence, compared to covalent solids, ionic compounds more easily excite phonons along transverse vibration whereas suppress longitudinal phonon modes excitations. Therefore, a material with stronger ionicity usually has a larger NTE.
In the covalent limit of group IV elemental semiconductors without ionic characters, as the bond length gets longer from Si through Ge to Sn, the bond-rotation effect and bond strength are both getting weaker. Consequently, the mode Grüneisen parameters of the TA modes become more negative from Si through Ge to Sn. Additionally, the weakness in the bond strength (both k r and k θ become smaller) softens TA modes and makes them excited more easily. Therefore, even in the covalent limit, a L M is different in different solids with its magnitude raises as increasing the bond length of materials (table 1 and figure 4).
Despite the overall agreement, figure 4 also exhibits deviation of the first-principles DFT calculated a L M from the VFF curve, which occurs mostly in groups IB-VII and IIB-VI semiconductors. We should note that, in figure 4, we use the first-principles calculated a L M (solid dots) instead of experimental measured values, since there are no available experimental data for CuBr, CuI, CdS and CdSe and due to potentially inconsistent in different experimental measurements. This deviation is mainly due to underestimation of the first-principles calculation in a L M as we have briefly mentioned above. We can remarkably improve the agreement if we adopt the experimental data such as for CuCl given in table 1 or figure 4 (hollow dot). We attribute this underestimation to the first-principles DFT predicted occupied d orbitals placed too high in energy [55,56]. which leads to the d electrons to be more delocalized than in the real situation. The delocalization in d electrons artificially enhances the covalency of the materials and enlarges the bond-rotation effect. Therefore, we obtain smaller a L M in the first-principles DFT calculations. In reality, the outermost occupied d orbitals are highest in group IB-VII cuprous halides and are getting deeper in IIB-VI, through III-V, to IV semiconductors [55,56]. Consequently, the incorrectness of DFT in the d levels in group IB-VII semiconductors has most remarkably impact on a L M , and then in group IIB-VI semiconductors. Considering a substantial underestimation of a = -´--3.96 10 K

Conclusions
In conclusion, we have investigated systematically and quantitatively the fundamental origin of the anomalous NTE at low temperature in ZB and diamond-like semiconductors. We clarified that the material ionicity controls fundamentally the NTE of those semiconductors. We demonstrated that the NTE becomes significantly larger as the ionicity gets stronger from the groups IV, III-V, IIB-VI to IB-VII ZB semiconductors due to the weakness in the bond-rotation effect, which makes the relative motion between cation and anion transverse to the direction of the bond more feasible and the mode Grüneisen parameters of the TA modes more negative. Therefore, CuCl shows the largest NTE coefficient among all the ZB semiconductors because it has the highest ionicity.