Measuring sub-surface spatially varying thermal conductivity of silicon implanted with krypton

Ion irradiation of materials is well known to impact critical material properties. For example, irradiation of semiconductors is an intentional procedure for doping^1,2 and exfoliation of films.³ Conversely, unintentional irradiation effects on thermal, chemical, and structural properties of materials can contribute to catastrophic failure in the areas of nuclear power generation⁴ and outer-space transport.⁵ Ion irradiation of a crystalline material leads to structural disorder ranging from point defects with crystallinity preserved to full amorphization.^6–9 In nuclear materials, radiation can also result in the fission of atoms and creation of new atomic species,¹⁰ which yields compositional disorder in addition to structural. All of these effects will result in a reduction in thermal conductivity.^11–18 For this reason, further insight into spatial variations in a nanoscale structure and accompanying changes in thermal transport processes as a product of ion bombardment are critical toward the engineering of thermally relevant materials.

Various mechanisms of damage will occur when crystalline materials are bombarded with ions. First, ions may knock atoms in the target material off their lattice sites. If the primary knock-on atom (pka) acquires enough momentum, it can further displace other atoms, resulting in a cascade of defects. If the ion continues into the medium, the local chemical composition remains unchanged, despite crystalline disorder being induced. The target may also be heated as a result of bombardment, which can yield some degree of recrystallization or assist in the combining of defects.^14,19 Structural defects of the host material may take many forms. Individually displaced atoms that are unable to return to their initial position may form interstitial/vacancy pairs (Frenkel defects). Cascades between host atoms may also create a dense region of defects, resulting in amorphous pockets. Annealing effects may also serve to combine individual defects, resulting in networks of dislocations and stacking faults.¹⁹ The incident ions may also come to a stop either between lattice sites or in voids left within the host material, forming interstitial or substitutional defects. Regardless of the specific form, all of these defects will serve as additional scattering sites for phonons.^11–16 

The length scales and distribution of ions and defects can be conceptualized as follows. First, the energy and mass of the bombarded ions generally control implantation depth, where heavier atoms require more energy to implant deeper into the surface. Typical depths range from a few nanometers (surface level) to tens of micrometers within normal semiconductor processing conditions. Similarly, bond stiffness controls the distribution of damage. Ions do not implant at all the same depth (referred to as “straggle”), and the concentration of damage is not uniform with depth. Stiffer bonds tend to limit this straggle and result in more localized ion stopping and damage regions. Both ion and defect distributions are often approximated as gaussian;^{1–3,20–26} however, some exceptions apply. For example, ions may experience “channeling” in crystalline targets.²⁷ Depending on the angle of irradiation,²⁰ ions may easily travel between columns of atoms, extending the depth of the tail of the otherwise Gaussian distribution. Additionally, accumulation effects have been observed. Initial damage and implanted ions change how subsequent ions interact with the target. Above a certain threshold, higher doses under the same fluence may result in a shift of the overall defected region to shallower.^25,26 Conversely, if an increased dose is obtained via a higher rate of bombardment rather than by simply increasing the exposure time, this will result in increased heating of the sample and may yield an increased level of recrystallization or combining of defects with dose.¹⁴ 

The software package The Stopping and Range of Ions in Matter (SRIM)²⁸ is frequently used to statistically predict the final ion distribution and the vacancies produced due to ion bombardment; however, it does not account for all effects mentioned above. SRIM accounts for instantaneous initial interactions between the ions’ and target atoms’ electrons and nuclei (termed electronic and nuclear stopping, respectively) but assumes the target material is amorphous. This means it is unable to account for channeling or fluence-dependent annealing or defect recombination effects. To negate the effects of channeling, single-crystal or large-grain materials are typically bombarded with a slight angle applied between the incoming ion beam and the lattice.^20,27,29 SRIM also poorly predicts the sputter yield in certain situations³⁰ and has been shown to yield damage and ion depth distributions that are off by as much as 33%^31,32 for heavy low-energy ions. In most cases, however, it is not off by more than 10%.³³ 

The thermal properties of ion-bombarded materials have been studied in the past in various material systems. For example, Scott et al.¹⁴ used Time Domain Thermoreflectance (TDTR) to measure the thermal conductivity of silicon substrates bombarded with silicon ions. The authors observed a reduction in thermal conductivity with dose, generally following a trend for $ω4$ phonon scattering rates. Their lowest measured silicon thermal conductivity was roughly 1/3 that of their pristine reference. Given the negligible mass difference between the silicon ions and host atoms, the authors also suggested that it is primarily structural changes responsible for their great reduction in thermal conductivity; a point later confirmed through annealing studies with germanium and silicon ions in silicon.^17,29 A similar study was performed on irradiated diamond,¹⁵ showing a 40$×$ reduction in thermal conductivity when irradiated with carbon, oxygen, or nitrogen. A similar reduction was seen regardless of ion, implying that structural damage was the primary mechanism.

When considering thermal property measurements of ion bombarded crystals, the depth-dependent distribution of damage and final ion concentration may complicate the analysis because the spatially varying defect profile may lead to spatially varying thermal conductivities. For nanoscale thermal measurement techniques, the analysis typically assumes uniform material properties. This is true for 3$ω$ or thermoreflectance-based techniques such as Time Domain (TDTR), Frequency Domain (FDTR), or Steady State Thermoreflectance (SSTR),^34–37 where the analysis is performed for thin films or stacks of thin films, with uniform properties within each layer. The measured volumes using these techniques can be on the order of typical ion implantation depths and defect gradient profiles, which pose a challenge for thermal measurements of the ion-irradiated region.

Others have accounted for depth-varying properties³⁸ or heat deposition^39,40 within the TDTR thermal analysis by discretizing the affected region in the through-plane direction. In the former case, this was used to explain frequency-dependent trends in silicon thermal conductivity, where interfacial scattering effects serve to reduce the near-interface thermal conductivity. In the latter, volumetric heating occurs due to a finite optical penetration depth of the pump, which was approximated in the model via surface-applied heat across a series of interfaces.

In this study, we present a modification to the thermal analysis of time-domain thermoreflectance data that account for continuous spatial gradation of thermal conductivity within the measurement volume. We use this to measure krypton-irradiated silicon and explore the spatial distribution of ion-induced defects and their impact on thermal conductivity. We irradiate silicon samples with krypton ions with doses up to and beyond that required to partially amorphize the silicon at the Kr⁺ ions' end of range. The result is a highly varying continuous gradation of irradiation-induced damage, centered a few hundred nanometers below the silicon surface, and fully within the TDTR measurement volume. This invalidates the typical assumption of a spatially invariant thermal conductivity typical for TDTR. The highly varying defect profile also prevents the simple approximation of the irradiated region as a few individual layers. Instead, by discretizing the region and fitting for thermal conductivity as a function of depth (as opposed to the properties of individual layers), we are able to quantify the magnitude and location of the minimum resulting thermal conductivity, along with the net increase in thermal resistance. We further propose the use of this analysis method with thermoreflectance measurements as a non-contact/non-destructive method for sub-surface measurement of ion-induced defects and damage distributions.

Most thermoreflectance techniques, such as Time Domain (TDTR),³⁵ Frequency Domain (FDTR),³⁶ or Steady State Thermoreflectance (SSTR),³⁷ operate under similar fundamental principles. A pump laser beam heats the sample, while a second probe beam is used to measure the changes in reflectivity of the sample due to the pump-induced temperature rise. The pump beam’s intensity is modulated at a set frequency that allows for lock-in detection of relatively small temperature changes (<1 K⁴¹) and even smaller resulting changes in reflectivity (⁠$10−4−10−6K−1$⁴²). If the two beams are pulsed, the time delay between the heating event and the probe measurement can be varied (as is the case for TDTR). For both pulsed and continuous laser cases, the modulation rate can be varied and the frequency response recorded (as in the case of FDTR). Alternatively, the pump power can be varied and the magnitude of the temperature response recorded (as with SSTR). The measured thermoreflectance response is then fit to a thermal model that accounts for an arbitrary stack of spatially uniform layers, each with their own thermal properties and interfacial conductances between them. The full works of Cahill,³⁵ Schmidt,³⁶ and Braun^37,41,43 should be consulted for the complete math for TDTR, FDTR, and SSTR, respectively; however, it should be noted that this model is based on an analytical solution to the heat equation in cylindrical coordinates under periodic heating, with matched boundary conditions between layers, and an adiabatic condition at the back side, so as to relate changes in temperature at the surface to a surface-applied heat flux. This can be represented as a matrix equation [Eq. (1)] relating the surface and backside temperature (T) and heat flux (Q) in Hankel space to each other. A single matrix (M) represents a single homogeneous layer, and a stack of layers is represented by simply taking the product of all layers’ matrices. Matrix elements A,B,C, and D capture the thermal properties of each layer or interface,

An understanding of the volume of the measured region is critical for the measurement of ion-irradiated samples, where the effects of the ions may vary with depth. The size of the measured region at megahertz modulation frequencies is often expressed as

where K is thermal conductivity, $fmod$ is the modulation frequency, and C is the heat capacity.⁴⁴ This is based on the depth at which the temperature decays to $1e$ the surface temperature. Equation (2) should only be considered as a first approximation⁴³ but is useful under normal circumstances with TDTR. One has decreased sensitivity to the properties of materials beyond this depth, which may allow or disallow certain assumptions.

For the case of ion-bombarded samples, if the expected ion stopping range is outside the measured region, one may be able to neglect the spatial variation resulting from the bombardment process¹⁵ and simply treat the sample as spatially uniform in the model used. If the expected stopping range approaches or is within the measured region, one might use two or three discrete layers for their thermal model: a region where ions pass through, a region where ions come to a stop, and an unaffected region.¹⁶ While this may yield acceptable fits of the multilayer thermal model to the experimental data in certain situations (e.g., for low bombardment doses), it should be noted there there is still some amount of variation in thermal properties within each of these regions.

In order to account for the highly variable thermal properties of the bombarded samples within the measured region, we discretize the ion-affected region in our thermal model into an arbitrary number of sub-layers. The math behind the thermal model mentioned previously is identical; however, many more matrices are used to represent the many discretized layers. Each layer’s thermal properties are set according to a user-selected function, with an infinite thermal boundary conductance (zero thermal resistance) between each, as seen in Fig. 1. We can then fit for the function constants instead of or in addition to individual layer properties.

Care must be taken to ensure that one chooses a reasonable function for thermal conductivity and specific heat; however, if the best fitted constants still yield a poor fit (of the thermal model to the collected data), this can be an indication of issues. As a demonstration of this, we have included the raw TDTR data and fit at three modulation frequencies in Fig. 2 below. When we treat the ion-affected region as homogeneous (fitting for an effective average thermal conductivity), we find significant deviation between the curves and data, suggesting this homogeneous assumption is invalid. Using the gradient fitting technique however, the quality of fit of the model to data is greatly improved, finding a residual of below 0.5%. This quality of fit is similar to those on our calibration reference samples (sapphire, silicon dioxide, and silicon, with aluminum transducers).

This strategy should be universally applicable for any cases where there is a significant spatial variation in thermal properties with depth, so long as (1) in-plane uniformity is preserved and (2) an appropriate form of the fitted function can be chosen. In-plane variation was minimized by rastering the ion beam, and any remaining variation is on the order of the ion beam diameter (typically millimeters). As this is significantly larger than the in-plane length scales of our measurement technique (tens of micrometers), these variations can be neglected. Similarly, in-plane variation was observed in TEM (to be discussed later); however, this was on length scales of tens of nanometers, too small to affect the thermal model for the experiment specifically. Robust uncertainty analysis is also required to check for the presence of other acceptable conductivity functions, as will be discussed later. For the case of ion bombarded samples, we chose a Gaussian function [Eq. (3)] based on the extensive experimental^{1–3,20–26} and modeling^25,28 work showing that ion and damage distributions can both usually be treated as Gaussian. Furthermore, we use SRIM to predict ion and damage distributions and do not qualitatively see significant deviations from Gaussian (Fig. 3). We also make the assumption that ion-induced damage will have an additive effect on thermal resistivity based on the premise that (1) bombardment adds impurities in the form of structural and/or substitutional defects and (2) phonon impurity scattering follows Matthiessen’s rule.¹⁴ This yields the following parameterized expression that we expect thermal conductivity to follow

with constants $Rmin$⁠, $zcenter$⁠, and w controlling the center thermal conductivity, location, and through-plane spread, respectively. We are able to fit for these constants rather than, or in addition to, fitting for individual layers’ thermal properties. We also assume the heat capacity is unaffected by ion bombardment, as the final percentage of ions is quite small (<0.04%, see SI), and given the relative small difference in volumetric heat capacity between amorphous and crystalline silicon.^45,46 We also note that given similar sensitivity to the substrate’s thermal conductivity and heat capacity, any slight reduction in volumetric heat capacity due to a lower amorphous density would result in a roughly proportional increase in calculated thermal conductivity.

Note that while fitting for three or more thermal parameters is atypical, a hybrid FDTR-TDTR approach can be taken; TDTR measurements are taken at multiple modulation frequencies and fitted simultaneously,^47–49 allowing the fitting of additional parameters, or fitting with higher precision. In our case, this also serves as a check on our chosen function; if a poor residual is seen at a specific frequency, it can be an indication that the function used is inappropriate. Our approach to simultaneously fit the TDTR data collected using multiple modulation frequencies is described in our prior work.⁴⁹ 

We also compare our discretization approach to the treatment of the ion bombarded system as a three-layer system (aluminum transducer, an ion-affected region, and atop pristine silicon). This three-layer system requires either the fitting for or assumption of at least three unknown parameters: the Al/Si TBC, the affected region’s thickness, and thermal conductivity. By discretizing the system in the manner discussed however, only one additional parameter is introduced, which improves the quality of fit greatly, as seen in Fig. 2.

We measure crystalline silicon (001) that has been bombarded with krypton ions at an energy of 500 keV, at doses of $108$ through $1014$ ions $cm−2$⁠. Bombardment was done on a 3 MeV Pelletron implanter and performed at a few degrees angle off normal, in order to negate channeling effects. All samples were subsequently coated with an 80 nm layer of aluminum to serve as a transducer for thermoreflectance measurements. No surface preparation was performed prior to implantation; however, the surface was cleaned prior to aluminum deposition using our standard procedure (washed with methanol, acetone, IPA, and followed by 30 min $O2$ plasma cleaning). For the analysis, we assume our aluminum to have a volumetric heat capacity of $2.42MJm−3K−1$⁠,⁴⁶ and we measure its thermal conductivity to be $∼$120 W $m−1K−1$ via four point probe sheet resistance measurements. Time Domain Thermoreflectance was performed at three modulation frequencies (8.4, 4.2, and 2.1 MHz), with pump and probe spot diameters of approximately 20 and 10 $μ$m, respectively. Our laser is pulsed at 80 MHz with an 800 nm wavelength. We use SRIM²⁸ software package to estimate the stopping distances of ions, predicting a stopping range of 290 nm with a longitudinal straggle of 65 nm. This estimation is performed using the built-in material library using the Kinchen–Pease formalism. This, along with thermal penetration depth considerations above, suggests that the entire damage region and the ion stopping region fall well within our measured region (400 nm depth).

Transmission Electron Microscopy (TEM) was performed on highest dose samples to validate the damage profile measured thermally via TDTR and gain a fuller understanding of the exact effects of bombardment. Real-space images give qualitative insight into the structure at varying points within the sample, showing features such as atomic density, crystallinity, strain, and the presence of voids. Diffraction images lend insight into the crystallinity of the sample, where a crystalline structure has sharp diffraction spots as the lattice diffracts the electron beam to discrete points, whereas an amorphous material will appear as diffuse rings as electrons are psuedo-randomly scattered.

Cross section TEM samples were made using either a Thermo Fisher Helios or an FEI Nova 600 Nanolab Dual Beam FIB. A platinum protective layer is used to prevent damage due to the gallium ion beam itself, as cross sections are milled to 100 nm or less. Milling was performed at 30 keV initially and cleaning done at 5 keV. We also prepare and image a pristine silicon sample in an identical manner to ensure that the FIB milling procedure itself is not responsible for the defects observed.

TEM was done using a Thermo Fisher Scientific Themis Z-STEM operating at 200 kV and equipped with a Ceta detector, and HRTEM images were taken using an FEI Titan and double tilt stage at 300 kV.

We use TDTR to collect data at three pump modulation frequencies and fit for a gradient in thermal properties at each frequency simultaneously. Where possible, we report the minimum thermal conductivity and the location of this minimum, and we also take the integral of the inverse of the conductivity profile so as to report the net thermal resistivity as added by ion bombardment (Table I). We present an example fit for the thermal conductivity as a function of depth (⁠$1014ionscm−2$ sample), along with SRIM predictions of the ion distribution and damage profiles. We compare this against a TEM image of the sample [Fig. 3(b)], noting the location of the minimum thermal conductivity closely agrees with the centers of the damaged regions as predicted by SRIM or qualitatively observed via TEM. Diffraction images also indicate some degree of amorphization within the ion stopping region, warranting further investigation via high-resolution TEM.

The predicted dose required to amorphize silicon via krypton ion bombardment at 500 keV is $6×1013ionscm−2$ (Ref. 9) (see supplementary material). Based on this and the amorphous rings seen in diffraction images, we perform high-resolution TEM (HRTEM) on the ion end-of-range region of the $1014ionscm−2$ dose sample. We are able to see discrete amorphous pockets surrounded by crystalline regions within this area [outlined in red, Fig. 4(a)]. By performing Fourier transforms on regions of the image, we confirm whether the material is crystalline or amorphous via the presence [Fig. 4(b)] or absence [Fig. 4(c)] of reciprocal lattice peaks. While the development of amorphous pockets has been observed before,^6–8,24,50 the thermal ramifications of this have not been explored to the best of our knowledge, which we will discuss later in this work.

Having successfully fitted for a function of thermal conductivity and noting its apparent agreement with imaging, we must also explore the range of functions that all yield acceptable fits. This is a similar premise to the quantification of uncertainty when fitting for an individual layer or interfacial properties via standard TDTR or FDTR.

There are four things to consider when quantifying uncertainty in thermoreflectance measurements generally.

First, spot-to-spot variability of measurements should be considered, and a large enough set of measurements should be taken to prevent the possibility of random physical variation across the sample from affecting the results. Error from this can be found by simply taking the standard deviation of the best-fit results (x) across multiple measurements, per Eq. (5) for N measurements,

Uncertainty due to random noise in the measured data should also be considered, i.e., statistical noise added to a dataset may result in different results for the fitted parameters. This can be found by simply taking the diagonal of the variance-covariance matrix returned by the least squares fitting algorithm,⁵¹ and in the case of exceptionally clean TDTR data, it is likely to be negligible.

The influence of the assumptions made in the thermal model must also be evaluated. For example, if our measured transducer thermal conductivity is 10% off (or we consider our four-point probe measurements to have an uncertainty of 10%), this error in the model will affect the fitted conductivity of other layers. This can be explored via the Monte Carlo analysis, where a normal distribution for each fixed parameter is set up, and the data are re-fit.⁵¹ Alternatively, and much more computationally efficient, each fitted parameter can simply be perturbed once.³⁷ This is done for each measurement for P parameters, per Eq. (6), and can be combined across measurements through an arithmetic average, as per Eq. (7),

Finally, there is concern of over-fitting if too many free parameters are used.⁵² Just as changing the transducer thermal conductivity affects the fitted results in the example above, it may be the case where manually changing one fitted parameter significantly and re-fitting for the remaining still yields acceptable fits between the model and experimental data. This can be explored by simply testing many combinations of fitted values and checking whether the quality of fit remains below a given threshold.⁵³ The threshold chosen is dependent on the quality of data based on the level at which the model would be in obvious disagreement with the experimental data. This strategy can also be combined across measurements by averaging (Fig. 5).

Finally, multiple methods of computing uncertainty can be combined, as per Eq. (8), for M methods used,

This final approach can be visually represented as a contour plot, where the two parameters form x and y axes, with color and/or contours denoting the quality of fit at each combination of values. All points within a chosen contour represent an acceptable fit of the model to the data, or stated differently, the extremities of the contour represent the measurement’s uncertainty. There may also be combinations of parameters that yield acceptable fits at one modulation frequency, but that may not yield good fits universally. This can be represented as multiple overlapped contour plots, where the overlap of two or more regions can serve to narrow the uncertainty.^53–57 

While this analysis is traditionally performed with the fitted homogeneous material properties as the axes of the contour plots, we can instead use the fitted function constants if fitting a function for thermal conductivity. Neglecting thermal boundary conductance across the aluminum/silicon interface for the moment, we extend this to 3D for our three gaussian parameters. We then consider the within-threshold 3D volume rather than the area of a 2D contour. For analyzing TDTR data taken at multiple frequencies, we then consider the Boolean intersection of each frequency’s within-threshold volume. We present this 3D overlapped contour for our $1014ionscm−2$ dose sample in Fig. 6.

In order to expand this to four or more parameters (e.g., thermal boundary conductance for our system), it may become impractical to construct a full map of the N-dimensional parameter space. Traditionally, we would discretize a range of values for each parameter, trying all combinations in this 2D or 3D grid [$O(nN)$ computation time for N dimensions for n discretized values in each dimension]. This will severely limit the resolution of the map as the number of fitted parameters grows. Instead, we simply perturb one parameter at a time and fit for the rest [$O(n)$ computation time for n discretized values], checking whether the best fit’s residual falls within the same threshold we set previously. While this does not allow clean visualization of uncertainty, it is able to robustly account for an arbitrary number of fitted parameters. Performing this more comprehensive analysis, we are then able to plot all conductivity vs depth functions that yield acceptable fits to the data (Fig. 6).

We observe a nearly 50-fold reduction in thermal conductivity in our $1014ionscm−2$ dose sample, from $120$ to $2.46Wm−1K−1$⁠, and also note the presence of amorphous pockets. This massive reduction in thermal conductivity warrants discussion of the various phonon-scattering mechanisms that could be involved.

Given the presence of pockets of amorphous silicon, we first consider the modified Effective Medium Approximation (EMA) presented by Minnich and Chen.⁵⁸ This gives an approximation for the effective thermal conductivity for a host material filled with nanoparticles of another material and accounts for both the resistance associated with the heat flow between host and particles and boundary scattering within the host material due to the presence of nanoparticles. In general, this approximation is needed if nanoparticles are of similar or smaller size as compared to the mean free paths of the host (>100 nm for silicon^59–61).

We apply the EMA to our system of crystalline silicon interspersed with amorphous silicon pockets in order to relate the effective thermal conductivity (measured) to the expected thermal conductivity of the crystalline silicon regions if the amorphous pockets were not present. We perform this calculation using the following parameters: (1) the volume fraction of amorphous pockets does not exceed 33%, (2) pockets range in size from 5 to 15 nm, as roughly approximated from HRTEM, (3) the thermal conductivity of the amorphous region is $1.5Wm−1K−1$^45,62,63 and (4) the amorphous-crystalline thermal boundary conductance is 1 GW $m−2K−1$⁠.⁶⁰ We also note our insensitivity to TBC in this regime, where TBC as low as $150MWm−2K−1$ or as high as $2GWm−2K−1$ changes predictions from EMA by only $∼$10%. Similarly, a 100% change in the amorphous silicon thermal conductivity in the model changes our result by only $∼$30%. Insensitivity to TBC and low sensitivity to the amorphous silicon conductivity in the EMA model both suggest the heat flow in this regime is dominated by boundary scattering due to nanoparticles, as opposed to the heat flow between and through the host and particles, both mechanisms being captured by EMA.

We calculate that an upper limit of $∼$10 $Wm−1K−1$ for the defected crystalline region is required in order to find $2.46Wm−1K−1$ effective thermal conductivity measured. This implies that our thermal conductivity reduction, from 120 to 2.46 W $m−1K−1$⁠, is not solely due to the presence of nanoparticles and the boundary scattering they introduce; however, this is still a significant contribution. We must explore other scattering mechanisms by which the thermal conductivity might be reduced from $120$ to 10 $Wm−1K−1$ or below. It should be noted that while low crystalline silicon thermal conductivity has been seen before in nanowires and other constrained geometries,^59,64 the EMA calculation already captures size effects within the resulting structure of a host material interspersed by defects.

We next explore the effects of mass and bond strength due to the introduction of krypton ions. It is well known that the introduction of differing masses into a parent crystal serves as point defects, increasing phonon scattering rates, according to

where the scattering rates scale with the square of mass difference $Δ$m and proportionally to the concentration c.⁶⁵ While the mass difference between silicon (28.0855 $gmol−1$⁠) and krypton (83.798 $gmol−1$⁠) is high, the concentration of Kr ions at the highest dose (⁠$1014ionscm−2$⁠) is still only predicted to be around 0.04% (atomic percent). Lacking a comprehensive study on the thermal effects of alloying of silicon and krypton specifically, we instead turn to literature on the well-studied alloying of silicon and germanium (72.64 g $mol−1$⁠). In Si–Ge systems studied both experimentally⁶⁶ and computationally,⁶⁷ tenths⁶⁶ or hundredths^68,69 of a percent concentrations of Ge still yield thermal conductivities above 100 W $m−1K−1$⁠. Despite the slight increase in mass between krypton and germanium, our ion concentrations are still extremely low, and we reject the notion that mass effects could play a significant role in reducing the thermal conductivity of crystalline regions to 10 W $m−1K−1$ or less.

One may also consider the effects of bond strength, where the bonding between Si–Si and Si–Ge atom pairs is expected to be significantly stronger than the interactions between Si and Kr (a noble gas). Ratsifaritana and Klemens⁷⁰ considered the removal of atomic linkages due to vacancies as being equivalent to a point defect with a factor of 2 mass difference, following the expression from Eq. (9), and we take this as an extreme lower limit for the weakened bonding between silicon and krypton. Even considering this extreme limit, concentrations must be on the order of 1% in order to achieve a thermal conductivity reduction to 10 W $m−1K−1$⁠.⁶¹ We again see that significantly higher ion concentrations would be required for us to attribute our thermal conductivity reductions to these effects.

We next consider the influence of damage to the crystalline lattice structure as induced by displacements and/or cascades. SRIM predicts that up to 40% of atoms may be dislocated as a product of bombardment; however, this should be taken as an upper bound due to the possibility of recombination of self-interstitials and vacancies during the bombardment process. Furthermore, a level of 30% structural damage can be taken as the threshold for amorphization of crystalline silicon.⁷ Taking the effects of mass point defects as an approximate analogy, and without going into depth as to specific types of structural defects present, we note that a 30% concentration of point defects could lead to a thermal conductivity as low as $5Wm−1K−1$ for bulk. This effect is far more significant than either mass or bond-strength scattering effects alone, as the passage of each individual ion can introduce far more structural disorder than would result from the mere presence of the ion itself. We also refer to past studies that have made similar observations; post-annealed ion-bombarded samples have seen a near-complete restoration of high thermal conductivity.¹⁷ Similarly, studies in which the ion and the target were the same material (e.g., silicon bombarded with silicon¹⁴ or diamond bombarded with carbon¹⁵) have seen large reductions in thermal conductivity, despite the negligible mass difference between the ion and the target material.

We, thus, attribute our gross reduction in silicon thermal conductivity to both structural defect scattering (from 120 W $m−1K−1$ down to 10 or below for the crystalline regions) and interface scattering effects between our crystalline and amorphous regions (from 10 or less down to our measured 2.46 W $m−1K−1$⁠). EMA’s insensitivity to TBC in this regime also suggests that interfacial scattering dominates rather than heat exchange between the host and nanoparticles. We also reject the notion that mass-scattering or differences due to the effect of Si–Kr bond strength play a significant role in measured thermal conductivity reduction based on the low concentration of krypton ions predicted. Finally, we note that this upper bound for highly defected but still crystalline silicon represents the lowest measured to our knowledge.

In order to test these hypotheses further, we annealed the $1014ionscm−2$ dose sample to $700°C$ for 30 min in air. We do not expect the loss of krypton ions until above $900°C$⁠;⁷¹ however, this temperature should be sufficient to remove the purely structural disorder in silicon^6,24 and potentially result in the creation of Kr bubbles.⁷¹ TEM on the annealed sample (Fig. 7) shows the restoration of crystallinity with some sparse defects. Diffraction images taken at locations that previously contained amorphous pockets also confirm the removal of the majority of structural disorder and boundaries. By removing structural disorder without the loss of ions, we, thus, explore the relative contribution of mass and bonding defects as compared to structural defects and boundary scattering. Thermal measurements demonstrate a complete restoration of thermal conductivity to that of pristine silicon within uncertainty. We are also able to model the system as two homogeneous layers for our thermal analysis (aluminum on crystalline silicon) and achieve exceptional fits. This restoration of thermal conductivity, despite the preservation of ions, confirms the negligible contribution of mass and bonding effects to our reduction in thermal conductivity and supports boundary scattering as the predominant scattering mechanism leading to gross reduction in thermal conductivity observed. While we do not know the exact nature of sparse defects seen via TEM, or whether they are from recrystallization dynamics or segregation, they appear sparse enough to not measurably affect the thermal conductivity.

In this work, we present a development on the Time Domain Thermoreflectance technique, wherein we discretize a region and fit a depth-dependent function for the thermal conductivity rather than fitting for discrete thermal parameters. This allows TDTR measurements on materials containing a steep gradient in properties, as is the case in ion bombardment. We use this technique to measure crystalline silicon bombarded with krypton ions at varying doses, with an ion end of range centered roughly 250  nm below the surface.

We are able to measure the functional, depth-dependent, thermal conductivity distribution, noting a lowest silicon thermal conductivity of 2.46 W $m−1K−1$ in our highest dose sample. We attribute this large reduction to two primary causes. TEM finds amorphous pockets surrounded by crystalline regions as well as the presence of structural disorder (point defects) within the crystalline regions. We, thus, attribute the great reduction in thermal conductivity to a combination of defect scattering (due to structural disorder alone) and surface-scattering effects between amorphous and crystalline regions (effective medium from amorphous nanoparticles). We further anneal the highest-dose sample so as to remove these effects while preserving the ions themselves, measuring a restoration to pristine thermal conductivity. This allows for rejection of mass and bonding differences from the ions themselves as significant contributors to the great reduction in thermal conductivity.

Finally, we explore the limitations of the analysis, finding we may take advantage of the established hybrid TDTR-FDTR approach for reducing our uncertainty and also finding that our sensitivity to the location and value of lowest thermal conductivity is reduced as the dose decreases. Interestingly, our ability to measure net thermal resistance is preserved, however. We also note that accurate resolution of the location of the end of the range region in the highest dose case suggests that this technique could be applicable outside of thermal engineering, where one is interested in measuring the ion-induced damage profile.

See the supplementary material for (1) the calculations of the critical ion bombardment amorphization threshold, (2) predictions of structural damage via SRIM, and (3) additional details on the combining of uncertainty contours for multi-frequency TDTR fitting.

This work was supported under a MURI program through the Office of Naval Research, Grant No. N00014-18-1-2429.

This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE’s National Nuclear Security Administration under Contract No. DE-NA-0003525. The views expressed in the article do not necessarily represent the views of the U.S. DOE or the United States Government.

Patrick E. Hopkins acknowledges support from the Alexander von Humboldt Foundation.

The authors have no conflicts to disclose.

Thomas W. Pfeifer: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Methodology (lead); Project administration (lead); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review and editing (equal). John A. Tomko: Conceptualization (supporting); Formal analysis (supporting); Methodology (supporting); Supervision (supporting); Visualization (supporting); Writing – review and editing (equal). Eric Hoglund: Data curation (supporting); Formal analysis (supporting); Methodology (supporting); Supervision (supporting); Writing – review and editing (supporting). Ethan A. Scott: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Methodology (supporting); Resources (lead); Writing – review and editing (supporting). Khalid Hattar: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Resources (supporting); Supervision (supporting); Writing – review and editing (equal). Kenny Huynh: Data curation (supporting); Formal analysis (supporting); Methodology (supporting); Visualization (supporting); Writing – review and editing (supporting). Michael Liao: Data curation (supporting); Formal analysis (supporting); Methodology (supporting); Visualization (supporting); Writing – review and editing (equal). Mark Goorsky: Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Project administration (supporting); Supervision (supporting); Writing – review and editing (equal). Patrick E. Hopkins: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Methodology (supporting); Project administration (supporting); Resources (supporting); Software (supporting); Supervision (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review and editing (equal).

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