On this point, one can see the great value of a TCAD model to the conversation, as with a TCAD model, one can capture these physical effects in a direct way—one simply inserts a spatially varying potential with the correct properties into the simulation. Thus, no ad hoc assumptions or hand-wavy insertions of percolation results from resistive networks or effective media are needed. Instead, the microscopic TCAD model presented here is fundamentally agnostic to these often difficult-to-justify approximations.
2.1. Correlated Gaussian Random Field Mobility Edges
The fundamental task at hand is to generate a spatially varying band-edge, , subject to the following two conditions:
Globally Gaussian distributed: Globally, the finite set of band-edge values,
, of the mesh, corresponding to certain values of the band-edge at each vertex
i and integrated over the Voronoi or box volume associated with
i, must follow the Gaussian distribution:
with standard deviation
and a mean of
.
Locally spatially correlated: Band-edge values of nearby mesh points should not vary independently but, rather, exhibit a spatial (auto)correlation such that:
where
represents the position of vertex
i, and
is a correlation function that must be given explicitly. By assuming that the correlation function only takes in a scalar variable representing the distance between points
i and
j, it is being assumed that the correlations are
homogeneous and, thus, have no explicit spatial dependence (i.e., there are no “special” places and the correlations are the same everywhere) and are
isotropic and, thus, have no explicit dependence on orientation or angle. The possible specific forms of the correlation function,
, are discussed in
Appendix A.
These combined properties are demonstrated in
Figure 2.
Note that, in general, we would like an approach based on unstructured meshes, but the approach developed here relies on the use of Fourier transforms, which can be efficiently generated—via the use of the fast Fourier transform (FFT) algorithm—only on regular grids. One could attempt to construct an unstructured algorithm, but here we take the simpler approach of generating a field on a regular, finely meshed grid and using interpolation to map the values back to the unstructured one. To us, this approach is quite justifiable, as this band-edge generation step must only be done once at the beginning of the TCAD simulation and, thus, even though it relies on fine-gridding that is then potentially discarded, the total contribution to the numerical burden of the full simulation is typically minimal.
The first requirement stated above amounts to generating a so-called Gaussian random field (GRF), which is a well-known object in the field of statistics. The second requirement can be obtained by exploiting the Wiener–Khinchin (WK) theorem [
20,
21]. As both these aspects are well-studied topics in statistics, we will make no attempt to prove or justify them here. Instead, for completeness, in
Appendix A, we provide a derivation of and the motivation for an algorithmic approach for correlated GRF generation, but in order to keep the main body of this manuscript concise, here we will only give the final resulting algorithm.
In summary, the procedure is as follows (shown for 3D, but the approach is also suitable for 2D):
Generate a field of
values either according to Equation (
A9) directly or through generating random Gaussian numbers of unit variance,
, in real-space and performing a discrete Fourier transform of the values. We call the resulting field
. (Input:
(if Equation (
A9) used),
; Output:
).
Generate an grid of k-values, where , where has been shifted to include negative and positive values centered at . Many numerical environments have built-in functions for this (for example, 2 * numpy.pi * numpy.fft.fftfreq((Nx,Ny,Nz)) in numpy). We call this grid or simply k though it is three-dimensional. (Input:, and ; Output: k).
Scale the field
by
, where
is either Equation (
A3) or (
A4) (depending on the desired correlation) evaluated at each point
k. This scaled field we call
. (Input:
; Output:
).
Perform an inverse discrete Fourier transform of to yield . (Input: ; Output: ).
Rescale and shift
according to Equation (
A10) to give the final spatially correlated field. (Input:
; Output:
(normalized)).
Once the above algorithm is performed (again, see
Appendix A for a detailed derivation of this algorithm), the result is a field
with the desired properties.
Figure 3 and
Figure 4 demonstrate this for a sample GRF generated using this algorithm.
Figure 3a shows that the generated GRF does indeed manifest the desired spatial correlations: in this case, exponentially decaying correlations with a correlation length of 20 nm are used, and the expected behavior is shown as a dotted black line.
Figure 3b shows that the band-edge values at each vertex also follow the probability distribution dictated by Equation (1).
Figure 4a considers the differences in the means and standard deviations of the calculated mobility for 10 different randomly generated GRF samples for the case of exponentially correlated and uncorrelated spatially varying band-edges as a function of grid size for a fixed system size. The model parameters and procedure used for the determination of the mobility are discussed in
Section 2.3. It can be seen that a correlated model approaches saturation to a final finite value with minimal variation with sufficiently small gridding, whereas the uncorrelated approach saturates to zero mobility. Based on these results, it would seem that for mesh resolutions below approximately half the correlation length (10 nm here), the mobility reaches an acceptably saturated value with quite low variation. This implies that a mesh resolution of
or finer is desirable for simulation.
Figure 4b is the same data as (a), only plotted with mobility values and standard deviations normalized to the large grid size (i.e., right-most) value.
2.2. The Quasi-Fermi Level Transport (QFT) Model
In our model, we assume drift diffusion (DD) transport through the spatially varying system. Furthermore, we assume a simple constant field- and temperature-independent mobility,
, at all mesh points. Note that, as will be seen in
Section 2.3, this certainly does not mean that the a-IGZO layer as a whole will exhibit a constant field- and temperature-independent mobility. Rather, the spatially varying band-edge will inject complex transport physics driven by percolation effects into the transport characteristics of the device, and thus, the
emergent, aggregate mobility of the whole device will indeed exhibit field and temperature dependence.
Furthermore, the assumptions of DD transport and a constant mobility are not essential ingredients to the model we present, and thus, other transport models and more complex mobility models can easily be used along with the spatially varying potentials of
Section 2.1 and the band-tails that will be described in
Section 2.4 to expand the basic model presented here. In fact, the use of a more complex mobility model with mobility degradation near the oxide surface will be considered later in
Section 3.2.
However, even the use of the relatively straightforward and ubiquitous DD model requires special care when it comes to a-IGZO. In general, TCAD simulation in situations where carrier concentrations are very low presents challenges for numerical convergence. This is true for wide-band semiconductors such as a-IGZO and for systems at cryogenic temperatures where there are only very few thermally excited carriers. This issue is then especially bad if one wants to consider a wide-band-gap semiconductor like a-IGZO at cryogenic temperatures. Yet the study of the mobility of a-IGZO as a function of temperature is frequently a focus of experimental studies and thus of great interest for supportive TCAD simulations.
There are two main reasons for the numerical difficulties associated with lower carrier densities, and they can be understood by considering the fundamental current equation of the drift diffusion model:
where
is the current,
is the spatially varying band-edge,
q is the electron charge,
is the Boltzmann constant, and
T is the temperature. The value
is the carrier concentration, which may, in principle, be calculated using either Maxwell–Boltzmann or—as is more appropriate for a-IGZO, where the Fermi level frequently enters the conduction band—Fermi–Dirac statistics. However, in reality, in the usual implementation of drift diffusion, by using the so-called Scharfetter–Gummel (SG) scheme [
22], Boltzmann statistics become explicitly baked-in to the model.
The first issue is what has been called
catastrophic cancellation [
17], where the first drift term and the second diffusion term are comparatively large but very close in absolute value such that they agree for many decimal places and only differ in the deep trailing decimal values. This creates a problem as the mantissa (i.e., decimal portion) of the
double variables frequently used in numerical computation may not be sufficient to capture this finite, non-zero difference in nearly identical numbers as distinct from zero.
This first problem can be somewhat addressed by using a larger numerical datatype, such as a long double, and that is what is done in our model as well. However, there is an even greater issue, and that is that at low carrier concentrations, as the carrier concentration exponentially depends on the band-edge and/or Fermi level in a form , as T becomes small or becomes large, small changes in can lead to dramatic changes in carrier concentration. This creates a great problem for meshing in TCAD, as all carrier concentration gradients will become sharper and sharper and more abrupt and thus require a finer and finer mesh in order to be spatially resolved as the temperature becomes lower or the band-gap becomes wider. Such ultra-fine meshes are obviously numerically cumbersome.
This issue can be greatly alleviated by reformulating the basic drift diffusion equation into a form we call “Quasi-Fermi-Level Transport” or QFT. This approach is expanded upon in greater detail in a previous publication by some of the authors [
17] and so will not be fully described here. However, the salient idea of the QFT approach is to make the spatially varying quasi-Fermi level,
, the key unknown quantity to solve for:
This is in contrast to the regular DD Equation (
5), where the carrier concentration,
, and the band-edge (or electrical potential),
, are the unknown quantities. In DD, as there are two different quantities, one generally solves the equation iteratively until self-consistency is obtained between
and
. However, in QFT, where the Fermi level is the solved quantity, this self-consistency is traded in for the basic equation being explicitly non-linear and with the coefficient of the gradient term being some complex function of
(i.e., the Maxwell–Boltzmann or Fermi–Dirac equation weighted by some density of states).
There is also the subtle detail of how the current should be evaluated between mesh points. For DD simulation, as was said, the well-known Scharfetter–Gummel scheme is used and attempts to formulate the relevant current equations in a way that explicitly considers the exponential relationship between carrier concentration and potential, as assumed in Boltzmann statistics, in order to yield greatly enhanced numerical convergence. A similar scheme is also necessary in QFT, but luckily, such a scheme can be created with only small changes to the regular SG and can be extended to permit Fermi–Dirac statistics. Details about this and any further details about the properties of QFT are given in [
17].
Thus, our a-IGZO model assumes drift diffusion transport but casts the key equations in the QFT form. This, combined with the use of
long double precision, allows the model to achieve much better convergence than regular DD. This is demonstrated in
Figure 5, where the mobility (using the model parameters derived later in
Section 2.3 and given in
Table 1) is shown as a function of the inverse temperature for DD and QFT. As can be seen, at temperatures less than ∼200 K, the DD simulation fails to converge, whereas convergence can be obtained with QFT down to 100 K. Below 100 K, however, even the QFT results become increasingly noisy for the mesh size used (not shown). This could be improved by considering a higher mesh resolution, but experimental results on a-IGZO in the literature rarely go below 100 K, so we consider that limit to be sufficient for our purposes here.
2.3. Validation against Experimental Results
In order to demonstrate the physical accuracy of the TCAD model developed here, it must be validated against experimental results. Although there are a great number of experimental results available, here we make use of those from Germs et al. [
5], which have already previously been used in the literature for model validation: both in the original paper itself and in [
11].
As described in Germs et al. [
5], the measured device consists of a 10 nm-thick layer of a-IGZO atop a ∼200 nm-thick layer of SiO
2, making for a back-gated device. In the process of fitting, the value of the a-IGZO layer was held fixed, but we found it beneficial to allow for the value of the oxide thickness,
, to vary by some percent and to thus act as a strongly constrained fittable parameter. Alternatively, this allowed for small variations to be considered to encapsulate any intermixing of the a-IGZO and oxide to create a thin intermediate layer of different permittivity. The a-IGZO film was then contacted on both sides by 25 nm-thick gold contacts with 5 nm-thick titanium adhesion layers—though within the TCAD simulation, they are simply treated as perfect conductors (note that this means that they have zero-resistance and is independent of whether they contact the channel in an Ohmic or Schottky way). The experimental devices had lengths and widths in the 100s of micrometers; however, it is both unnecessary and computationally prohibitive to simulate such large films while maintaining microscopic detail at the scale of 10s of nanometers. Therefore, the simulated device has a much smaller length and width of 200 nm, with the grid resolved to a scale of 2.5 nm. Numerical experimentation shows (as will be shown in
Section 3.1) that this reduced device size has a negligible effect on the results provided that the device is still many times larger than the correlation length of the spatial fluctuations. A 2D cross-section of the final simulated device is shown in
Figure 6b, though, note that as percolation physics plays a central role in transport, a 2D device has fundamentally different percolative behavior than a 3D device, and all simulations are thus 3D. The question of whether an effective 2D or “compact” model can be adapted from this 3D model is considered to be beyond the scope of this publication.
In order to fit the data, there are a number of parameters that must be set. The first issue is that a-IGZO devices such as those in Germ et al. [
5] have no doped
n- or
p-junctions and, thus, are fundamentally junction-less or Schottky-based devices. Thus, Schottky contacts must be assumed to replicate accurate behavior, with the workfunction differences between the source (S) and drain (D)—which are assumed to have the same workfunction as they are the same material—and the a-IGZO channel (C) potentially playing an important role. We notate this value as
, and the case for which it is zero corresponds to perfect ideal Ohmic contacts.
Assigning a single authoritative value to
is very tricky for a number of reasons. The first of these is that in the literature, both the values of the band gap,
, and the electron affinity,
, of a-IGZO [
23,
24] fluctuate quite a bit, with typical values ranging from 3.0–3.5 eV and 4.0–4.5 eV, respectively. This uncertainty is compounded by the fact that the workfunction of the gold contacts, which is typically in the neighborhood of ∼5.0 eV, varies significantly depending on the process and Fermi-level pinning.
However, a greater issue is that even in the absence of Fermi pinning much of the “textbook” intuition around the expected band bending in a Schottky contact that is observed in conventional semiconductors like silicon does not apply in ultra-low-carrier-concentration materials like a-IGZO. For example, conventionally, one has a Schottky contact boundary dictated by the difference between the work function and the electron affinity of the metal and semiconductor, respectively, which then transitions as one moves deeper into the semiconductor to bulk semiconductor behavior. However, in a-IGZO, one will never observe bulk behavior.
The reason for this is because the length-scale of this transition from a Schottky boundary to bulk behavior is related to the Debye length of the semiconductor, which is based not only on the permittivity (here, a permittivity, , of 10.0 was assumed) but also on , where n is the electron concentration. In undoped a-IGZO with a band-gap of 3.2 eV—and, thus, an intrinsic carrier concentration that is times less than that of silicon—this Debye length at room temperature may actually be in the range of kilometers. Therefore, the difference between the average band-edge and the Fermi level everywhere in the device is entirely determined by the contacts in a-IGZO: the channel is extremely electrostatically flat (when neglecting the band-edge fluctuations), and any Schottky barrier mainly represents an “effective” value resulting from the competing electrostatics of the gate and source/drain contacts.
However, perhaps the biggest issue with assigning a definitive value to is the fact that the “doping concentration” of a-IGZO can be a somewhat nebulous concept: a-IGZO is generally not doped with conventional dopants but, rather, the oxygen vacancies present within the film as a result of fabrication processes dictate the unbiased intrinsic Fermi level. Furthermore, for an “undoped” a-IGZO film, even a small concentration of the dopant-behaving vacancies can cause a substantial upward-in-energy movement of the Fermi level towards the conduction band-edge due to the ultra-low intrinsic carrier concentration. Thus, the “intrinsic” Fermi level of the experimental film cannot reasonably be expected to actually be in the mid-gap if any oxygen vacancies may be present, and this will manifest as a reduction in as the average band-edge is pulled down (i.e., the Fermi level goes up due to dopant-behaving vacancies).
Thus, we assumed a fixed band gap value of 3.2 eV, and in light of all these considerations, we treat the contact workfunction as a fittable parameter encapsulating the physics of both the doping level and the material energetic difference and allow it to take a fairly broad range. Furthermore, due to the wide band gap and the system always only operating in the n-type regime, the simulation of holes was neglected.
In addition to the contact workfunctions, , a value for the workfunction difference between the doped-Si bottom gate (G) and the a-IGZO channel and/or source/drain contacts must be set, which will determine the threshold voltage, , of the “turn-on” of the device. Here, we choose to denote this quantity as and define it as the offset relative to the source/drain contacts rather than defining it relative to the channel due to the extremely weak electrostatic influence of the channel itself and the lack of any identifiable “bulk” region being present anywhere in the device. As before, this is simply treated as a fittable parameter, and, unlike the Schottky contact workfunction value, we can expect to potentially vary wildly in value depending on the number of traps present at both the gate and the a-IGZO sides and based on other process conditions.
As a result of these workfunctions, there is a certain non-uniform
baseline band-edge reflecting only the Schottky contacts and gate workfunction, onto which the additional Gaussian band-edge variation is added. A 1D cross-section from the source to the drain of these two band-edges is shown in
Figure 7.
As the mobility being modeled is an effective mobility for the entire system as controlled by the percolative impediment a carrier faces as it traverses the system, assigning a single-valued mobility to the system cannot be done by doing some sort of integral of the local mobility at all points. Rather, here we define the final TCAD mobility,
, in a manner directly identical to how it is extracted from actual experimental measurements on FET devices as:
where
L is the system length and
W is its width (here, in both cases, 200 nm), and
is the capacitance of the oxide layer.
and
are the drain current and bias, respectively (in all simulations, a
of 0.01 V was used), and
is the gate voltage, which varied from 0 V to 20 V.
Finally, with the geometry- and device-specific parameters set, what remains is to determine the microscopic model properties of the a-IGZO film itself. These include: the correlation length of the band-edge variation, ; the standard deviation of the band-edge variation, ; and the constant bulk mobility, . The average value of the band-edge, , is not treated as a fittable parameter but is instead dictated by the band gap, which, as mentioned previously, is given the fixed value of 3.2 eV.
Although the correlation length can be treated as a fittable parameter, instead, we fix its value to be 10 nm. This corresponds to the typical size of a spatially varying band “hill” or “valley” of ∼20 nm in diameter (i.e., the correlation length effectively defines the radius of influence about a given point). We believe this to be a reasonable value that is justifiable by, for example, some scanning tunneling microscopy (SPM) measurements [
25] on a-IGZO films, where characteristic structures of roughly that size were found. However, we also admit that studies that reveal such information are rare and can often vary wildly in their observed fluctuation size, and thus, it may be necessary to treat
as a fittable parameter in general.
Thus, the set of fixed values is
,
,
, and
, and the set of fittable parameters is
,
,
,
, and
. Experimental reference values were extracted directly from Figure 4 of [
5], and an optimization algorithm was run to optimize the TCAD model to match that data set. The final optimal values are listed in
Table 1, and the resulting fits are shown in
Figure 8a,b, which are the same data values plotted with different x-axes in order to highlight the gate voltage dependence and the temperature dependence of the mobility, respectively.
As can be seen, the TCAD model reproduces the experimental results very nicely for almost all the range of gate voltages (except, perhaps, the subthreshold region) and much of the temperature range (with some deviation at low temperatures for, again, the subthreshold region). However, the model considers none of the typical mechanisms for subthreshold degradation, such as interface traps—though such mechanisms could easily be added—and thus, the discrepancy is perhaps not surprising.
As has also been found in previous modeling attempts in the literature [
10,
11], an impressive fit to experimental mobility results can be obtained without considering the effect of bound band-tail states below the mobility band-edge. However, the existence of such band-tails are still a known reality of a-IGZO films, and thus, we now turn to the question of their inclusion in our microscopic TCAD model.
2.4. Inclusion of Band-Tails and Hypergeometric Functions
As was originally shown by Sir Nevill Mott in a seminal 1967 work, in the presence of heavy disorder—such as in non-crystalline, amorphous, impurity-heavy, or degenerately doped systems—under a fairly generic set of assumptions, the energy states of a physical system will decompose into two distinct energy ranges: a spectrum of non-localized mobile “free” states for energy above a certain cut-off value and a spectrum of non-conducting localized trap-like states for energies below that value [
26,
27]. This specific value that separates these two non-conducting and conducting spectral regions is then called the “mobility edge”. These states, like all trapping states, mean that only some fraction of electrons that comprise the charge density of a disordered system actually contribute to its conductivity.
These so-called “tail states” or “band-tails” of bound states that hang below the mobility band-edge are usually assumed to either have a Gaussian or exponentially decaying energy dependence for their density of states. In a-IGZO, the presence of these bands has been measured [
28,
29,
30,
31] and has been found to be mainly of the exponentially decaying type with a density of states (DOS) given by:
The sum represents the potential to have multiple exponential band-tails with a density-of-states constant of
, starting at an energy
, and having a characteristic decay length of energy of
. The value of
is defined relative to the band-edge
and, thus, spatially varies up and down throughout the system in a manner directly following the band-edge. Thus, to make this clear, we explicitly notate it as position-dependent. Alternatively, one can conceptualize
as the system temperature above which carriers can free themselves from the bound trap states to populate the mobile, delocalized states. Most experiments [
28,
29,
30,
31] show that a maximum of two bands may be present (assuming deep traps are neglected) but that often only a single band is necessary to adequately match experimental results.
Given these results, we specialize our interest here to the case of only a single band-tail of the exponentially decaying type, though the inclusion of Gaussian tails is similar, although the integrals involved cannot be directly analytically solved—as they will be in
Section 2.4.1 for the exponentially decaying case—and require either numerical integration or an analytical approximation, such as in [
32].
However, despite the existence of these band-tails in a-IGZO being beyond dispute, their relative importance in the modeling of the mobility in a-IGZO is somewhat contentious in comparison to that of other amorphous materials such as a-Si, where their effect is known to be crucial. Some theoretical studies have argued that their inclusion is necessary to match experimental results [
5]. However, often such studies treat the effect of the conduction band variation explored in
Section 2.1 using imperfect models or neglect it entirely and, thus, may need to anomalously enhance the effect of these bands to unphysical levels in order to compensate for this important omission. Conversely, the majority of modeling works that demonstrate experimental matching [
4,
10,
11]—including our own in the previous section—only pay lip service to the presence of such tails but then neglect them entirely when actually performing parameter fitting. This strongly suggests that their role is perturbative at best.
There is also a matter of dispute regarding how mobile these “bound” states should be. Many models [
4,
10,
11] assume that electrons in these states are entirely immobile and can only contribute to conduction through the process of so-called “multiple catch and release”, where they thermally excite upwards in energy to mobile states, drift a bit, then are recaptured. However, others allow for so-called variable-range hopping (VRH)—a concept and model also introduced by Mott [
33]—from trap-state to trap-state.
Thus, depending on the modeling objectives, it may be necessary to include multiple tails, with the carrier densities within these tails following their own mobility models. We now demonstrate how such tails can be integrated into the TCAD model considered here.
2.4.1. Analytical Evaluation of the Band-Tail Charge Density
In principle, the inclusion of exponential band-tails in a simulation is straightforward if numerical integration is performed for the relevant integral:
However, in a TCAD simulation, this numerical integration would need to be performed at each and every mesh point, which would dramatically increase the numerical burden. However, for the specific case of exponentially decaying band-tails using the Fermi–Dirac distribution, an analytical evaluation is possible by recasting Equation (
9) as:
where
is a so-called
Gaussian hypergeometric function, with
being the Gamma function.
Although this mathematical equality has been identified and exploited in many places in the literature [
34,
35], to our knowledge, this unexpected and somewhat esoteric equivalence has not actually been proven within the literature and has instead been simply stated without derivation. Thus, for conceptual clarity and convenience, in
Appendix B, we provide a proof of this equality.
The great value of this connection with Gaussian hypergeometric functions is that such functions are well-studied and have many known useful analytical properties. For example, their derivative
is also analytically known, which can be quite useful when, for example, one is using a Newton simulation loop, which benefits from knowing the derivative of the charge density with respect to the potential and/or Fermi level.
However, an even greater benefit for recasting the key integral in terms of a Gaussian hypergeometric function is that a number of common numerical libraries have dedicated, numerically efficient functions for their evaluation. For example, in
Python, there is
scipy.special.hyp2f1, and in
Boost C++, there is
hypergeometric.hpp. In fact, although in the literature there have been a number of
approximations to this integral developed [
35,
36,
37,
38] for use in, for example, compact circuit models, for the purposes of TCAD, we find that any form of approximation is completely unnecessary, as when one uses such pre-existing libraries, the evaluation of these integrals on every mesh point only contributes a tiny fraction to the overall simulation time.
To demonstrate this relative numerical lightweightness, we consider the case of a single trap band with a
of 770 K and a density-of-states constant of
of
(1.6
eV
). We simulate the device from
Section 2.3, but we vary the mesh resolution, thus increasing the number of mesh points,
n, and, thus, the number of band-tail integrals that must be evaluated. We track the computation time required for only the trap-band evaluation portion of the simulation as well as the discrepancy in value between a numerical integration versus an analytical evaluation of the band-tails using
scipy.special.hyp2f1. The results are shown in
Figure 9a. The left y-axis shows the percent of the total simulation time that is spent evaluating the analytical Gaussian hypergeometric functions. Never does this exceed 0.1% of the total simulation time. Furthermore, the right y-axis shows the relative speed-up of analytical evaluation over numerical evaluation, which is about a 650× speed-up for all system sizes. Finally, the inset shows that there is no detectable difference in the computed mobility values whether numerical or analytical approaches are used.
2.4.2. The Effect of Band-Tails on Mobility
Figure 9b shows how the mobility varies as a function of the inverse temperature at a high gate bias of 20 V in the case where there is no trapping band present versus the case where there is one with a density-of-states constant
of
m
(1.6
) and a
of 770 K. Additionally, the difference between a varying band-edge (using the parameters from
Section 2.3) versus a flat non-a-IGZO-like band-edge are shown. It is immediately clear that the presence of a varying band-edge significantly changes the quantitative and qualitative behavior of the temperature dependence of the mobility.
When there is no varying band-edge landscape, the effect of traps is to substantially erode the mobility at low temperatures as more and more carriers freeze into the non-conducting tail-states. However, when varying band-edges are present, the behavior is quite different. The mobility is eroded somewhat since some carriers now lie in the band-tail and thus do not contribute to the conductivity, but this seems to yield a constant shift that is effectively temperature-independent.
This difference in behavior is shown to be even more unusual when looking at
Figure 10a, where the average and maximum fractions of bound electrons (averaged over the full 10 nm film) are shown for the two cases in question. It seems that for the varying-band case, the average fraction of bound carriers actually goes to zero at low temperatures, while the maximum fraction goes to 100%. This is in great contrast to the case of an unvarying band-edge, where the average fraction of bound carriers increases with decreasing temperature.
This is because, whether the band-edge is varying or not, the electrostatic effect of the gate contact forces the system to produce a certain amount of charge. However, when the band-edge is varying in space, there is still only a single Fermi level across the whole system. When the band-edge is not varying, the system is forced to form a surface charge layer that is effectively uniform throughout the device (neglecting the effect of the Schottky contacts), but when it is varying, most of the charge required by electrostatics is given by the “valley” regions—i.e., those regions where the band-edge is lowest—whereas the “hill” regions of high contribute very little.
Thus, the overwhelming majority of the charge required to compensate for the gate electrostatics comes from the comparatively small volume fraction of the films with deep valleys, and in these valleys, most carriers are free, leading to the fraction of bound carriers being very small. This can be seen by considering a 1D cross-section of the film at two points—one where the band-edge is quite low (i.e., a “valley”) and one where it is quite high (a “hill”)—and comparing this to the case where the band-edge is not varying at all. This cross-section is shown in
Figure 10b for two specific points in the plane of the film that were chosen to be near the center (in order to minimize the effect of the contacts) but to have a band-edge of one standard deviation above the average band-edge (i.e.,
, a “typical” hill) at the exact semiconductor–insulator interface and that of a very low valley with a band-edge of two standard deviations below the average band-edge value (i.e.,
, a “deep” valley) at the interface.
As can be seen, the low valley is below the Fermi level and lower than the unvarying case for most of the film thickness, with the Fermi level being deep in the band. Thus, the overwhelming majority of carriers are free carriers.
Capturing effects like this is very difficult in a “compact” model and even in TCAD models that use only an effective uniform mobility model. Thus, we see the great value of using an explicitly microscopic TCAD model when considering trapping-band behavior.