Combined treatment of phonon scattering by electrons and point defects explains the thermal conductivity reduction in highly-doped Si

The mechanisms causing the reduction in lattice thermal conductivity in highly P- and B-doped Si are looked into in detail. Scattering rates of phonons by point defects, as well as by electrons, are calculated from first principles. Lattice thermal conductivities are calculated considering these scattering mechanisms both individually and together. It is found that at low carrier concentrations and temperatures phonon scattering by electrons is dominant and can reproduce the experimental thermal conductivity reduction. However, at higher doping concentrations the scattering rates of phonons by point defects dominate the ones by electrons except for the lowest phonon frequencies. Consequently, phonon scattering by point defects contributes substantially to the thermal conductivity reduction in Si at defect concentrations above $10^{19}$ cm$^{-3}$ even at room temperature. Only when, phonon scattering by both point defects and electrons are taken into account, excellent agreement is obtained with the experimental values at all temperatures.


INTRODUCTION
With the characteristic lengths of nanoscale devices approaching the mean free paths of the heat-carrying phonons, 1 the need for a detailed and predictive understanding of thermal conductivity is more stark today than ever. In the last few decades, the usage of highly doped (10 19 − 10 21 cm −3 ) semiconductor materials in electronics and energy devices has become prevalent, serving the purpose of achieving enhanced functional properties. 2 For these systems, predictive models are especially important because the increased phonon scattering due to point defects and additional charge carriers results in a substantial drop in thermal conductivity.
Owing to its high abundance, non-toxicity, and ease of dopability, Si continues to be the linchpin of the semiconductor industry. Highly P-and B-doped Si are routinely used as source/drain materials in transistors to avoid unwanted Schottky junctions. 3 Furthermore, highly-doped Si has also found usage in photovoltaics, 4 microelectromechanical systems, 5 and microelectronics, 6,7 to cite a few applications. In thermoelectric applications the advantage of using highly-doped Si is twofold: 8-10 the thermoelectric figure of merit is, on the one hand, proportional to the electronic power factor, which increases with increasing carrier concentration, and, on the other, inversely proportional to the thermal conductivity.
The lattice contribution to the thermal conductivity (κ ) dominates in Si. Phonon scattering by point defects (PDPS) and by electrons (EPS), were identified as the two main contributors to the κ reduction in highlydoped Si. [9][10][11][12] Zhu et al. 10 reported a ≈36% reduction in κ in fine-grained, highly P-doped Si, and asserted that EPS is the major contributor to the κ reduction. In contrast, Ohishi et al. 9 attributed the κ reduction in single-crystal, highly P-and B-doped Si solely to intrinsic anharmonic phonon-phonon scattering and PDPS.
The aforementioned disagreement highlights the problem of separating scattering mechanisms when employing fitted models based on strongly-simplified assumptions about the underlying phonon band structures and scattering mechanisms. Important progress towards a more predictive treatment of κ in doped Si was made recently by Liao et al., 12 who performed an ab-initio study of nand p-doped Si and showed that, EPS at a carrier concentration of p ≈ 10 21 cm −3 can result in a ≈45% reduction in κ at room temperature. The calculations reproduce how κ is lower in p-doped samples than in n-doped ones, in agreement with the experiments. 13,14 However, they do not capture the magnitude of the reduction observed in B-doped p-type single-crystal Si, which at a doping level of 5 × 10 21 cm −3 amounts to more than 70%. 13 In the present work, we investigate the precise mechanisms responsible for the κ reduction observed in highlydoped Si. We calculate κ by employing the Boltzmann transport equation (BTE) for phonons, using only inputs in the form of interatomic force constants (IFCs) and electron-phonon coupling (EPC) matrix elements obtained from density functional theory. We extend the earlier work on EPS 12 and include also the PDPS from first principles. At high defect concentrations, we find that the PDPS rates dominate the EPS rates at all frequencies except the lowest ones and contribute substantially to the κ reduction at all temperatures. On the other hand, EPS dominates at low defect concentrations due to a fundamentally different frequency and concentration dependence. As a result, a correct quantitative prediction of the κ dependence on defect concentration and temperature is obtained only when both EPS and PDPS are taken into account.

II. METHODOLOGY
Within the relaxation-time approximation, the lattice thermal conductivity tensor can be expressed as 15,16 where α and β run over the Cartesian axes, k B is the Boltzmann constant, V uc is the unit cell volume, n 0 iq , v α iq ,and ω iq are the Bose-Einstein occupancy, the group velocity, and the angular frequency of a phonon mode with wave-vector q and branch index i, respectively. τ 0 iq represents the relaxation time of mode iq and is obtained as: where the lifetime of a phonon of mode iq as limited by the scattering caused by: the three-phonon processes is given by τ 3ph iq , the mass disorder due to isotopes by τ iso iq , the point defects by τ pd iq , and the electrons by τ ep iq . The expressions for the phonon scattering by three-phonon processes and isotopes can be found in Ref. 15.

A. Point-defect phonon scattering
The point defect phonon scattering rates, 1/τ pd iq , can be calculated as: 17 where n def is the volumetric concentration of the point defects and e iq represents an incoming phonon mode. T is the matrix that relates the Green's function of the perturbed lattice to that of the unperturbed lattice. It can be represented in terms of the retarded Green's functions of the unperturbed host lattice, g + (ω 2 ), and the perturbation V as: 18 The matrix element of g + projected on the atom pairs lη and l η is given as: 18 where, e iq (lη) is the eigenvector of mode iq projected on the η-th atom in the l-th unit cell. Moreover, where V M and V K are the mass and IFC perturbation matrices, respectively. Their matrix elements are given by: where M is the mass of the defect atom which replaces a host atom of mass M 0 at site lη. Φ and Φ 0 are the IFC matrices of the defect-laden and host systems, respectively.

B. electron-phonon scattering
The electron-phonon scattering rates can be expressed as: 12,19 where g i mn (k, q) is the EPC matrix element of an interaction process involving a given phonon iq and two charge carriers with band indices m and n and wave-vectors k and k + q, respectively. f nk is the Fermi-Dirac distribution function and ε nk is the eigenenergy of an electron state nk. In non-spin-orbital-coupling or non-magnetic calculations, the formula above must be multiplied by a factor of two to take the electron spin degeneracy into account.
The EPC matrix element can be computed within density functional perturbation theory as: 20 where M η is the atomic mass of the η-th atom in the unit cell, α is the Cartesian direction, and ∂V KS (r)/∂u α iq (0η) is the perturbation of the Kohn-Sham potential with respect to the displacement u α iq (0η).

III. COMPUTATIONAL DETAILS
For calculating the PDPS rates, the total energy as well as the force calculations are carried out using the projector-augmented-wave method 21 as implemented in the VASP code, 22 within the local density approximation (LDA) 23,24 to the exchange and correlation energy. For completeness and comparison to LDA we also calculate the PDPS rates with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional. 25 The volume of the unit cell is relaxed until the energy is converged within 10 −8 eV. The 2 nd -and 3 rd -order IFCs are extracted using 5 × 5 × 5 supercell of the rhombohedral primitive cell containing 250 atoms using just the Γpoint. For the 2 nd -order IFC calculations we use the Phonopy 26 software package and for the 3 rd -order IFCs we use thirdorder.py from the ShengBTE package. 15 The same supercell size was used to calculate the IFC for the defect-laden structures. To compute the PDPS rates, the Greens functions are calculated on a 38×38×38 q-point mesh using the linear tetrahedron method 27 for integration over the Brillouin zone. The scattering rates are then calculated on 35 × 35 × 35 q-point mesh.
For the calculation of EPS rates, the EPC matrix elements are first computed on coarse grids and then interpolated to dense grids with the Wannier function interpolation method. 28,29 The interpolations are performed using Quantum Espresso 30 and the built-in EPW package 31 with norm-conserving pseudopotentials. Likewise, both the LDA and PBE exchange and correlation functionals are considered. The initial k and q grids are both 6×6×6, which are interpolated to 35 × 35 × 35 meshes needed for the thermal conductivity calculations. The energy conservation δ-function is treated by Gaussian function with self-adaptive broadening parameters. 20 Finally, the bulk thermal conductivity is also calculated using the 35 × 35 × 35 q-point mesh with the almaBTE 32 package. Due to the dense q-mesh the thermal conductivity is converged down to 40 K. Si defects, which are the most common doping elements and also those studied experimentally. The PDPS has a trivial dependence on the defect concentration [Eq. (3)] and is shown for only one defect concentration (10 21 cm −3 ). In a rigid band approximation, the EPS is dependent only on the carrier concentration through the distribution function in Eq. (8). The EPS is shown for three different carrier concentrations in Figure 1 and the corresponding chemical potentials are shown in Fig. 2 (inset). Our LDA EPS rates compare well with the PBE rates reported by Liao et al.. 12 The LDA rates are slightly higher than the PBE ones which results in a correspondingly higher κ reduction due to EPS in LDA as compared to PBE. However, this does not have any major effect on the interpretation of our results.

IV. RESULTS AND DISCUSSION
If we first consider the low frequency (ω < 12 rad/ps) behavior, the PDPS rates exhibit a simple Rayleigh ω 4 behavior and the EPS rates for a given carrier concentration are close to being independent of ω. Plotting the EPS rates as a function of the electronic density of states (DOS) at the electron chemical potential corresponding to a given doping shows an almost linear dependence (Fig. 2). The EPS rates thus behave in accordance with a simple τ −1 ∝ n(ε) model 33 at low frequencies. With a n(ε) ∝ ε 1/2 behavior of the electronic DOS, the lowfrequency EPS rates will scale approximately as n 1/3 def with the carrier/defect concentration as opposed to the linear scaling of the PDPS rates evident from Eq. (3). These simple relations are in accordance with the expectations that EPS will dominate over PDPS at low temperatures and defect concentrations while PDPS will become increasingly important at high defect concentrations.
At the same time, it is clear from Fig. 1 that for ω > 12 rad/ps the calculated rates deviate substantially from the aforementioned simple relations. The cumulative κ plots in the respective top panels in Fig. 1 (blue curves) illustrate that modes with ω > 12 rad/ps carry about two-thirds of the heat at 300 K. Simply extrapolating the low frequency behavior would lead to a strong overestimation of the scattering. This is especially so for the EPS rates where a simple extrapolation would result in a strong overestimation of the predicted κ suppression. The cumulative plots in the top panels of Fig. 1 also show for PDPS (dotted lines), at a large defect concentration, that the contribution to κ for frequencies higher than 20 rads/ps is only ≈15% in case of B (−1) Si defect and ≈40% in case of P (+1) Si . In contrast, the EPS causes a majority reduction in κ by frequencies below 20 rads/ps, for both P-and B-doping.
Next, in Fig. 3, we look into the individual contributions from the EPS (dashed lines) and PDPS (dotted lines) to the room temperature thermal conductivity reduction for increasing doping concentrations and compare them to the experimental κ data from Slack. 13 In accordance with the analysis of the scattering rates, EPS dominates the reduction of κ at low carrier concentrations. In the case of n-doped Si, EPS alone is almost enough to explain the experimental point at 2 × 10 19 cm −3 . However, at higher carrier concentrations EPS alone clearly underestimates both the absolute reduction of κ and the trend. Interestingly, PDPS captures the trend correctly for large defect concentrations, but also underestimates the absolute κ reduction. At a doping concentration of 10 21 cm −3 , the EPS and PDPS contribute almost equally in κ reduction for B doping. Even though both EPS and PDPS contribute substantially to κ reduction in Si, neither can explain the absolute reduction in κ on its own. Only when both are taken into consideration in Eq. (2) is the experimentally observed reduction of κ reproduced. This is shown by the black and red solid lines in Fig. 3. Apart from a slight underestimation of κ in case of B-doping, the calculated value of κ considering both EPS and PDPS agree very well with the experimental values available for concentrations ∼10 19 − 10 21 cm −3 both in value and trend, Fig. 3.
Besides the work of Slack et al. 13 used for Fig. 3, there is a general lack of systematic experimental data on the thermal conductivity of single-crystal Si with varying doping concentrations. In order to gain further confidence in our results, we compare them to the more recent  experimental data from Ref. 14, which were measured on single-crystal Si films. Even at low defect concentrations, n def = 10 17 − 10 18 cm −3 , these samples exhibit a substantially lower κ than the bulk samples. 14 However, a good agreement with the experimental curves can be obtained by adding a simple boundary scattering term, 1/τ B iq = |v iq |/L with L = 10 µm, to Eq. (2) to emulate the effect of a film, as seen in Fig. 4 (purple line). Adding now the effect of the PDPS and EPS we calculate the variation of κ as a function of temperature for both B-and P-doped Si films, shown in Fig. 4. For Bdoping, when we include the PDPS and EPS along with the boundary scattering, we see that there is only a slight reduction from the purple line for the 10 18 cm −3 doping 10 16   (10 20 and 10 21 cm −3 ) as there is no experimental thermal conductivity data available for such high doping levels. These are shown in Fig. 5. At room temperature, we find more than 60% reduction as compared to the bulk value at 10 20 cm −3 doping level and 90% at 10 21 cm −3 for B doping, and 40% and 80% for P doping, respectively. Fig. 5 also shows the individual contributions to the κ reduction by EPS and PDPS in the case of B doping at a doping concentration of 10 21 cm −3 , revealing their characteristic effects on the thermal conductivity over the temperature range. The reduction in κ caused by PDPS overtakes the one by EPS at ≈250 K. This behavior highlights the importance of PDPS at high doping concentrations in Si and also emphasizes that, at temperatures higher than 300 K, PDPS is the dominant scattering mechanism besides the intrinsic anharmonic scattering. We have also performed calculations including an L =10 µm boundary scattering term, however at such high doping concentrations the EPS and PDPS dominate the boundary scattering throughout the temperature range and only a small effect at very low temperatures was found on the calculated κ .
It is important to note how the present case is very different from our previous work on SiC. 34,35 In the present work, we observe a significant contribution of both EPS and PDPS to κ of a highly-doped system; whereas in our previous work we observed that PDPS was sufficient on its own to correctly predict the κ of B-doped cubic SiC owing to the resonant phonon scattering that boron causes. 34,35 The resonant scattering was at least one to two orders of magnitude higher than that caused by other defects and resulted in a drastic reduction (approximately two orders of magnitude at room temperature) in the thermal conductivity even at relatively mod- est defect concentrations (≈ 10 20 cm −3 ). Boron does not cause resonant scattering in Si and therefore the contribution of both EPS and PDPS are comparable.

V. CONCLUSIONS
We have calculated the ab-initio lattice thermal conductivity of highly P-and B-doped Si considering phonon scattering caused by three phonon processes, isotopes, point defects, and electrons. We illustrate that at low doping concentrations EPS causes a higher reduction in κ as compared to PDPS owing to a near constant behavior of EPS rates at low frequencies. At low concentrations EPS alone is sufficient to reproduce the absolute reduction in κ . However, at high doping concentrations it fails to reproduce the absolute values and the trend of κ reduction. PDPS has a substantial contribution to κ reduction at high doping concentrations and neither EPS nor PDPS is sufficient on its own to reproduce the experimental κ values. Only when the effect of all the scattering mechanisms are considered together, we get a good agreement with the experimental values across the temperature range. At a doping concentration of 10 21 cm 3 , we observe almost 90% reduction in the room temperature thermal conductivity of B-doped Si as compared to the bulk, whereas, 80% reduction in case of P doping is found. This is mainly because of the higher PDPS rates of the B defects than those of P defects. Neither EPS nor PDPS can be captured by conventional parameterized models. Together with our previous works on doping diamond, SiC, GaN and GaAs, 34,36-38 we show that, at temperatures above 300 K, PDPS is the most dominant scattering mechanism in highly B-and P-doped Si and cannot be neglected.

VI. DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

VII. ACKNOWLEDGMENTS
The authors acknowledge support from the European Union's Horizon 2020 Research and Innovation Ac-