We report on a method by which we can systematically extract spectroscopic information such as isotropic electron–nuclear hyperfine coupling constants from near-zero field magnetoresistance (NZFMR) spectra. The method utilizes a least squares fitting of models developed from the stochastic quantum Liouville equation. We applied our fitting algorithm to two distinct material systems: Si/SiO_{2} metal oxide semiconductor field effect transistors and a-Si:H metal insulator semiconductor capacitors. Our fitted results and hyperfine parameters are in reasonable agreement with existing knowledge of the defects present in the systems. Our work indicates that the NZFMR response and fitting of the NZFMR spectrum via models developed from the stochastic quantum Liouville equation could be a relatively simple yet powerful addition to the family of spin-based techniques used to explore the chemical and structural nature of point defects in semiconductor devices and insulators.

## I. INTRODUCTION

For decades, electron paramagnetic resonance (EPR) has been used to determine the physical and chemical nature of paramagnetic point defects in semiconductor devices and insulators.^{1–7} When studying fully formed micro- and nanoscale electronic devices, EPR is limited in both its sensitivity and its selectivity to defects directly involved in device performance. A closely related technique, electrically detected magnetic resonance (EDMR), takes advantage of spin-dependent currents in these devices to overcome these shortcomings.^{8–10} However, EDMR spectrometers are complex and costly. EDMR is also still limited in its potential to study devices below metallization layers in three-dimensional integrated circuits. A new technique, near-zero field magnetoresistance (NZFMR), has recently been investigated as a new type of spectroscopy that is capable of overcoming these issues to study point defects in fully processed devices.^{11}

The NZFMR measurement is simpler than both EPR and EDMR measurements. The measurement may provide much of the analytical power of the resonance measurement while also being a suitable technique to study point defects in fully processed three-dimensional integrated circuits.^{11} In this work, we discuss a theoretical approach to the analysis of the NZFMR spectrum that can provide values for electron-nuclear hyperfine interactions in multiple systems. These hyperfine values are sufficiently accurate to provide substantial physical insight with some simplifying assumptions. The approach utilizes modeling of the stochastic Liouville equation (SLE). We have developed a nonlinear least squares fitting algorithm based on modeling of the SLE to extract spectroscopic information from NZFMR spectra. We applied the algorithm to two distinct material systems: Si/SiO_{2} metal oxide field effect transistors (MOSFETs) and a-Si:H metal insulator semiconductor (MIS) capacitors.

To better understand NZFMR measurements, a basic understanding of EPR and EDMR is useful. First, consider a free electron in a slowly varying magnetic field, $B0$, and a perpendicular oscillating magnetic field, $B1$, with frequency $\nu $. $B0$ provides an energy splitting between the spin-up and spin-down states of the electron, the Zeeman splitting. The resonance condition is reached when the energy difference between the two spin states is equal to Planck's constant, *h*, times the frequency of the oscillating magnetic field. The resonance condition for the simplest case is of an unpaired electron otherwise isolated from its surroundings is given by

where $ge$ is the Landé g-factor for the free electron and $\mu B$ is the Bohr magneton.

The analytical power of EPR comes from deviations to this resonance response when the electron resides in paramagnetic defect sites in real material systems. The EPR response is altered by the local environment. In most cases, the two most important factors that alter the resonance condition are spin–orbit coupling and electron–nuclear hyperfine interactions. Spin–orbit coupling changes the isotropic free electron g-factor (2.0023, etc.) to an orientation-dependent factor generally expressed as a second rank tensor. Electron–nuclear hyperfine interactions between the electron in the paramagnetic defect site and nearby nuclei with magnetic moments also alter the resonance response. Considering both spin–orbit coupling and electron–nuclear hyperfine interactions with nearby magnetic nuclei, the resonance condition becomes

The free electron $ge$ is replaced by an orientation-dependent *g* that is usually expressed as a second rank tensor, $mi$ is the nuclear spin quantum number of the *i*^{th} nearby magnetic nuclei, and $Ai$ is the electron–nuclear hyperfine coupling due to the nucleus usually expressed as a second rank tensor.

In the EPR measurement, a sample containing paramagnetic defects is placed in a high-Q microwave cavity that is tuned to the cavity's resonance frequency. The cavity sits in an electromagnet that slowly varies the magnetic field across the resonance condition. In classical EPR, at the resonance condition, an absorption of power detected allows for the evaluation of g-tensor and hyperfine tensors.

Conventional EPR measurements are sensitive to about 10^{10} total paramagnetic defects,^{12} greatly limiting its application to the study of defects in state-of-the-art micro- and nanoscale devices. EPR is not necessarily able to directly determine which defects are directly involved in device performance because the technique is sensitive to all paramagnetic centers within the sample. Additionally, the observations of such defects often require measurements in fully processed devices with far fewer than 10^{10} total defects. EDMR is a variation of EPR which eliminates these shortcomings. EDMR, which is measured via changes in a current or voltage, within a device structure has sensitivity more than ten million times greater than that of classical EPR allowing observations of 10^{3} and possibly fewer defects. Since the measurement involves voltage or current within the active volume of the device structure under study, it is exclusively sensitive to defects within the active volume of the device. Since this volume can often be varied via the application of a voltage or voltages, it is also often possible to alter and crudely identify the physical locations of the defects under study. The majority of EDMR studies have utilized one of the two mechanisms: spin-dependent recombination (SDR)^{8,13–21} and spin-dependent trap assisted tunneling (SDTAT).^{10,22–25} Each of these methods will be discussed in detail in the Modeling section of this paper. Optically detected magnetic resonance (ODMR) can offer even higher sensitivities; under some circumstances, single spin detection is possible with ODMR. Although ODMR offers higher sensitivity, it does not provide as direct a connection between paramagnetic defects and the performance of solid-state devices.

In trying to understand defects present in three-dimensional integrated circuits, both classical EPR and EDMR are limited by their use of a microwave or radio frequency magnetic field. The penetration depth of an oscillating magnetic field is a strong function of frequency. As frequencies are increased, penetration depth decreases. Performing EPR and EDMR measurements at lower frequencies could overcome this problem. However, in the construction of three-dimensional integrated circuits, metallic interlayers placed in between the layers of the devices of interest will completely shield these devices from any oscillating magnetic field. Without the perpendicular oscillating magnetic field, both techniques are impossible.^{11}

During a standard EDMR measurement, a change in current can also be measured when the magnetic field is swept through zero: the NZFMR response.^{10,11,23,24,26–30} The NZFMR response occurs with or without the presence of the oscillating magnetic field. The elimination of the oscillating magnetic field makes the measurement much simpler than conventional EPR and EDMR. The small magnetic fields utilized in NZFMR also simplify the measurement. In some systems, the NZFMR response has been shown to closely scale with the amplitude of the EDMR response to the EDMR response in studies of device stressing,^{27} radiation damage,^{11,24,30} changes in temperature,^{23} and bias.^{10,11,23} An experimental setup for the NZFMR measurement is shown in Fig. 1. The device under study is placed between a set of electromagnetic coils that provide slowly varying magnetic field that is swept through zero. A second, smaller set of coils provide magnetic field modulation. The device is biased to create a current that is fed to a current to voltage preamplifier and then to a lock-in amplifier coupled to the modulation circuitry. The spin-dependent current through the devices is measured as a function of the magnetic field: the NZFMR response. Due to the use of lock-in amplification, the observed NZFMR spectrum is an approximate derivative of the change in spin-dependent current. The NZFMR experimental setup is very similar to that of EDMR spectrometers but lacks the resonant microwave or radio frequency source.^{11}

Although some analytical results have been extracted from the NZFMR spectrum, no systematic approach has been available with which to obtain information on the strength of the interaction given a defined spin cluster. We report on a systematic method to extract information about electron–nuclear hyperfine interactions from the NZFMR via a least squares fitting method. We do this using solutions of the SLE. We show that it is possible to extract information about hyperfine parameters utilizing this approach in the two different material systems: MOSFETs and MIS capacitors.

## II. MODELING NZFMR RESPONSES

Several mechanisms have been explored to explain similar magnetic-field effects on current in organic semiconductors.^{31–36} The proposed models for these similar effects in organics focus on the possible spin–spin interactions present in bipolaron transport,^{33,37} electron–hole (e–h) pair recombination,^{34,35} or detrapping of triplet excitons.^{36} These models have had considerable success in describing the spin–spin interactions present in these organic systems. Unfortunately, organic systems are disordered and contain an abundance of nuclear spins, and multiple mechanisms could conceivably be involved in many cases. This makes it difficult to develop models from which NZFMR line shapes can be directly linked to defect structures. Inorganic crystalline semiconductors and some inorganic amorphous semiconductors offer substantial advantages in the analysis. First, the long range order is present in crystalline semiconductors. Second, the short range order is often present in amorphous inorganic semiconductors. Third, popular constituent elements such as silicon and carbon contain few nuclear magnetic moments. Additionally, transport phenomena in inorganic semiconductor devices are generally well understood. We will focus on two models that best align with our understanding of the spin–spin interactions happening in inorganic semiconductor materials and devices in our NZFMR experiments. These models include an adaptation of the e–h pair model^{34,35} for the case of the SDR and an adaptation of the two-site model^{38} for the case of SDTAT. In order to explain these models, we first need to understand the SLE first proposed by Scully and Lamb.^{39}

### A. The stochastic Liouville equation

Employment of the density-matrix formalism is instrumental to describe the statistical nature of an ensemble of spin states. Consider an ensemble of N spin-pairs in their respective spin states, $\Psi n$. The dynamics of the ensemble can be described by the density operator,

In a complete orthonormal basis, $\varphi i$, we can then define the density-matrix as

where elements of $\rho $ describe the probability of finding the spin system in a corresponding basis state. The density-matrix, ** ρ**, completely describes the state of the ensemble that includes both electronic and nuclear spins. The time evolution of that ensemble can be described through the Liouville equation (LE),

The square brackets represent the commutator operation, while $H$ represents the Hamiltonian of the spin system.

The LE represents the non-dissipative, coherent evolution of the spin ensemble through time, without taking interactions with the environment into account. In order to properly describe both the coherent evolution and the environmentally dissipative, incoherent evolution of the spin ensemble, we turn to the SLE,

The first term on the right side of the SLE describes the coherent evolution of spin pairs in the absence of interactions with the environment, exactly as the LE would. The second term introduces a dissipative effect that selectively removes spin pairs from the system at a rate of *k*. The projection operator, $\Lambda $, projects onto the spin subspace from which the projected interactions are removed. The braces in the second term represent the anticommutator. The third term represents a source term that adds spin pairs in random orientations to the system at a rate of *p* and is typically proportional to an identity matrix, $\Gamma $, in non-magnetic systems.^{31}

Equation (6) is able to describe the generation, annihilation, coherent, and non-evolution of spin–spin interactions through time.^{31} This equation also serves as the basis from which the adaptation of the e–h pair model^{34,35} and the adaptation of the two-site model^{38} are formed that best describe SDR and SDTAT, respectively.

### B. Adaptation of the electron–hole pair model

The e–h pair model put forth by Prigodin *et al.*^{35} and the trap-induced magnetoresistance model of Harmon and Flatté^{34} serve as the basis from which we constructed our model that best aligns with SDR taking place at a deep level defect. In this model, as depicted in Fig. 2, we consider an electron in the conduction band, a hole in the valence band, and an unpaired electron in a defect such as a dangling bond associated with a nuclear spin and hyperfine field present. The electron in the conduction band drops to an intermediate state. In the absence of a magnetic field, the singlet and triplet spin state between the electron in the intermediate state and the electron in the deep level defect site will be mixed by the hyperfine field of the nucleus. In the case of a singlet pair, the electron in the intermediate state will fall to the defect and recombine with the hole in the valence band. Recombination removes these carriers, and therefore spins, from the ensemble and reduces device current. In the case of a triplet pair, the electron in the conduction band has only the option to dissociate from the intermediate state and return to the conduction band. Since the electron is unable to recombine with the hole, device current is not altered in the system. The application of a modest external magnetic field suppresses the mixing caused by the hyperfine fields.^{31} As the modest external field is reduced, the recombination rate increases.

The time evolution of a density-matrix that can best describe this spin system can be represented by an adaptation of the SLE proposed by Hansen and Pedersen,^{40}

The first term on the right-hand side describes the coherent evolution of the spin pairs in the absence of the interaction with the environment—the mixing of the states caused by the hyperfine fields. The second term is a dissipative effect that selectively removes singlets (e–h pair) from the system at a rate of $kS$ and the third term is also a dissipative term that represents the selective removal of triplets that have dissociated to the conduction band at a rate of $kD$. $\Lambda S$ and $\Lambda T$ denote singlet and triplet projection operators, respectively. The fourth term represents a source of carriers to the system that are generated at a rate of *p*, while $\Gamma $ represents an identity matrix. Since carriers are being generated and annihilated at constant rates, we are only concerned with the steady state form of the SLE $(d\rho /dt=0)$. The resulting form of the SLE is then

The factor of $18$ in the last term exists because the dimension of * Γ* is 8 such that the total generation rate of the spin pairs is $18pTr(\Gamma )=p$, which is independent of spin dimension.

The Hamiltonian for a single deep level recombination site within the bandgap can be represented by

where the first term represents the Zeeman splitting of electron spins, $S1$ (electron in the conduction band) and $S2$ (electron in the defect site), in the magnetic field, while the second term represents an isotropic electron–nuclear hyperfine interaction, of strength *a*, between the electron, $S2$, and the nucleus, $I$. Explicit representation of the spin operators can be found in the Appendix. This Hamiltonian considers isotropic interactions but can be adapted to accommodate anisotropic interactions. In this model, we assume that the electron in the intermediate state is not interacting with any nuclear spin or hyperfine field.

### C. Adaptation of the two-site model

The two-site model for organic semiconductors put forth by Wagemans *et al.* that considers only two characteristic sites best aligns with our current understanding of SDTAT in dielectrics and insulators.^{38} SDTAT is enabled through variable range hopping and exploits the conservation of spin angular moment from trap-to-trap tunneling events.^{10} Despite the existence of many sites that a carrier may occupy while hopping through the disordered system, the resistance is dominated by one bottleneck pair for which the two-site model appropriately describes.^{33} Figure 3 depicts the spin-dependence of the trap-assisted tunneling event. In order for the electron to successfully tunnel or hop to the other site, the two electron spins must be in a singlet configuration. If they are in a triplet configuration, the hopping event is forbidden. As the resonance condition is reached and SDTAT is detected through EDMR, the transition from triplet to singlet configurations is forced, allowing for previously forbidden tunneling events to take place, and thus, a change in leakage current is measured. SDTAT detected through NZFMR does not involve a resonance induced transition of the spin state but rather a mixing of states.^{11,23,41}

The adapted two-site model, as depicted in Fig. 4, consists of two nuclear sites, $\alpha $ and $\beta $. Both sites have their own respective nuclear spins (both spins assumed $I=12$), and hyperfine fields and thus coupling constants, $a\alpha $ and $a\beta $, when singly occupied with an electron. If site $\alpha $ is unoccupied, an electron can hop from the environment, *e*, to site $\alpha $ at a rate of $rea$. Upon the occupation of site $\alpha $, the electron will precess about the hyperfine field of the present nucleus. At this point, the electron can either hop directly to the environment at a rate of $rae$ or hop to site $\beta $ at a rate of $Psrab$, where $Ps$ is the probability of the electron at site $\alpha $ and the electron at site $\beta $ being in a singlet configuration. After occupation of site $\beta $, the electron can dissociate back into the environment at a rate of $r\beta e$. In this model, the current flow is defined as the total flow of electrons from the environment to site $\alpha $. In order to greatly simplify the arithmetic, we assume that the current limiting rate is not electron hopping from site $\beta $ into the environment; therefore, we can assume $r\beta e=\u221e$.^{31}

The density-matrix representation of this ensemble can be represented by an adapted SLE as originally proposed by Wagemans *et al.*,^{38}

where $\Lambda S$ and $\Lambda T$ are singlet and triplet projection operators, respectively, and $\Gamma $ is an identity matrix that corresponds to a generation of spins at random orientation. The current through the system is the current through site $\alpha $ and is defined as $I=rea(1\u2212Tr(\rho ))$. Since product states (e.g., the doubly occupied state of $\beta $) are not accounted for, $Tr(\rho )\u22641$.^{42–46} It should be noted that since we are interested in steady state conditions, the current through the system is constant, and therefore, the current, *I*, is treated as a constant scalar value. Assuming steady state conditions $(d\rho /dt=0)$ and that traps remain mostly unoccupied $[1\u226bTr(\rho )]$, Eq. (10) becomes

Since there are now four spins included, the dimension of $\Gamma $ is 16 which leads to the factor of $116$ in the last term in the same way we had for spin-dependent recombination involving three spins.

The Hamiltonian of this system can be represented by

where the first term represents the Zeeman splitting of the electron spin, $S1$ and $S2$, while the second and third terms represent an isotropic electron–nuclear hyperfine interaction, of strength *a*, between the nucleus and the electron occupying the site. Explicit representation of the spin operators can be found in the Appendix. Again, this Hamiltonian considers isotropic interactions but can be adapted to accommodate anisotropic interactions.

### D. Solution of the SLE

In both models, we solve for the steady state $(d\rho /dt=0)$ where the equations can be algebraically manipulated to take the form

where

Here, $c1$ and $c2$ represent the respective rate constants that depend on the form of the SLE being solved and *g* represents the generation rate at which spins are being introduced into the system. Equation (13) is in the form of the time-continuous Lyapunov equation. The solution of Eq. (13) and the density-matrix at steady state can be represented as

In both models, the measurable-observable change in current is due to singlet-related events, which makes us concerned with the singlet projection of the density-matrix. The measurable SDR current is, therefore, proportional to the product of the probability of the formation of a singlet and the trace of the singlet population of the density-matrix,

Here, the constant *Z* describes the relative contribution due to a signal in terms of an amplitude scaling factor. The measurable SDTAT current is proportional to the product of the rate at which the sites are filled and the probability of occupation of the sites,

## III. LEAST SQUARES FITTING OF THE SLE

When solving for the density-matrix in Eq. (15), $\rho $ is an inherent function of the rate constants, $c1$ and $c2$, the generation rate, *g*, and the Hamiltonian, $H$. The Hamiltonian is a function of the externally applied magnetic field and hyperfine coupling constant(s). We know the externally applied magnetic field to the system at each point in our measurement, so the Hamiltonian simply becomes a function of the hyperfine coupling constant ($a$ in the case of SDR or $a\alpha $ and $a\beta $ in the case of SDTAT). If we consider the simplest case of a single recombination center with a hyperfine coupling constant, *a*, in the case of SDR, or if we consider two sites with the same hyperfine coupling constant, $a=a\alpha =a\beta $, in the case of SDTAT, then the density-matrix becomes a function of hyperfine coupling constant,$a$, rates $c1$ and $c2$, and generation rate, *g*, so

Since the measurable-observable current, *I*, through the system is proportional to the density-matrix, we find that the measurable-observable current is a function of those same constants,

Given a measured NZFMR spectrum as a function of the magnetic field, we can then fit a simulated spectrum to the given dataset using the appropriate form of the SLE via a non-linear weighted least squares method. The residual error, $r2$, of the least squares method can be defined as

where $wi$ is a weighting factor. Here, the residual square error between the measured spectrum, $IM$, and the simulated spectrum, $IS$, is a function of the hyperfine coupling constant *a*, rates $c1$ and $c2$, and generation, *g*. When the residual error is minimized at some set of values, that is, when

we can gain the information about the defects present as a result of the value corresponding to hyperfine coupling constant(s) at the minimum.

In our fitting algorithm, we defined the weights, $wi$, on the curvature of the spectra (second derivative of the absorption spectra) such that

The second derivative of the NZFMR spectrum highlights areas where important features generally appear. Weighting from the second derivative places higher residual error, or effectively higher importance, on these features for fitting. This forces the minimization routine to fit these sections of the experimental data with greater accuracy instead of neglecting them.

Equation (20) is minimized through an interior point minimization routine outlined by Byrd *et al.*^{47} This type of algorithm is also built into the function fmincon within the global optimization toolbox in MATLAB. The minimization routine is also passed a window of constraints for the values of the hyperfine interactions. We know that these values generally fall in between 0 and 100 mT. Anything larger than 100 mT for the material systems we choose to investigate is physically unrealistic. The rates $c1$ and $c2$, and generation rate, *g*, were left unconstrained.

## IV. EXPERIMENTAL DETAILS

### A. Experimental measurements

Our fitting algorithm was applied to the NZFMR spectrum of two distinct material systems: Si/SiO_{2} MOSFETs and a-Si:H MIS capacitors. Experimental details and biasing schemes varied for each material system.

The NZFMR measurement of the Si/SiO_{2} device was conducted via the DCIV (gated-diode) biasing scheme where the source and drain are shorted together and forward biased and the gate is biased so that the channel resides in depletion.^{48} When the gate bias is swept through depletion, a peak in substrate current is measured. This peak in the substrate current is due to recombination at interfacial defects and is observed when there are roughly equal populations of electrons and holes are present at the interface. The DCIV NZFMR response is measured via the substrate current, with the gate bias held at peak recombination current. This biasing scheme has been used in the past to explore the Si/SiO_{2} interface via DCIV EDMR.^{49} In our case, the gate bias was 0.3 V, and the source and drain bias was −0.33 V. This device was irradiated with a ^{60}Co source to a dose of 1 Mrad(Si) in order to generate interface states. The NZFMR DCIV measurement reported was made after irradiation. No detectable NZFMR DCIV signal was present before irradiation. The Si/SiO_{2} MOSFET used had an oxide thickness of 7.5 nm and a gate area of 41 000 *μ*m^{2}.

The NZFMR measurements on a-Si:H MIS capacitor were made via SDTAT. The a-Si:H was biased with −2 V and the SDTAT NZFMR response was measured via leakage current across the insulator region.^{50} The MIS stack consisted of Ti/10 nm of a-Si:H/p-Si.

## V. RESULTS AND DISCUSSION

### A. Spin-dependent recombination

The top of Fig. 5 depicts a gated diode DCIV NZFMR measurement of a gamma irradiated Si/SiO_{2} MOSFET and its fitted spectrum. Based on magnetic resonance measurements, the MOS interface defects involved are silicon dangling bonds, P_{b} centers.^{3,4,6,7,51–53} We would also expect a fair amount of hydrogen in the interface region due to irradiation damage,^{54} so we would also expect the possibility of an interaction with nearby hydrogen. We constructed our fitting algorithm to consider two separate recombination paths each with their own isotropic hyperfine coupling constant acting as an adjustable parameter with which to fit the experimental data. The first recombination path describes an electron falling from the intermediate state to a deep level P_{b} center that is interacting with distant ^{29}Si, while the second path describes an electron falling from the intermediate state to a deep level P_{b} center that is interacting with nearby hydrogen. The recombination and dissociation rates were tethered together since we assume that these rates should be the same at all P_{b} centers regardless of the interaction taking place with the P_{b} center. The scaling factors that relate to relative contributions were left as independent adjustable parameters. This brought the total number of adjustable parameters to six. The algorithm then fits these six parameters to the experimental data. The adjustable parameters of each path are represented in Table I. Table II contains the isotropic hyperfine coupling constants, *a*, singlet recombination rates, $kS$, pair dissociation rates, $kD$, and relative contribution scaling factor, *Z*, the fitting algorithm produced for the two recombination paths.

P_{b} interactions with
. | a (mT)[ a (MHz)]
. | k_{s} (GHz)
. | k_{D} (GHz)
. | Z (arb. units)
. |
---|---|---|---|---|

Distant ^{29}Si | v_{1} | v_{3} | v_{4} | v_{5} |

Nearby hydrogen | v_{2} | v_{3} | v_{4} | v_{6} |

P_{b} interactions with
. | a (mT)[ a (MHz)]
. | k_{s} (GHz)
. | k_{D} (GHz)
. | Z (arb. units)
. |
---|---|---|---|---|

Distant ^{29}Si | v_{1} | v_{3} | v_{4} | v_{5} |

Nearby hydrogen | v_{2} | v_{3} | v_{4} | v_{6} |

P_{b} interactions with
. | a (mT)[ a (MHz)]
. | k_{s} (GHz)
. | k_{D} (GHz)
. | Z (arb. units)
. | Total residual error of the normalized spectra . |
---|---|---|---|---|---|

Distant ^{29}Si | 1.401 [39.263] | 0.174 | 0.037 | 24.333 | 1.958 |

Nearby hydrogen | 0.206 [5.773] | 0.174 | 0.037 | 1.704 |

P_{b} interactions with
. | a (mT)[ a (MHz)]
. | k_{s} (GHz)
. | k_{D} (GHz)
. | Z (arb. units)
. | Total residual error of the normalized spectra . |
---|---|---|---|---|---|

Distant ^{29}Si | 1.401 [39.263] | 0.174 | 0.037 | 24.333 | 1.958 |

Nearby hydrogen | 0.206 [5.773] | 0.174 | 0.037 | 1.704 |

Our fitting algorithm generated the hyperfine coupling constants to be 1.4 mT and 0.2 mT. A distant ^{29}Si interaction with a P_{b} center has an isotropic hyperfine coupling constant that ranges between 1.4 and 1.6 mT.^{55–57} We find the 1.4 mT electron–nuclear hyperfine constant generated as a result of the fit to be consistent with interactions between P_{b} centers with distant ^{29}Si, while we attribute 0.2 mT to an interaction with nearby hydrogen.

We are also confident of the parameters listed in Table I to be within roughly ±20%. Figure 6 displays the percent difference from the minimum square error between the generated spectrum and the experimental spectrum as a range in colors within the parameter space. The axes of plot represent the parameter space and range from 80% to 120% of the values of best fit for the singlet recombination rates, $kS$, pair dissociation rates, $kD$, and isotropic hyperfine coupling constants, *a*, for the x, y, and z axes, respectively. Two plots were generated such that one hyperfine parameter of the two was fixed at its value of best fit, while the second hyperfine parameters varied. The centers of the figures (100%, 100%, 100%) are the values of best fit as presented in Table I and display the minimum square error. Points further from the center display increased the minimum square error. The outer most corners of the plot on the top and bottom planes represent the extremes and of the parameter ranges and display differences from minimum square errors of 100% of the minimum. We attribute the larger range in the difference from the minimum square error to the modest signal-to-noise ratio (SNR) of the experimental spectrum. The algorithm's primary task is to reduce the residual error between the experimental data and its fit. Lower SNR inherently forces more residual error into the fit, decreasing the confidence of the values of best fit.

From the ratio of the relative contribution scaling factors, we found the relative abundances of the distant ^{29}Si and nearby hydrogen to be about 7% and 93%, respectively. This is consistent with what we would expect at the interface; while the amount of hydrogen is not precisely known, it is reasonable to assume it to be much greater than the number of ^{29}Si sites.

In a recent publication,^{58} Si/SiO_{2} NZFMR and EDMR were examined within a mean field model where hyperfine interactions were assumed to be due to gaussian distribution of classical nuclear spin vectors at each of the electron sites (i.e., shallow and deep states). The work found a typical hyperfine field magnitude at the deep level defect to be 0.5 mT which falls within the range of values determined here. A complete comparison between the two approaches will require the model presented herein to be extended to include EDMR.

### B. Spin-dependent trap assisted tunneling

The top of Fig. 7 shows an NZFMR SDTAT measurement in a-Si:H along with a fitted spectrum. Based on previous high-field EDMR measurements as well as EPR and EDMR literature, the dominating defect is the silicon dangling bond. As reported in the literature, a distant ^{29}Si interaction with a silicon dangling bond in a-Si:H has an isotropic hyperfine coupling constant of about 1.5 mT with an integrated intensity corresponding to a relative ratio of abundance of about 14%.^{59} These particular a-Si:H samples also had a very high concentration of hydrogen of about 20%–30%, so we expect the dangling bond to have some hyperfine interaction with nearby hydrogen. The coupling constant for hydrogen ranges between 0.35 mT and 1.0 mT (potentially higher) depending upon the average distance between hydrogen and the dangling bond.^{60,61} Since interactions with nearby hydrogen are much more likely than interactions with distant ^{29}Si, the fitting algorithm was constructed to consider two separate hopping paths. The first path (labeled path 1 in Tables III and IV) describes electron hopping from a dangling bond interacting with nearby hydrogen to another dangling bond also interacting nearby hydrogen. We would expect that the isotropic hyperfine coupling constant at both sites interacting with hydrogen would be same; that is, we expect that $a\alpha =a\beta $. The second path (labeled path 2 in Table III and Table IV) describes electron hopping from a dangling bond interacting with nearby hydrogen to another dangling bond experiencing an interaction with distant ^{29}Si. These two sites in this path are experiencing two different hyperfine interactions, so we cannot assume their isotropic coupling constants would the same, that is, $a\alpha \u2260a\beta $. Since all of the interactions with nearby hydrogen must have the same isotropic hyperfine coupling constants, all three parameters were tethered together as one adjustable parameter. The isotropic hyperfine coupling constant for distant ^{29}Si was left as another adjustable parameter. We also assume the probability of hopping rates must be the same between paths. That is, $rae$ and $rab$ must be the same for both paths, so $rae$ values were tethered together as one adjustable parameter, while $rab$ values were tethered together as another adjustable parameter. The rate at which the two sites are filled, $rea$, were left as separate independent adjustable parameters. This brought the total number of adjustable parameters to six. The adjustable parameters of each path are represented in Table III. The algorithm then fits these six parameters to the experimental data. Table IV contains the isotropic hyperfine coupling constants, $a\alpha $ and $a\beta $, hopping rates, $rae$ and $rab$, and the rate at which the two sites were filled, $rea$, that the fitting algorithm produced for the two hopping paths. Our fitting algorithm generated the hyperfine coupling constants for the two interactions present in the hopping paths to be about 1.7 mT and 0.3 mT. We find these isotropic electron–nuclear hyperfine interactions to be consistent with dangling bond interactions with distant ^{29}Si and nearby hydrogen.

Hopping paths . | $a\alpha $(mT) [$a\alpha $(MHz)] . | $a\beta $ (mT) [$a\beta $(MHz)] . | r_{ab} (GHz)
. | r_{ae} (GHz)
. | r_{ea} (GHz)
. |
---|---|---|---|---|---|

Path 1 | v_{1} | v_{1} | v_{3} | v_{4} | v_{5} |

Path 2 | v_{1} | v_{2} | v_{3} | v_{4} | v_{6} |

Hopping paths . | $a\alpha $(mT) [$a\alpha $(MHz)] . | $a\beta $ (mT) [$a\beta $(MHz)] . | r_{ab} (GHz)
. | r_{ae} (GHz)
. | r_{ea} (GHz)
. |
---|---|---|---|---|---|

Path 1 | v_{1} | v_{1} | v_{3} | v_{4} | v_{5} |

Path 2 | v_{1} | v_{2} | v_{3} | v_{4} | v_{6} |

Hopping paths . | $a\alpha $(mT) [$a\alpha $(MHz)] . | $a\beta $ (mT) [$a\beta $(MHz)] . | r_{ab} (GHz)
. | r_{ae} (GHz)
. | r_{ea} (GHz)
. | Total residual error of the normalized spectra . |
---|---|---|---|---|---|---|

Path 1 | 0.332 [9.304] | 0.332 [9.304] | 0.128 | 0.048 | 616.488 | 0.990 |

Path 2 | 0.332 [9.304] | 1.721 [48.231] | 0.128 | 0.048 | 238.244 |

Hopping paths . | $a\alpha $(mT) [$a\alpha $(MHz)] . | $a\beta $ (mT) [$a\beta $(MHz)] . | r_{ab} (GHz)
. | r_{ae} (GHz)
. | r_{ea} (GHz)
. | Total residual error of the normalized spectra . |
---|---|---|---|---|---|---|

Path 1 | 0.332 [9.304] | 0.332 [9.304] | 0.128 | 0.048 | 616.488 | 0.990 |

Path 2 | 0.332 [9.304] | 1.721 [48.231] | 0.128 | 0.048 | 238.244 |

Given that the electron starts at a site interacting with nearby hydrogen, the probability that the electron will hop to a site interacting with distant ^{29}Si is about 15%. The probability that the electron will hop to another site interacting nearby hydrogen is 20%–30%. Therefore, the probability of hopping to a site experiencing an interaction with nearby hydrogen is about 1.3–2 times more likely than hopping to a site experiencing an interaction with distant ^{29}Si. The ratio of rates at which the sites are filled from Table IV, $(rea)Path1/(rea)Path2$, indicates that the hop to a site interacting with nearby hydrogen is 2.6 times more likely than a hop to a site experiencing an interaction with distant ^{29}Si. We find agreement between the ratio of the rates at which the sites are filled to be consistent with the two most likely hopping paths.

We are confident of the parameters listed in Table IV to be within at least ±10%. Figure 8 displays the percent difference from the minimum square error of a generated spectrum and the experimental data as a range in colors within the parameter space. The axes of plot represent the parameter space and range from −10% to +10% of hopping rates, $rae$ and $rab$, isotropic hyperfine coupling constants, *a*, for the x, y, and z axes, respectively. Two plots were generated such that one hyperfine parameter of the two was fixed at its value of best fit, while the second hyperfine parameters varied. The centers of the figures (100%, 100%, 100%) are the values of best fit as presented in Table I and display the minimum square error. Points further from the center display the increased minimum square error. The outermost corners of the plot on the top and bottom planes represent the extremes and of the parameter ranges and display differences from minimum square errors off 100% of the minimum. We attribute the smaller range in the difference from the minimum square error to the higher signal-to-noise ratio (SNR) of the experimental spectrum. Since SNR was significantly higher for this experimental data, the confidence of the fit was higher as well.

## VI. CONCLUSIONS

In this work, we show that the NZFMR response contains at least some of the analytical power that has been demonstrated in EPR and EDMR measurements when determining the chemical and physical nature of point defects in semiconductor devices and insulator thin films. NZFMR measurements are straightforward to conduct, and they do not require complex microwave circuity. The NZFMR measurement is also expected to be possible on devices deep within integrated circuits, where EPR and EDMR would not be possible. We discuss the use of two models based on the SLE used to describe the NZFMR response in semiconductors and insulators. With relatively small computational effort, we fit the experimental NZFMR response from two distinct material systems by utilizing knowledge of the present spin systems, the SLE, and a least squares fitting approach. Spin–spin interaction strengths, such as isotropic electron–nuclear hyperfine coupling constants, are extracted from each system and are in good agreement with prior knowledge of the defects present in the systems. Our work indicates that the NZFMR response and fitting of the NZFMR spectrum via the SLE could be a relatively simple yet powerful addition to the family of spin-based techniques used to explore the chemical and physical nature of point defects in semiconductor devices and insulator thin films. We note that performing a simultaneous fit of the results obtained with NZFMR with those from EDMR on the same device may provide still more insight into the effectiveness of NZFMR as an independent tool.

## ACKNOWLEDGMENTS

This work was supported by the Defense Threat Reduction Agency (DTRA) under Award No. HDTRA1-18-0012. The content of the information does not necessarily reflect the position or the policy of the federal government, and no official endorsement should be inferred.

## DATA AVAILABILITY

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

### APPENDIX: OUTLINED CONSTRUCTION OF OPERATORS

The construction of the spin operators, $S1$, $S2$, and $I$, that are described in the section outlining the adaptation of the electron–hole pair model to describe SDR is shown in Eqs. (A1)–(A3), respectively,

The construction of the spin operators, $S1$, $S2$, $I\alpha $, and $I\beta $, that are described in the section outlining the adaptation of the two-site model to describe SDTAT is shown in Eqs. (A4)–(A7), respectively,

In Eqs. (A1)–(A7), $\sigma k$ represents the Pauli matrix of the spin along the k-axis and $\Gamma 2$ represents an identity matrix of dimension 2. The symbol $\u2297$ between two matrices represents the operation of the Kronecker product.

The definitions of singlet, $\Lambda s$, and triplet, $\Lambda T$, operators are defined in Eqs. (A8) and (A9), respectively,

Here, $\Gamma n$ represents an identity matrix of dimension *n* which matches the dimension of the product of two spins. The product of two spins within Eqs. (A8) and (A9) is defined in Eq. (A10),

This definition of the product of two spins also applies to the definition of the product of the electron and nuclear spin as described in Eq. (A11),