Theoretical Study of Abnormal Thermal Expansion of CuSCN and Effect on Electronic Structure

CuSCN, as a new type of inorganic hole-transporting semiconductor with a wide bandgap (>3.4 eV), is attracting much attention in the fabrication of perovskite solar cells. In this article, by using first-principles density functional theory (DFT) and the quasi-harmonic approximation (QHA) approach, we have studied lattice dynamics and abnormal thermal expansion of the system, including α- and β-CuSCN phases. The influence of the abnormal thermal expansion of the lattice on the electronic structure, especially on the bandgap of the system, was explored and discussed. We found that due to the shearing modes and the three acoustic modes along the direction of the c-axis, the α- and β-CuSCN show a negative thermal expansion (NTE) in the direction of the c-axis. The torsion modes of the Cu–N–C–S atomic chains in the α-CuSCN may lead to an NTE in the directions of the a, b-axes of the α-phase. As a result, our theoretical results demonstrated that the α-CuSCN exhibits an anisotropic bulk NTE. While the β-CuSCN displays a strong uniaxial negative thermal expansion in the direction of the c-axis, in the directions of the a, b-axes, it exhibits positive thermal expansion. Our DFT calculations also predicted that the α-CuSCN has a direct bandgap, which increases slightly with increasing temperature. However, the β-CuSCN has an indirect bandgap at low temperature, which converts to a direct bandgap near the temperature of 375 K due to the strong positive expansion in the ab plane of the phase. Our work revealed the mechanisms of the abnormal thermal expansion of the two phases and a strong coupling between the anisotropic thermal expansion and the electronic structures of the system.


INTRODUCTION
Positive thermal expansion (PTE) is considered to be one of the basic properties of matter in nature. However, some materials have also been found to exhibit abnormal behavior, that is, negative thermal expansion (NTE), within a certain temperature range (Mary et al., 1996;Mittal et al., 2018). The NTE mechanisms could be divided into two categories: phonon-driven NTE or electronic origin NTE. The first is mainly contributed by some low-frequency phonons. This type of mechanism mainly occurs in compounds with framework structures (Gao et al., 2017;Liang et al., 2021), such as oxides (Khosrovani et al., 1997;Wei et al., 2020), fluorides (Greve et al., 2010), MOFs (Wei et al., 2013;Liu et al., 2018), and cyanide/Prussian blue analogs (Goodwin and Kepert, 2005;Gao et al., 2020a). The electronic origin NTE is mainly caused by the temperature dependence of charge transfer (Azuma et al., 2011), magnetic orders, and ferroelectric effects (Pan et al., 2019), such as in perovskite manganese nitride (Takenaka et al., 2008;Sun et al., 2010) and some alloys (Song et al., 2020).
Based on the crystal structure analysis of the NTE materials with framework structures, it is generally believed that the driving force of NTE is the transverse vibration of bridge chain atoms. A typical example is the transverse vibration of O atoms in the Zr-O-W chain of ZrW 2 O 8 (Mary et al., 1996). The NTE of cyanide was first discovered by William et al. in 1997 in Zn(CN) 2 (Collings et al., 2013;Gao et al., 2020a). The NTE comes from the flexibility of the structure with a large number of low-frequency phonon modes, involving the transverse vibrations of the CN groups. In addition, for cyanides, such as AuCN, AgCN, and HT-CuCN, the atomic chains in their crystal structure are arranged in layers, and the thermal expansion shows strong anisotropy in low temperature ranges (Hibble et al., 2010).
Recently, Gao et al. suggested a simple rule based on the concept of average atomic volume (AAV) according to the statistical data of NTE materials, where they proposed a critical AAV of 16 Å 3 for the framework compounds with NTE property, and they believed that the flexibility of the structure is the key factor for the NTE characteristics (Gao et al., 2020c). It provides a simple way to search for and discover new NTE materials. According to the AAV rule, we found that the layer of cyanide-copper thiocyanate (CuSCN), a transparent semiconductor with high hole mobility, has an AVV larger than the value of 16 Å 3 . Because the wide-bandgap semiconductor CuSCN is one of the most ideal inorganic hole transport materials used in perovskite photovoltaic devices and has received extensive attention in recent years (Wijeyasinghe and Anthopoulos, 2015;Arora et al., 2017), it is important and interesting to study the abnormal thermal expansion properties of the material, especially the effect of the abnormal thermal expansion on the electronic structures.
In this work, abnormal thermal expansion and the effect of the abnormal thermal expansion on the electronic structures were investigated for the two phases of the CuSCN, that is, α-CuSCN (Pbca) and β-CuSCN (P63mc), using the combined DFT and QHA approaches. The lattice dynamics calculations indicated that the α-CuSCN shows volume NTE behavior due to the shearing modes between the CuS layers, the acoustic modes, and the torsion modes of the Cu-N-C-S atomic chains, while for the β-CuSCN, the NTE only appears in the direction of the c-axis due to the shearing mode between the CuS layers, the acoustic modes, and the transverse vibration of cyanide (CN), but the lattice parameters in the directions of the a, b-axes show strong PTE and lead to the volume CTE being positive. The calculated electronic property shows that the α-CuSCN displays a direct bandgap semiconductor, while the β-CuSCN has an indirect bandgap at low temperatures. But interestingly, the indirect bandgap of the β-CuSCN converts to a direct bandgap when the temperature increases to 375 K. To our knowledge, the report about the effect of abnormal thermal expansion on the electronic property of CuSCN is yet to be published. The present results revealed the mechanisms of the abnormal thermal expansion of the materials and gave a hint to resolve the puzzle about the direct or indirect bandgap of the β-CuSCN.

COMPUTATION METHODS
In this work, the calculations of the structural and the electronic properties were performed using the first-principles density functional theory (DFT), as implemented in the VASP code (Kresse and Furthmüller, 1996). In the calculations, the projector-augmented wave method (Kresse and Joubert, 1999), with the Perdew-Burke-Ernzerh (PBE) functional (Perdew et al., 1999), was used. The screened hybrid functional (HSE06) (Krukau et al., 2006) was also used to calculate band structures. The first Brillouin zone was integrated with a Monkhost-Pack (Monkhorst and Pack, 1976) k-point mesh (4 × 4 × 3 or 8 × 8 × 3) for the αor β-CuSCN phase, respectively. The cutoff energy of the plane wave basis set is taken to be 500 eV. In the calculations, the atomic positions were relaxed until the force on an atom was smaller than 10 −5 eV/Å and the total energy was converged to 10 −9 eV/atom.
In our calculations, the α-CuSCN primitive cell contains 32 atoms, and the β-CuSCN primitive cell contains eight atoms. The density functional perturbation theory (DFPT) method is used to calculate the force constants of the systems (Baroni et al., 2001), and the phonon spectra are obtained using the PHONOPY code (Parlinski et al., 1997), which is a software package for the calculation of phonons and thermal expansion of crystals under harmonic and quasi-harmonic approximation (QHA). In QHA, the volume dependence of the phonon frequency, that is, part of the anharmonic effect, is considered. Detailed discussion of the QHA can be found in the references (Kuzkin and Krivtsov, 2015). In our previous work, the QHA provided a reasonable description of the dynamic properties of the compound (Liu et al., 2015;Chang et al., 2018;Li et al., 2020).
In order to obtain the coefficient of thermal expansion (CTE), 10 different volumes were taken to calculate the free energies of the 10 volumes around the equilibrium position of systems at a temperature of absolute zero. For each volume which was taken by altering lattice constants by a factor of 0.002, the shape of the unit cell was optimized while keeping the volume of the unit cell fixed, which was implemented by the VASP code. The Vinet equation of state (EOS) (Vinet et al., 1987) was employed to fit these data to get the volumetric CTE of the systems.

Structure and Thermal Expansion Properties
CuSCN is a member of the cyanide family and is a coordination compound formed by thiocyanate ions and Cu metal ions. It has a variety of different structures, including three-dimensional polycrystalline structures and two-dimensional (2D) polycrystalline structures (Tsetseris, 2016a), among which threedimensional polycrystalline structures have received extensive attention as inorganic hole transport materials (Wijeyasinghe and Anthopoulos, 2015;Arora et al., 2017). The three-dimensional polycrystalline structure of CuSCN includes several phases, of which two phases, the α-CuSCN orthogonal phase (Pbca) and the β-CuSCN hexagonal phase (P63mc), are considered here.
Constructed with an atomic chain of the Cu-N≡C-S module, through the twisting of the module and the bonding between the Cu and S atoms, one can achieve the orthogonal α-CuSCN, as shown in Figure 1A. The hexagonal β-CuSCN is obtained through the translation of the module and the connection between the Cu and S atoms, as shown in Figure 1B. For the two phases, each Cu has four coordination atoms, that is, three S and one N. Each S also has four coordination atoms, that is, three Cu and one C.
The optimized lattice parameters and the AAV of the two phases are shown in Table 1. Available experimental and FIGURE 1 | Crystal structure of (A) α-CuSCN and (B) β-CuSCN. The first Brillouin zone and the special k points used in band structure calculations are also presented. Temperature dependence of the cell volume and the CTE of (C) α-CuSCN and (D) β-CuSCN are given.
TABLE 1 | Calculated lattice parameters a, b, and c, bandgaps, and AAV of α-CuSCN and β-CuSCN. Available experimental and theoretical (PBE and HSE06 functionals used) results were also listed for comparison. theoretical results are also listed for comparison. Our theoretical calculations agree well with previous results and show that the β-CuSCN is more stable than the α-CuSCN by about 53 meV per CuSCN chemical unit, which agrees well with the previous theoretical result (Tsetseris, 2016b). The Gibbs free energies of the two phases vs. temperature were calculated and given in Supplementary Figure S1 in the electronic supplementary information (ESI), and they show that the β-CuSCN is more stable than the α-CuSCN in the whole temperature range (0−400 K). The volume thermal expansion behavior of α-CuSCN and β-CuSCN was calculated, as shown in Figure 1C and Figure 1D, respectively. The α-CuSCN shows obvious volumetric NTE characteristics, and its average volumetric coefficient of thermal expansion (CTE) approaches −6.68 × 10 −6 K −1 (0−400 K), while the β-CuSCN shows strong volume PTE behavior (α v 42.09 × 10 −6 K −1 , 0-400 K). It should be noted that though the β-CuSCN has a larger AAV than the α-CuSCN, only the α-CuSCN that was predicted to do so in our calculations exhibits NTE behavior. Seemingly, the criterion of NTE based on the AAV value is not strict, especially when comparing the thermal expansion behavior of two kinds of open-framework structures with different symmetry, because there are other factors, such as electronegativity, symmetry, chemical bond, and so on, to influence thermal expansion of materials. Anyhow, if the AAV of a system is larger than 16 Å 3 , it is very worth it to investigate its thermal expansion behavior and speed up to find new NTE compounds.

Phonon Dispersion and Grüneisen Parameters
In order to investigate the mechanisms of the abnormal thermal expansion of α-CuSCN and β-CuSCN, the phonon dispersion and Grüneisen parameters of different phonon modes were calculated. The phonon spectra and phonon density of states (PDOS) of the two phases are shown in Figures 2A,B, respectively. The phonon spectra of the two phases are roughly divided into four regions: two low-frequency regions (below 300 cm −1 ), one intermediate-frequency region (400-430 cm −1 and 740-760 cm −1 ), and one high-frequency region (2130-2180cm −1 ). It can be seen from the partial PDOS and the eigenvectors of the represented normal modes shown in Figure 3B (α-phase) and Figure 4B (β-phase) that the vibrations in the low-frequency region (below 300 cm −1 ) are mainly contributed from the vibration of Cu, N, and S. In the intermediate-frequency region, the modes are originated from the opposite transverse movement of C and N (400-430 cm −1 ) and the stretching vibration of the S and CN groups along the direction of the Cu-N≡C-S atomic chain (740-760 cm −1 ). In the high-frequency region, it is caused by the stretching vibration of the C≡N triplet bond. These phonon frequency regions can be well explained by the structure and chemical bonds of the system (Liu et al., 2015). According to previous  studies (Dove and Fang, 2016;Mittal et al., 2018;Gao et al., 2020b), we know that phonons in the low-frequency region have a great influence on the thermal expansion characteristics at low temperatures. Different from other cyanides, such as Prussian blue analogs, where the vibrations of metals are weak, the Cu atoms in CuSCN vibrate drastically at low frequencies, and then the vibrations of the S, N, and C atoms are dominated as the frequency increases. The Grüneisen parameters of the two phases were calculated to further understand the influence of phonons in the low-frequency region on NTE, as shown in Figures  2C,D. It can be clearly seen from the figures that the negative values of the Grüneisen parameters of the two phases appear mainly in the frequency range of less than 200 cm −1 . In the α-CuSCN, the number of phonon modes with negative values of Grüneisen parameters is greater than that of those with positive values; thus, it shows the volumetric NTE behavior. In contrast, in the β-CuSCN, it has relatively fewer negative values of Grüneisen parameters in the low-frequency region and a large number of positive values in a wide frequency range, so it displays PTE behavior. The difference in the Grüneisen parameters should be the reason for the obvious difference in the thermal expansion property of the two phases.

Analysis of Phonon Vibration Mode Based on Lattice Dynamics
In order to reveal the mechanism of the abnormal thermal expansion of the two phases derived using the low phonon modes, some representative low phonon modes with negative Grüneisen parameters are presented in Figures 3B, 4B. The lattice constants, as a function of temperature, are calculated and given in Figures 3A, 4A for the two phases, respectively. For the α-CuSCN, eight representative phonon modes (16 cm −1 -416 cm −1 ) with negative Grüneisen parameters are shown in Figure 3B. The vibration mode (−26.1 cm −1 ) with the largest negative Grüneisen parameter (−15.6) comes from the shearing motion between the CuS layers. The three acoustic modes [two transverse acoustic (TA) modes and one longitudinal acoustic (LA) mode] along the direction of the c-axis also have large negative Grüneisen parameters (−14.2, −11.0, and −10.3). These modes contribute NTE in the direction of the c-axis of the system. The torsion vibration modes (46.7 cm −1 and 48.5 cm −1 ) with negative Grüneisen parameters (−7.0 and −7.4) may lead to NTE in the ab plane, which may be the reason for the negative expansion of a ,b below 100 K shown in Figure 3A. The vibration modes for C and N (in a CN group) in phase (143.3 cm −1 ) and opposite phase (412.3 cm −1 ) also have negative Grüneisen parameters (−1.2 and −0.4), which also contribute to the NTE in the direction of the c-axis. For the β-CuSCN, six representative vibration modes (7 cm −1 -416 cm −1 ) with negative Grüneisen parameters are shown in Figure 4B. Similar to the α-phase, the vibration mode (−33.9 cm −1 ) with the largest negative Grüneisen parameter (−17.6) comes from the shearing motion between the CuS layers. The three acoustic modes [two transverse acoustic (TA) modes and one longitudinal acoustic (LA) mode] along the direction of the c-axis also have large negative Grüneisen parameters (−10.4, −11.7, and −10.5). These modes contribute to NTE in the direction of the c-axis of the system just like the case in the α-phase. The lack of the torsion vibration modes may be the reason for the positive expansion of a, b of the β-phase in the temperature range of 0-400 K, as shown in Figure 4A. The vibration modes for C and N (in a CN group) in phase (143.5 cm −1 ) and opposite phase (416.4 cm −1 ) also have negative Grüneisen parameters (−1.0 and −0.5), which also contribute to the NTE in the direction of the c-axis.
In conclusion, due to the shearing and acoustic phonon modes, both the αand β-CuSCN exhibit NTE in the direction of the c-axis. For α-CuSCN, the a, b-axes also show negative thermal expansion in the low-temperature region, which may be due to the contributions from the torsion modes which are not found in the β-CuSCN. Therefore, the strong anisotropic thermal expansion is shown for the β-CuSCN, which is similar to the cyanide AuCN reported previously (Hibble et al., 2010). It should be pointed out that the acoustic modes play an important role in the NTE of the system. A more detailed discussion about the role of acoustic modes in NTE of materials can be found in the recent review of Liang et al. (2021).

Electronic Structure and Effect of Abnormal Thermal Expansion
In this section, we investigated the influence of the abnormal thermal expansion on the electronic structures of the two phases. First, we calculated the energy band structure and electronic partial density of states (EPDOS) of the α-CuSCN and β-CuSCN, as shown in Figures  5A,B. The calculated bandgap of the two phases is given in Table 1, and the previous theoretical and experimental results are also listed for comparison. It can be seen from the energy band structures and EPDOS that the α-CuSCN is a direct bandgap semiconductor with a bandgap of 2.35 eV. The valence band maximum (VBM) at the Γ point is mainly contributed from Cu 3d orbitals, while the conduction band minimum (CBM) at the Γ point is from C and N 2p orbitals. The bandgap of the β-CuSCN is 2.15 eV, which is an indirect bandgap where the VBM at the Γ point is also mainly contributed from Cu 3d orbitals, but the CBM is located at the K point, which is originated from the antibonding molecular state of (CN) −1 formed by 2p orbitals (Jaffe et al., 2010). If a screened hybrid functional (HSE06) was used to calculate the band structure of the β-CuSCN, the obtained bandgap can reach 3.38 eV (see Table 1), but it was still predicted to be an indirect bandgap. Our calculated band structures of the two phases are in good agreement with the previous theoretical results (Smith and Saunders, 1982;Jaffe et al., 2010;Tsetseris, 2016a;Ji et al., 2012). In fact, the previous theoretical calculations predict that the β-CuSCN is an indirect band semiconductor, which seems to disagree with the experimental observation (Ji et al., 2012;Pattanasattayavong et al., 2017). According to our following results about the effect of the abnormal thermal expansion on the band structure of the β-CuSCN, the disagreement may be explained by an indirect-to-direct band transition induced by the strong thermal expansion in the ab plane.
Due to the GGA, the bandgaps of the two phases were underestimated. Compared with the experimental values, the calculated gap values are smaller, but they are consistent with the previous theoretical results (Ji et al., 2012;Pattanasattayavong et al., 2017). If the screened hybrid functional (HSE06) were used, the calculated bandgap would be very close to the experiment, but it was still predicted to be an indirect bandgap.
In order to clearly understand the bond properties between atoms in the two phases, we also calculated and plotted the total charge density in the plane involving the Cu-N≡C-S chain, as shown in Figures 6A,B, for the two phases, respectively. From the calculated total charge density, the bond features can be found, where C and N form a strong covalent bond, S and C form a weak covalent bond, N and Cu form a dative bond, and S and Cu form an ionic bond. These bond features are closely related to the character of different phonon modes. In order to further explain the bond nature, the calculated net charge of atoms was obtained according to the Bader analysis and is given in Supplementary Table S1 in the ESI. From the calculations, it can be seen that for both phases, the net charge of N, C, S, and Cu is −1.2|e|, 0.6|e|, 0.0|e|, and 0.6|e|, indicating that Cu and C lose electrons, while N gains electrons (|e| is charge of an electron). Born effective charges were also calculated by moving a small displacement (0.001 × c Å) of the atoms along the c direction, which are also listed in Supplementary Table S1. For the α-CuSCN, the calculated Born effective charges are Z Cu * 0.12|e|, Z S * −0.73|e|, Z N * −0.32|e|, and Z C * 0.95|e|. For the   β-CuSCN, they are Z Cu * 0.11|e|, Z S * −0.39|e|, Z N * −0.67|e|, and Z C * 0.95|e|. Thermal expansion alters the average distance of atoms in materials, and consequently, the electronic structure of materials is perturbed. Thermal expansion is one of the temperature effects. The electron-phonon coupling interactions due to atom displacement are another important temperature effect on the electronic structure (Zacharias et al., 2020). Here, only the effect of thermal expansion on the electronic structures was considered. According to the calculated lattice constants at different temperatures, the bandgaps at different temperatures have been calculated and are given in Figure 7 for the αand β-CuSCN. The band structures of the β-CuSCN at 0 and 375 K are shown in Figure 8. The band structures of the α-CuSCN at different temperatures are given in Supplementary Figure S2 in the ESI. From the calculated temperature dependence of bandgaps, it can be found that for the α-CuSCN, the gap increases a little (Δ 7 meV) as the temperature increases from 0 to 500 K. Meanwhile, for the β-CuSCN, the gap increases much more (Δ 57 meV) as the temperature increases from 0 to 375 K and then decreases as the temperature further increases. It is interesting to find that for the β-CuSCN, the indirect band structure is converted to be a direct band structure as the temperature increases from 0 to 375 K (Figure 8). Comparing the band structure of the β-CuSCN at low temperature to that at high temperature, we can find the CBM at the K point shifts upward and becomes higher than the minima level of the conduct band at the G point at 375 K. The CBM at the K point is mainly derived from the antibonding states of C and N 2p orbitals (see Figure 5), and because the overlap of the antibonding states between neighboring cells becomes weaker as temperature increases (the distance between the Cu-N-C-S chains increases), the CBM at the K point shifts upward. From the viewpoint of total energy, the total energy of the β-CuSCN increases by about 0.04 eV/unit cell from the structure at 0 K (indirect bandgap) to the structure at 375 K (direct bandgap). This thermal expansion-induced indirect-to-direct bandgap transition may be similar to the indirect-to-direct band transition induced by element substitution, where Se substitution (CuSeCN) leads to the indirect bandgap of CuSCN converting to the direct bandgap of CuSeCN due to the larger lattice constants in the ab plane in CuSeCN (Smith and Saunders, 1982).

CONCLUSION
In this work, abnormal thermal expansion and its effect on the electronic structures of αand β-CuSCN have been studied by combining the first-principles DFT and QHA approaches. Due to the different crystal structures of the two phases, the α-phase exhibits anisotropic volumetric NTE property and the β-phase shows strong anisotropic PTE property where the ab plane (CuS layer) exhibits strong PTE, while there is linear NTE in the direction of the c-axis (perpendicular to the CuS layer). The shearing mode of CuS layers and the acoustic modes play key roles in the anisotropic NTE of the two phases. The torsion modes of the α-CuSCN may be responsible for the NTE of the ab plane at low temperatures. Our results further demonstrated that the α-CuSCN is a direct bandgap semiconductor with a bandgap of 2.35 eV, and the gap becomes a little smaller as the temperature increases from 0 to 500 K. The β-CuSCN was predicted to be an indirect bandgap semiconductor with a bandgap of 2.15 eV at 0 K, where the CBM is at the K point. Due to the strong PTE in the ab plane, the distance between the parallel CN groups increases, which weakens the overlap interactions between the (CN) −1 antibonding states and leads to the states shifting upward, and an indirect-to-direct band transition occurs. We predicted that the system converts to a direct band semiconductor at about 375 K. The present study revealed the mechanisms of the abnormal thermal expansion in the α-CuSCN and the β-CuSCN and shed light on the puzzle regarding the direct or indirect band structures of the β-CuSCN.