First Principles Investigation of Anomalous Pressure-Dependent Thermal Conductivity of Chalcopyrites

The effect of compression on the thermal conductivity of CuGaS2, CuInS2, CuInTe2, and AgInTe2 chalcopyrites (space group I-42d) was studied at 300 K using phonon Boltzmann transport equation (BTE) calculations. The thermal conductivity was evaluated by solving the BTE with harmonic and third-order interatomic force constants. The thermal conductivity of CuGaS2 increases with pressure, which is a common behavior. Striking differences occur for the other three compounds. CuInTe2 and AgInTe2 exhibit a drop in the thermal conductivity upon increasing pressure, which is anomalous. AgInTe2 reaches a very low thermal conductivity of 0.2 W·m−1·K−1 at 2.6 GPa, being beneficial for many energy devices, such as thermoelectrics. CuInS2 is an intermediate case. Based on the phonon dispersion data, the phonon frequencies of the acoustic modes for CuInTe2 and AgInTe2 decrease with increasing pressure, thereby driving the anomaly, while there is no significant pressure effect for CuGaS2. This leads to the negative Grüneisen parameter for CuInTe2 and AgInTe2, a decreased phonon relaxation time, and a decreased thermal conductivity. This softening of the acoustic modes upon compression is suggested to be due to a rotational motion of the chalcopyrite building blocks rather than a compressive oscillation. The negative Grüneisen parameters and the anomalous phonon behavior yield a negative thermal expansion coefficient at lower temperatures, based on the Grüneisen vibrational theory.


Introduction
Chalcopyrite compounds (A I B III C 2 VI , A I = IB elements (Cu, Ag), B III = IIIA elements (Al, Ga, In), C VI = VIA elements (S, Se, Te), space group I-42d, as shown in Figure 1) are well-known semiconductors with a band gap in the range of 0.1 to 1 eV [1][2][3][4]. The calculated band gap for the selected systems is 1.085, 0.364, 0.469, and 0.967 eV for CuGaS 2 , CuInS 2 , AgInS 2 , and AgInTe 2 , respectively [1][2][3][4][5], being consistent with common density functional theory deviations [6]. Their structure can be derived from zincblende (ZnS) by alternating the A I and B III constituents at the Zn site [7]. In the zincblende structure, a Zn atom is located in a center of a tetrahedron span by S, which is equivalent to an A I -or B III -based tetrahedron span by C VI [7]. Since the Grüneisen parameter and volumetric thermal expansion coefficient can be related to κ (see below for more details), it appears that pressure effects on κ are considerable. This is consistent with an experimental study reporting a decrease in κ by 30% for CuInTe2 under pressure up to 2.3 GPa [24]. Generally, κ should increase under compression [22]; however, the compounds in the current study exhibit a mixture of monotonic decrease and non-monotonic dependence under pressure. The latter phenomenon has been intensively studied recently [18,19,21]. This implies that the behavior of CuInTe2 and related compounds is of great interest due to non-monotonicity. The behavior of CuInTe2 is anomalous. Two possible mechanisms have been proposed based on experiments: (i) anharmonic behavior of lattice vibrations [24] and (ii) structural modifications under high pressure (e.g., stacking faults) [25]. The underlying physics of the κ reduction under compression of CuInTe2, and possibly other A I B III C2 VI compounds, is not fully understood.
In this work, we devise a strategy to identify the physical origin of the anomalous behavior of κ of A I B III C2 VI compounds under compression. Using phonon calculations, the atomic-level understanding of this anomaly is obtained by analyzing the vibrational modes and correlating these to the macroscopic observables, such as κ. To systematically explore the pressure effect on κ, CuInTe2 is taken as a reference, and the influence of mass, being decisive for lattice vibrations, is considered by replacing Cu with Ag, In with Ga, and Te with S within this isostructural and isoelectronic A I B III C2 VI system. Hence, CuGaS2, CuInS2, CuInTe2, and AgInTe2 chalcopyrites are explored.

Methods
In order to maximize ZT, σ and κ should be maximized and minimized, respectively. An increase in σ directly affects the total κ value, since = + , where κe is the electronic thermal conductivity (charge carriers also conduct heat) and κph is the lattice thermal conductivity. Hence, minimizing κph is the major route to minimize the total κ value. Furthermore, κph is likely the largest contribution for chalcopyrites since they are semiconductors. It can be obtained as follows: = These A I B III C 2 VI compounds have been extensively studied for photovoltaic and thermoelectric applications [8][9][10][11][12][13]. The efficiency of thermoelectric devices depends on the thermoelectric figure of merit ZT = TS 2 σ/κ, where T is the absolute temperature, S designates the Seebeck coefficient, σ is the electric conductivity, and κ stands for the thermal conductivity [14]. Other chalcopyrite compounds such as CuFeS 2 are currently under study for their mechanical, electronic, and thermodynamic properties, being promising for thermoelectric applications [15,16]. It is well established that κ decreases at elevated temperatures, which is beneficial for thermoelectric applications [17]. Different studies showed that pressure has an anomalous effect on A II B IV zincblende compounds [18,19]. Using the Slack model, Gui et al. [2] showed that ZT of CuInC 2 VI (C VI = S, Se, and Te) uniformly increases at elevated temperatures up to 850 K. While the effect of pressure on κ and its relation to other thermal properties has been explored, the thermal expansion coefficient [20][21][22] has not been thoroughly studied. Furthermore, atomic vibrations driving anomalous thermal behavior of CuInC 2 VI compounds are not known and cannot simply be deduced from other systems. Using the quasi-harmonic Debye model, Sharma et al. have evaluated electronic, thermal, and mechanical properties of AgInC 2 VI (C VI = S, Se, and Te) under pressure and reported a noticeable reduction in the Grüneisen parameter and volumetric thermal expansion coefficient, as well as the bulk modulus [23].
Since the Grüneisen parameter and volumetric thermal expansion coefficient can be related to κ (see below for more details), it appears that pressure effects on κ are considerable. This is consistent with an experimental study reporting a decrease in κ by 30% for CuInTe 2 under pressure up to 2.3 GPa [24]. Generally, κ should increase under compression [22]; however, the compounds in the current study exhibit a mixture of monotonic decrease and non-monotonic dependence under pressure. The latter phenomenon has been intensively studied recently [18,19,21]. This implies that the behavior of CuInTe 2 and related compounds is of great interest due to non-monotonicity. The behavior of CuInTe 2 is anomalous. Two possible mechanisms have been proposed based on experiments: (i) anharmonic behavior of lattice vibrations [24] and (ii) structural modifications under high pressure (e.g., stacking faults) [25]. The underlying physics of the κ reduction under compression of CuInTe 2 , and possibly other A I B III C 2 VI compounds, is not fully understood.
In this work, we devise a strategy to identify the physical origin of the anomalous behavior of κ of A I B III C 2 VI compounds under compression. Using phonon calculations, the atomic-level understanding of this anomaly is obtained by analyzing the vibrational modes and correlating these to the macroscopic observables, such as κ. To systematically explore the pressure effect on κ, CuInTe 2 is taken as a reference, and the influence of mass, being decisive for lattice vibrations, is considered by replacing Cu with Ag, In with Ga, and Te with S within this isostructural and isoelectronic A I B III C 2 VI system.

Methods
In order to maximize ZT, σ and κ should be maximized and minimized, respectively. An increase in σ directly affects the total κ value, since κ = κ e + κ ph , where κ e is the electronic thermal conductivity (charge carriers also conduct heat) and κ ph is the lattice thermal conductivity. Hence, minimizing κ ph is the major route to minimize the total κ value. Furthermore, κ ph is likely the largest contribution for chalcopyrites since they are semiconductors. It can be obtained as follows: where v g is the group velocity of phonons, c v designates the heat capacity, and τ is the phonon relaxation time [26]. The former two values were calculated herein from the phonon dispersion curves using the Phonopy package [27], while τ was obtained by solving the Boltzmann transport equation, as implemented in the ShengBTE package [28]. The isotropic approximation was applied due to isotropic pressure dependence evaluated in the current study. Harmonic and third-order interatomic force constants within three coordination shells were used as input in both packages. All interatomic force constants were generated using the Vienna Ab-initio Simulation Package (VASP) [29][30][31][32][33]. The exchange-correlation functionals were treated within the local density approximation [34], including phonon calculations. The all-electron projector augmented wave method [35] was utilized to evaluate electronic wave functions with a plane wave cutoff of 800 eV and the total energy convergence of 10 −7 eV. No configurations were spin polarized. All structures (internal free parameters) and unit cell sizes were optimized within a force convergence condition of 10 −6 eV·Å −1 . A 4 × 4 × 2 k-mesh Monkhorst-Pack [36] was used to sample the Brillouin zone (BZ) of the 2 × 2 × 1 supercell (64 atoms) constructed from the conventional A I B III C 2 VI unit cell (16 atoms, six coordination shells). Convergence tests were conducted for the k-mesh (from 2 × 2 × 1 to 6 × 6 × 3). At the chosen conditions (4 × 4 × 2), phonon band structure was converged and stable under pressure up to 9 GPa. Since a 10 × 10 × 10 q-mesh leads to a fluctuation in κ ph on the order of 10 −1 W·m −1 ·K −1 (348 atomic displacements per chalcopyrite configuration accumulating Hellmann-Feynman forces for the phonon calculations), such changes are acceptable for the calculated κ ph values in the range of several W·m −1 ·K −1 . As stated in the introduction, all A I B III C 2 VI chalcopyrite compounds exhibit a band gap of >0.4 eV (see Table 1) so that only phonons are considered for the evaluation of the transport properties. The pressure dependence was modelled in the form of isotropic compressive strains in steps of 3% of the equilibrium unit cell volume (a minimum of 5 strains were considered per chalcopyrite configuration). The bulk moduli were calculated using the Rose-Vinet equation of state to evaluate the pressure associated with each strain [37]. The calculated bulk moduli for CuGaS 2 , CuInS 2 , CuInTe 2 , and AgInTe 2 were 76.0, 64.6, 41.3, and 41.2 GPa, respectively, which is consistent with the literature [38][39][40][41][42] (see Table 1). The obtained lattice parameters were a = 5.387 Å and c/a = 1.981 for CuGaS 2 , a = 5.597 Å and c/a = 2.015 for CuInS 2 , a = 6.303 Å and c/a = 2.007 for CuInTe 2 , and a = 6.582 Å and c/a = 1.978 for AgInTe 2 . These lattice constants deviate max. 2.7% from the experimental data [43,44] (see Table 1), which is acceptable for the employed exchange-correlation functionals [6]. Finally, the quasi-harmonic approximation [45] was utilized to calculate the volumetric thermal expansion coefficient due to its relationship with κ ph , as detailed below.

Lattice Thermal Conductivity
In Figure 2, the lattice thermal conductivity of all A I B III C 2 VI compounds studied in this work is plotted against an increasing pressure. First, the ambient pressure conditions (0 GPa) are discussed. CuGaS 2 exhibits a relatively large κ ph value of 8.2 W·m −1 ·K −1 , which deviates 12% from the measured one (9.3 W·m −1 ·K −1 ) [46]. In the case of CuInS 2 , κ ph reaches 4.6 W·m −1 ·K −1 . The κ ph value of AgInTe 2 is 2.9 W·m −1 ·K −1 , which is only 7% offset from the experimental value of 2.7 W·m −1 ·K −1 [47]. CuInTe 2 attains 7.6 W·m −1 ·K −1 , which is an 18% deviation compared to the experimentally obtained value of 6.2 W·m −1 ·K −1 [48]. Hence, these differences are acceptable based on the κ ph deviations of other theoretical studies from measurements [46,49,50]. See Table 1 for comparison. Furthermore, based on a comparison of predicted and experimental κ ph values for diamond, SiC, GaN, Si, GaAs, InSb, SrTiO 3 , and PbTe, theoretical data at elevated temperatures commonly overestimate the measured values [51], so that the calculated low κ ph data in this study should be even lower under ordinary experimental conditions (e.g., presence of defects). It is clear that trends are properly captured within the methodology used in the current study.
Materials 2019, 12, x FOR PEER REVIEW 4 of 10 attains 7.6 W·m −1 ·K −1 , which is an 18% deviation compared to the experimentally obtained value of 6.2 W·m −1 ·K −1 [48]. Hence, these differences are acceptable based on the κph deviations of other theoretical studies from measurements [46,49,50]. See Table 1 for comparison. Furthermore, based on a comparison of predicted and experimental κph values for diamond, SiC, GaN, Si, GaAs, InSb, SrTiO3, and PbTe, theoretical data at elevated temperatures commonly overestimate the measured values [51], so that the calculated low κph data in this study should be even lower under ordinary experimental conditions (e.g., presence of defects). It is clear that trends are properly captured within the methodology used in the current study. Striking differences were obtained between the κph behavior of CuGaS2, CuInS2, CuInTe2, and AgInTe2 under compression (see Figure 2). The pressure range explored herein can be reached in samples synthesized by non-equilibrium vapor phase condensation processes [52]. The κph value of CuGaS2 always increases with pressure up to 9.5 GPa. This is a common behavior for most compounds [22]. A drastically different dependence is obtained as soon as the heavier In is considered instead of Ga (B III in A I B III C2 VI ). Under pressure up to 6 GPa, the κph of CuInTe2 decreases from 7.6 to 4.1 W·m −1 ·K −1 , which is anomalous. This is consistent with the experimental data [24,25], implying that important physics is captured within the methodology employed herein and structural modulations are not indispensable to drive the anomaly. By exchanging Te with lighter S (C VI in A I B III C2 VI ) and hence forming CuInS2, κph increases up to 2 GPa, which is again a common behavior and equivalent to that of CuGaS2. Upon a further pressure increase, κph begins to decrease and reaches a slightly lower value at 8 GPa than that at 0 GPa. To account for the effect of the transition metal constituent (A I in A I B III C2 VI ), Cu in CuInTe2 is exchanged with the heavier Ag. AgInTe2 exhibits a significantly lower κph value and a steeper decrease in κph under pressure, reaching 0.2 W·m −1 ·K −1 at 2.6 GPa. This is a very low value for κph and is comparable to that of some polymers, such as polytetrafluoroethylene [53]. Such low κph values under pressure should lead to an enhanced thermoelectric performance, according to several studies [54][55][56]. Striking differences were obtained between the κ ph behavior of CuGaS 2 , CuInS 2 , CuInTe 2 , and AgInTe 2 under compression (see Figure 2). The pressure range explored herein can be reached in samples synthesized by non-equilibrium vapor phase condensation processes [52]. The κ ph value of CuGaS 2 always increases with pressure up to 9.5 GPa. This is a common behavior for most compounds [22]. A drastically different dependence is obtained as soon as the heavier In is considered instead of Ga (B III in A I B III C 2 VI ). Under pressure up to 6 GPa, the κ ph of CuInTe 2 decreases from 7.6 to 4.1 W·m −1 ·K −1 , which is anomalous. This is consistent with the experimental data [24,25], implying that important physics is captured within the methodology employed herein and structural modulations are not indispensable to drive the anomaly. By exchanging Te with lighter S (C VI in A I B III C 2 VI ) and

Acoustic Phonon Dispersion
hence forming CuInS 2 , κ ph increases up to 2 GPa, which is again a common behavior and equivalent to that of CuGaS 2 . Upon a further pressure increase, κ ph begins to decrease and reaches a slightly lower value at 8 GPa than that at 0 GPa. To account for the effect of the transition metal constituent (A I in A I B III C 2 VI ), Cu in CuInTe 2 is exchanged with the heavier Ag. AgInTe 2 exhibits a significantly lower κ ph value and a steeper decrease in κ ph under pressure, reaching 0.2 W·m −1 ·K −1 at 2.6 GPa. This is a very low value for κ ph and is comparable to that of some polymers, such as polytetrafluoroethylene [53].
Such low κ ph values under pressure should lead to an enhanced thermoelectric performance, according to several studies [54][55][56].

Acoustic Phonon Dispersion
Since κ ph and all the corresponding factors, v g , c v , and τ, are governed by phonons, phonon dispersion curves are further discussed to explain the anomalous behavior. The acoustic phonon modes at 300 K contribute 85%, 70%, 60%, and 80% of κ ph for CuGaS 2 , CuInS 2 , CuInTe 2 , and AgInTe 2 , respectively, so that non-locality of the exchange-correlation functional is of less significance since it would mainly contribute to optical phonon frequencies [57]. Therefore, in this study, the behavior of acoustic phonon modes and the effect of pressure thereon are considered in detail. Since AgInTe 2 and CuInTe 2 exhibit the same κ ph behavior, whereby the former undergoes a more drastic change, AgInTe 2 is taken as a representative. Figure 3   The atomic vibrations were analyzed from the phonon dispersion curves using a tool developed by Miranda et al. [58]. For the phonon modes with decreasing frequency upon compression, e.g., LA of AgInTe2 starting at (0, 0, 0.3) on the Г-X and X-P path (30% of the path length), as shown in Figure  3c, the vibrations are unconventional. The metal-centered tetrahedra (A I in A I B III C2 VI ) oscillate in the manner that the angles between neighboring units are changing, while the bond lengths are fixed. However, the conventional thermal vibrations occur in the form of a compression wave, where the bonds are mostly stretched, as in the case of CuGaS2. The anomalous behavior is due to the circular motion of the metal atoms (A I = Cu or Ag) around the equilibrium position. The anomalous chalcopyrites tend to keep the bond length fixed and the excitations appear in the form of bond bending. Such behavior is consistent with the tension effect introduced by Dove et al. [59], showing that in the case where the energy required for stretching is too high, a finite transverse displacement in the centers of polyhedra occurs, resulting in rotation.

Phonon Relaxation Time
A decrease in the phonon frequency upon compression (softening of the acoustic phonon modes) implies the negative Grüneisen parameter (γ), which can be related to τ (1/τ = γ 2 within the Debye-Callaway model [60]) one of the three physical variables describing κph. Since vg 2 decreases uniformly with pressure for all A I B III C2 VI compounds explored in this work and cv is affected by only 0.1%, it appears that τ is the major constituent responsible for the anomaly. Hence, τ is further explored for CuGaS2, CuInS2, and AgInTe2 as a function of pressure. The corresponding τ data for the acoustic phonon modes are shown in Figure 4. For CuGaS2 (Figure 4a), τ increases up to 300% with pressure. This implies that the absolute value of γ decreases. For CuInS2 (Figure 4b), τ increases up to 2 GPa in the case of TA and further to 4 GPa in the case of LA and TAs and then decreases, which is supported by the fact that γ changes the sign for a number of low-frequency phonons as appearing The atomic vibrations were analyzed from the phonon dispersion curves using a tool developed by Miranda et al. [58]. For the phonon modes with decreasing frequency upon compression, e.g., LA of AgInTe 2 starting at (0, 0, 0.3) on the Γ-X and X-P path (30% of the path length), as shown in Figure 3c [59], showing that in the case where the energy required for stretching is too high, a finite transverse displacement in the centers of polyhedra occurs, resulting in rotation.

Phonon Relaxation Time
A decrease in the phonon frequency upon compression (softening of the acoustic phonon modes) implies the negative Grüneisen parameter (γ), which can be related to τ (1/τ = γ 2 within the Debye-Callaway model [60]) one of the three physical variables describing κ ph . Since v g 2 decreases uniformly with pressure for all A I B III C 2 VI compounds explored in this work and c v is affected by only 0.1%, it appears that τ is the major constituent responsible for the anomaly. Hence, τ is further explored for CuGaS 2 , CuInS 2 , and AgInTe 2 as a function of pressure. The corresponding τ data for the acoustic phonon modes are shown in Figure 4. For CuGaS 2 (Figure 4a), τ increases up to 300% with pressure. This implies that the absolute value of γ decreases. For CuInS 2 (Figure 4b), τ increases up to 2 GPa in the case of TA and further to 4 GPa in the case of LA and TAs and then decreases, which is supported by the fact that γ changes the sign for a number of low-frequency phonons as appearing in the dispersion relation in Figure 3b. For AgInTe 2 (Figure 4c), τ decreases in the whole pressure range and hence γ 2 is increasing, while the values are negative (softening of the acoustic phonon modes). Therefore, an increasing τ with pressure gives rise to an increasing κ ph , as in the case of CuGaS 2 . Anomalous κ ph of AgInTe 2 exhibits a small τ value and negative γ. The notion of negative γ values and softening of the acoustic phonon modes is consistent with the literature on AgGaS 2 [61].
In the present work, an important step is made towards a relationship with κ ph .

Thermal Expansion Coeficient
According to the Grüneisen vibrational theory of thermal expansion [59], negative γ yields a negative volumetric thermal expansion coefficient (α) since = ⁄ , where V is the volume and B stands for the bulk modulus. As α can be obtained independently from γ (quasi-harmonic approximation was used herein), probing α is important. Furthermore, this may help other experimentalists, besides those focusing on thermoelectric devices, to critically appraise the results obtained in this study. Therefore, in Figure 5 the α value is plotted at different temperatures. For CuGaS2, having an increasing κph value under pressure, is always positive. CuInS2 possesses a slightly negative α value. However, both CuInTe2 and AgInTe2 exhibit negative α at low temperatures. The pressure dependence on α is not shown here since the major effect is an offset of the negative α region to higher temperatures, conserving the trends between these chalcopyrites. The anomalous behavior of CuInTe2 and AgInTe2 is thus driven by softening of the acoustic phonon modes under compression, leading to negative α and γ.

Thermal Expansion Coeficient
According to the Grüneisen vibrational theory of thermal expansion [59], negative γ yields a negative volumetric thermal expansion coefficient (α) since γ = αV/Bc v , where V is the volume and B stands for the bulk modulus. As α can be obtained independently from γ (quasi-harmonic approximation was used herein), probing α is important. Furthermore, this may help other experimentalists, besides those focusing on thermoelectric devices, to critically appraise the results obtained in this study. Therefore, in Figure 5 the α value is plotted at different temperatures. For CuGaS 2 , having an increasing κ ph value under pressure, α is always positive. CuInS 2 possesses a slightly negative α value. However, both CuInTe 2 and AgInTe 2 exhibit negative α at low temperatures. The pressure dependence on α is not shown here since the major effect is an offset of the negative α region to higher temperatures, conserving the trends between these chalcopyrites. The anomalous behavior of CuInTe 2 and AgInTe 2 is thus driven by softening of the acoustic phonon modes under compression, leading to negative α and γ.
CuGaS2, having an increasing κph value under pressure, is always positive. CuInS2 possesses a slightly negative α value. However, both CuInTe2 and AgInTe2 exhibit negative α at low temperatures. The pressure dependence on α is not shown here since the major effect is an offset of the negative α region to higher temperatures, conserving the trends between these chalcopyrites. The anomalous behavior of CuInTe2 and AgInTe2 is thus driven by softening of the acoustic phonon modes under compression, leading to negative α and γ.

Conclusions
The thermal conductivity of CuGaS2, CuInS2, CuInTe2, and AgInTe2 A I B III C2 VI chalcopyrites exhibits a different behavior under pressure, ranging from increasing (CuGaS2), as for most compounds, alternating (increasing and decreasing for CuInS2), and to decreasing, as in the case for CuInTe2 and AgInTe2, which is anomalous. This can be understood based on the phonon dispersion curves. Softening of the acoustic phonon modes occurs for these anomalous chalcopyrites. This leads to the negative Grüneisen parameter and negative volumetric thermal expansion coefficient. The

Conclusions
The thermal conductivity of CuGaS 2 , CuInS 2 , CuInTe 2 , and AgInTe 2 A I B III C 2 VI chalcopyrites exhibits a different behavior under pressure, ranging from increasing (CuGaS 2 ), as for most compounds, alternating (increasing and decreasing for CuInS 2 ), and to decreasing, as in the case for CuInTe 2 and AgInTe 2 , which is anomalous. This can be understood based on the phonon dispersion curves. Softening of the acoustic phonon modes occurs for these anomalous chalcopyrites. This leads to the negative Grüneisen parameter and negative volumetric thermal expansion coefficient. The decrease in phonon frequency upon compression is suggested to be due to the phonon oscillations in the form of a rotational motion rather than compressive waves. The physical origin of the anomalous thermal conductivity is thus identified in this work in terms of higher-order phonon-phonon interactions, and AgInTe 2 with a very low thermal conductivity of 0.2 W·m −1 ·K −1 at 2.6 GPa is proposed to be a promising thermoelectric compound.