The Temperature-Dependent Thermal Conductivity of C-and O-Doped Si 3 N 4 : First-Principles Calculations

: Silicon nitride (Si 3 N 4 ) possesses excellent mechanical properties and high thermal conductivity, which is an important feature in many applications. However, achieving the theoretically high thermal conductivity of Si 3 N 4 in practice is challenging. In this study, we adopted a first-principles calculation method to assess the effects of doping β -Si 3 N 4 and γ -Si 3 N 4 with carbon and oxygen atoms. Applying geometric structure optimization combined with calculation of the electronic phonon properties generated a stable doped structure. The results revealed that carbon and oxygen doping have little effect on the Si 3 N 4 unit cell size, but that oxygen doping increases the unit cell volume. Energy band structure and state density calculation results showed that carbon doping reduces the nitride band gap width, whereas oxygen doping results in an n-type Si 3 N 4 semiconductor. The findings from this study are significant in establishing a basis for targeted increase of the thermal conductivity of Si 3 N 4 .


Introduction
Silicon nitride (Si 3 N 4 ) exhibits excellent mechanical properties, high chemical stability, and superior wear and thermal shock resistance [1,2].Furthermore, Si 3 N 4 also exhibits high thermal conductivity.Its intrinsic thermal conductivity has been reported to be about 200-320 W/m•K at room temperature [3].Currently, silicon nitride ceramics are mainly employed as insulating substrates for high-power electronic devices [4].However, the high thermal conductivity of as-received Si 3 N 4 is not easy to attain in reality.Sintered polycrystalline Si3N4 ceramics have much lower thermal conductivity than the intrinsic thermal conductivity of the compound.It is believed that residual pores, as well as intergranular secondary phase and lattice imperfections including impurity atoms, vacancies, dislocations, and stacking faults, could affect thermal conductivity [5][6][7].
Because of this, many studies have attempted to improve the thermal conductivity of Si 3 N 4 .Horng-Hwa Lu and Jow-Lay Huang [8] systematically explored the effects of raw powders, especially phase composition, on the final microstructure of sintered Si 3 N 4 .And then, Hiroshi Yokota and Masahiro Ibukiyama studied the effects of the raw particle size distribution of Si 3 N 4 powder using Yb 2 O 3 and ZrO 2 as sintering additives [9].They found that it is possible for Si 3 N 4 ceramics to achieve high thermal conductivity and high strength simultaneously by optimizing the particle size.In another study, Y Zhou et al. proposed using sintered reaction-bonded silicon nitride (SRBSN).The new technological method not only had a lower cost but also a higher possibility of attaining high thermal conductivity [10].
Notably, most of these works focused on microstructural features, including the relative density, the grain boundary, and the intergranular secondary phase, but few works have addressed impurity atoms.Oxygen and carbon, which can enter the Si 3 N 4 lattice to Crystals 2024, 14, 549 2 of 9 form a solid solution phase during fabrication, are widely recognized as being the most detrimental impurities to the thermal conductivity of Si 3 N 4 ceramics [11,12].It is said that the impurity atoms can form silicon vacancies that scatter phonons [13,14].Some studies have demonstrated that the elimination of lattice oxygen is crucial for improving the thermal conductivity of Si 3 N 4 ceramics by adding oxide additives [15].
However, with regard to commercial applications, various kinds of silicon nitride powders with different high oxygen contents should be considered as the starting material.It is costly to optimize additivities for different raw powders or choose the proper raw powder.Therefore, a deep understanding of the effect of impurity atoms on thermal conductivity will be helpful for real-world industrial application.Huimin Xiang et al. [16] showed that first-principles theory and its related numerical methods are effective tools for predicting the thermal conductivity of Si 3 N 4 with different crystal lattices.Based on these results, the effects of impurities on the thermal conductivity of different Si 3 N 4 lattices were studied using a similar method.

Theoretical Model
The theoretical analysis considered two crystalline forms of Si 3 N 4 , i.e., β-Si 3 N 4 and γ-Si 3 N 4 .The space group for hexagonal β-Si 3 N 4 is P63/m, as shown in Figure 1.This structure is composed of slightly twisted corner-congruent silicon-nitrogen tetrahedra, forming a hexagonal ring layer.Each β-Si have addressed impurity atoms.Oxygen and carbon, which can enter the Si3N4 lattice to form a solid solution phase during fabrication, are widely recognized as being the most detrimental impurities to the thermal conductivity of Si3N4 ceramics [11,12].It is said that the impurity atoms can form silicon vacancies that scatter phonons [13,14].Some studies have demonstrated that the elimination of lattice oxygen is crucial for improving the thermal conductivity of Si3N4 ceramics by adding oxide additives [15].
However, with regard to commercial applications, various kinds of silicon nitride powders with different high oxygen contents should be considered as the starting material.It is costly to optimize additivities for different raw powders or choose the proper raw powder.Therefore, a deep understanding of the effect of impurity atoms on thermal conductivity will be helpful for real-world industrial application.Huimin Xiang et al. [16] showed that first-principles theory and its related numerical methods are effective tools for predicting the thermal conductivity of Si3N4 with different crystal lattices.Based on these results, the effects of impurities on the thermal conductivity of different Si3N4 lattices were studied using a similar method.

Theoretical Model
The theoretical analysis considered two crystalline forms of Si3N4, i.e., β-Si3N4 and γ-Si3N4.The space group for hexagonal β-Si3N4 is P63/m, as shown in Figure 1.This structure is composed of slightly twisted corner-congruent silicon-nitrogen tetrahedra, forming a hexagonal ring layer.Each β-Si3N4 unit cell contains 14 atoms, and the unit cell parameters are a = b = 7.607 Å, c = 2.911 Å, α = β = 90°, and γ = 120°.The calculations performed in this study used a 2 × 2 × 2 supercell model with a total of 112 atoms, obtained by extending the β-Si3N4 primitive cell by one unit along the basis vector x, y, and z directions.Carbon doping replaces Si atoms in the supercell, whereas oxygen doping replaces N atoms.

Calculations
Two density functional theory (DFT)-based codes, the Vienna Ab initio Simulation Package (VASP) [3,4] and the orthogonalized linear combination of atomic orbitals (OL-CAO) code [5,6], were used in this study.The optimized VASP structures were used as inputs to calculate the electronic structure and interatomic bonding properties using an in-house OLCAO package.The OLCAO package uses the local density approximation (LDA) as the exchange and correlation potential.Atomic orbitals expanded in Gaussiantype orbitals (GTOs) are used for the basis set expansion of the wave functions.The use of localized atomic orbitals in the basis expansion contrasts with the plane-wave expansion used in VASP.
In terms of geometric structure optimization, a generalized gradient approximation (GGA) function was used to approximate the exchange correlation potential [7].Taking into account the subsequent phonon calculation, the two Si3N4 crystal models, β-Si3N4 and γ-Si3N4, were continuously optimized until the energy difference and force difference were less than 1 × 10 −8 eV/atom and 0.001 eV/Å, respectively, and the cutoff energy

Calculations
Two density functional theory (DFT)-based codes, the Vienna Ab initio Simulation Package (VASP) [3,4] and the orthogonalized linear combination of atomic orbitals (OLCAO) code [5,6], were used in this study.The optimized VASP structures were used as inputs to calculate the electronic structure and interatomic bonding properties using an in-house OLCAO package.The OLCAO package uses the local density approximation (LDA) as the exchange and correlation potential.Atomic orbitals expanded in Gaussian-type orbitals (GTOs) are used for the basis set expansion of the wave functions.The use of localized atomic orbitals in the basis expansion contrasts with the plane-wave expansion used in VASP.
In terms of geometric structure optimization, a generalized gradient approximation (GGA) function was used to approximate the exchange correlation potential [7].Taking into account the subsequent phonon calculation, the two Si 3 N 4 crystal models, β-Si 3 N 4 and γ-Si 3 N 4 , were continuously optimized until the energy difference and force difference were less than 1 × 10 −8 eV/atom and 0.001 eV/Å, respectively, and the cutoff energy was set at 500 eV.In addition, 2 × 2 × 4 k mesh point sampling was used for both the β-Si 3 N 4 and γ-Si 3 N 4 systems.The atomic pseudo-potentials described using the projection enhanced plane wave (PAW) method [8] are Si-3s2 p2, N-2s2 p3, O-2s2 2p4, and C-2s2 2p2.The bond order (BO) values are given by Equation (1), also called the overlap population ραβ between any pair of atoms (α, β).
Moreover, Phonopy software (v 2.24.1), based on quasi-harmonic approximation, was used to calculate the second-order force constant and phonon spectrum, applying the density functional perturbation method [9].As phonon thermal conduction is an anharmonic process, the third-order interatomic interaction force constant was considered in this study.

Results and Discussion
The geometry of the two crystal forms of Si 3 N 4 and the configurations of doped C and O atoms were optimized to obtain stable β-Si 3 N 4 , C@β-Si 3 N 4 , O@β-Si 3 N 4 , γ-Si 3 N 4 , C@γ-Si 3 N 4 , and O@γ-Si 3 N 4 ; subsequent calculations were based on the fully relaxed structure.The optimized lattice geometric parameters for the six systems after cell expansion are presented in Table 1.The incorporation of C or O atoms had a minimal effect on the unit cell size.The lattice constant for C@β-Si 3 N 4 was slightly reduced, while the lattice constant of the two Si 3 N 4 crystal forms was slightly increased with the incorporation of O atoms, and the overall effect of Al-N codoping increased.The ionic radius of O is 1.2 Å, which is significantly greater than the ionic radius of N (0.13 Å).The replacement of N by O was accompanied by an increase in the unit cell volume.

Band Structure and Density of States
The application of density functional theory (DFT) provides a calculation of the ground state configuration, i.e., the structure at 0 K, where the GGA approximation underestimates the interaction between excited electrons.This results in an underestimation of the band gap.Following optimization using VASP software (6.3.2), the OLCAO method was used to calculate the energy band and the structural density of states.
The energy band structures of the six systems calculated along the high symmetry point of the Brillouin zone are presented in Figures 2 and 3.As can be seen from Figure 2a, the direct band gap width of intrinsic β-Si 3 N 4 is ca.4.47 eV, which is close to the value (4.33 eV) calculated by Lu et al. [10].Doping with C (Figure 2b) reduced the band gap of β-Si 3 N 4 , where the indirect band gap of 3.13 eV was also close to the value (2.90 eV) reported by Lu et al. [8].The incorporation of C resulted in electrons occupying an increased number of energy levels.When compared with the intrinsic band structure of β-Si 3 N 4 , the Fermi level entered the conduction band after O doping (Figure 2c), making the doped nitride an n-type semiconductor material.The number of energy levels increased, and the bandgap width decreased (E g = 3.20 eV).In the case of γ-Si 3 N 4 , the direct band gap was ca.3.35 eV.Carbon doping (Figure 3b) slightly increased the band gap to 3.47 eV and significantly increased the number of energy levels occupied by electrons.Following the incorporation of O atoms (Figure 3c), the Fermi level entered the conduction band, γ-Si 3 N 4 exhibited the significantly increased the number of energy levels occupied by electrons.Following the incorporation of O atoms (Figure 3c), the Fermi level entered the conduction band, γ-Si3N4 exhibited the characteristics of an n-type semiconductor, and the energy band accommodated electrons.The number of energy levels also increased with a smaller band gap width (Eg = 3.21 eV) (Figures 4 and 5).significantly increased the number of energy levels occupied by electrons.Following the incorporation of O atoms (Figure 3c), the Fermi level entered the conduction band, γ-Si3N4 exhibited the characteristics of an n-type semiconductor, and the energy band accommodated electrons.The number of energy levels also increased with a smaller band gap width (Eg = 3.21 eV) (Figures 4 and 5).significantly increased the number of energy levels occupied by electrons.Following the incorporation of O atoms (Figure 3c), the Fermi level entered the conduction band, γ-Si3N4 exhibited the characteristics of an n-type semiconductor, and the energy band accommodated electrons.The number of energy levels also increased with a smaller band gap width (Eg = 3.21 eV) (Figures 4 and 5).

Bond Order
The bond order, or bond number, provides a numerical expression of the bonding strength of two adjacent atoms in the molecular orbital.It can be derived by calculating the number of electron cloud overlaps between atoms.The higher the bond level, the stronger the atomic bond, and the greater the crystal strength and stability.The bond levels of the six systems were calculated using OLCAO, and the results are shown in Figures 6-8.

Bond Order
The bond order, or bond number, provides a numerical expression of the bonding strength of two adjacent atoms in the molecular orbital.It can be derived by calculating the number of electron cloud overlaps between atoms.The higher the bond level, the stronger the atomic bond, and the greater the crystal strength and stability.The bond levels of the six systems were calculated using OLCAO, and the results are shown in Figures 6-8.

Bond Order
The bond order, or bond number, provides a numerical expression of the bonding strength of two adjacent atoms in the molecular orbital.It can be derived by calculating the number of electron cloud overlaps between atoms.The higher the bond level, the stronger the atomic bond, and the greater the crystal strength and stability.The bond levels of the six systems were calculated using OLCAO, and the results are shown in Figures 6-8.It can be seen from the bond length and bond order distribution plots that the contribution to the bond level of pure Si 3 N 4 (Figure 6) was mainly due to Si-N bonds.The Si-N Crystals 2024, 14, 549 6 of 9 bonds in the β-form (Figure 7) were mainly distributed in the length range of 17.4-1.76,and the bond order fell between 0.325 and 0.375.The Si-N bond lengths in γ-Si 3 N 4 (Figure 8) were ca.1.79 and ca.1.89, with a bond order of ca.0.36 and ca.0.23, associated with siliconnitrogen tetrahedra and silicon-nitrogen octahedra.After doping with C and O atoms, the Si-O bond order was close to 0 for both the βand γ-phases.The results reveal that the presence of lattice oxygen in the Si 3 N 4 crystal results in a loss of bond order, reducing strength and stability.The incorporation of oxygen into the lattice results in a decrease Si-N bond order in the two crystal forms to varying degrees, indicating differing levels of influence on bonding strength.
After the introduction of C atoms, although the bond order is changed slightly, the strength of the C-N bond that is formed is slightly lower than that of the Si-N bond, and the overall bonding strength is reduced.In terms of stability, the introduction of lattice oxygen should be avoided during the preparation of Si 3 N 4 crystals.Moreover, the length of the C-N bond is lower than that of the Si-N bond, and the Si-O bond is longer than the Si-N bond.This can account for the changes in the unit cell parameters after the incorporation of C and O atoms.

Phonon Properties
In order to assess the dynamic stability and thermal conductivity of the six configurations under consideration, the second-order force constant was calculated using first principles, and the associated phonon spectra and phonon density of states were obtained, as shown in Figures 9-12.The phonon spectra did not include any imaginary frequencies, indicating that the systems exhibited dynamic stability.It can be seen from the bond length and bond order distribution plots that the contribution to the bond level of pure Si3N4 (Figure 6) was mainly due to Si-N bonds.The Si-N bonds in the β-form (Figure 7) were mainly distributed in the length range of 17.4-1.76,and the bond order fell between 0.325 and 0.375.The Si-N bond lengths in γ-Si3N4 (Figure 8) were ca.1.79 and ca.1.89, with a bond order of ca.0.36 and ca.0.23, associated with silicon-nitrogen tetrahedra and silicon-nitrogen octahedra.After doping with C and O atoms, the Si-O bond order was close to 0 for both the β-and γ-phases.The results reveal that the presence of lattice oxygen in the Si3N4 crystal results in a loss of bond order, reducing strength and stability.The incorporation of oxygen into the lattice results in a decrease in Si-N bond order in the two crystal forms to varying degrees, indicating differing levels of influence on bonding strength.
After the introduction of C atoms, although the bond order is changed slightly, the strength of the C-N bond that is formed is slightly lower than that of the Si-N bond, and the overall bonding strength is reduced.In terms of stability, the introduction of lattice oxygen should be avoided during the preparation of Si3N4 crystals.Moreover, the length of the C-N bond is lower than that of the Si-N bond, and the Si-O bond is longer than the Si-N bond.This can account for the changes in the unit cell parameters after the incorporation of C and O atoms.

Phonon Properties
In order to assess the dynamic stability and thermal conductivity of the six configurations under consideration, the second-order force constant was calculated using first principles, and the associated phonon spectra and phonon density of states were obtained, as shown in Figures 9-12.The phonon spectra did not include any imaginary frequencies, indicating that the systems exhibited dynamic stability.There were three acoustic phonon branches in the phonon scattering profiles.Doping both β-Si 3 N 4 and γ-Si 3 N 4 with C and O atoms resulted in a reduction in the original acoustic branch frequency due to thermal conduction.The low-frequency phonons were the main contributors to the spectra with the participation of more phonons, resulting in a decrease in the thermal conductivity of the doped crystal lattice relative to the intrinsic β-Si 3 N 4 and γ-Si 3 N 4 crystals.The phonon density profiles for the β-Si 3 N 4 (Figure 11) and γ-Si 3 N 4 (Figure 12) systems exhibited contributions to the acoustic phonons and some low-frequency optical phonons due to the Si and N atoms.The introduction of C and O atoms had little impact on the overall phonon state density, a negligible effect on the low-frequency phonon branch, and a slight effect on the high-frequency region.Overall, the phonon spectra and phonon density of states analyses showed that the thermal conductivity of the intrinsic β-Si 3 N 4 and γ-Si 3 N 4 materials were reduced by C and O atom doping.There were three acoustic phonon branches in the phonon scattering profiles.Doping both β-Si3N4 and γ-Si3N4 with C and O atoms resulted in a reduction in the original acoustic branch frequency due to thermal conduction.The low-frequency phonons were the main contributors to the spectra with the participation of more phonons, resulting in a decrease in the thermal conductivity of the doped crystal lattice relative to the intrinsic β-Si3N4 and γ-Si3N4 crystals.The phonon density profiles for the β-Si3N4 (Figure 11) and γ-Si3N4 (Figure 12) systems exhibited contributions to the acoustic phonons and some lowfrequency optical phonons due to the Si and N atoms.The introduction of C and O atoms had little impact on the overall phonon state density, a negligible effect on the low-frequency phonon branch, and a slight effect on the high-frequency region.Overall, the phonon spectra and phonon density of states analyses showed that the thermal conductivity of the intrinsic β-Si3N4 and γ-Si3N4 materials were reduced by C and O atom doping.

Thermal Properties
The total thermal conductivity (κt) of ceramic materials is equal to the sum of their electrical (κe) and lattice (κL) thermal conductivity, as follows: Here, κe is calculated using the Wiedemann-Franz law, σ is electrical conductivity and T is temperature.L usually takes the Sommerfeld value (L = 2.44 × 10 −8 WΩK −2 ).However, the value of L is not a constant due to the influence of internal phonon-electron interactions.To rigorously evaluate heat conduction in more depth, the following sophisticated theoretical expression developed by Makinson was employed:

Thermal Properties
The total thermal conductivity (κ t ) of ceramic materials is equal to the sum of their electrical (κ e ) and lattice (κ L ) thermal conductivity, as follows: Here, κ e is calculated using the Wiedemann-Franz law, σ is electrical conductivity and T is temperature.L usually takes the Sommerfeld value (L = 2.44 × 10 −8 WΩK −2 ).However, the value of L is not a constant due to the influence of internal phonon-electron interactions.To rigorously evaluate heat conduction in more depth, the following sophisticated theoretical expression developed by Makinson was employed: where k F is the Fermi wave vector, θ D is the Debye temperature, q D is the Debye wave vector, and k F q D = 2 − 1 3 is derived from free electron theory.J n θ D T is related to Debye and is defined as follows: x n e x (e x − 1) 2 dx where R 0 reflects the extent of intrinsic disorder scattering and is the residual resistance at zero temperature, and R e−ph comes from electron-phonon scattering.
3 N 4 unit cell contains 14 atoms, and the unit cell parameters are a = b = 7.607 Å, c = 2.911 Å, α = β = 90 • , and γ = 120 • .The calculations performed in this study used a 2 × 2 × 2 supercell model with a total of 112 atoms, obtained by extending the β-Si 3 N 4 primitive cell by one unit along the basis vector x, y, and z directions.Carbon doping replaces Si atoms in the supercell, whereas oxygen doping replaces N atoms.