of 8
PHYSICAL REVIEW B
108
, 235207 (2023)
Hot-hole transport and noise phenomena in silicon at cryogenic temperatures from first principles
David S. Catherall
and Austin J. Minnich
*
Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California 91125, USA
(Received 28 June 2023; accepted 21 November 2023; published 19 December 2023)
The transport properties of hot holes in silicon at cryogenic temperatures exhibit several anomalous features,
including the emergence of two distinct saturated drift velocity regimes and a nonmonotonic trend of the current
noise versus electric field at microwave frequencies. Despite prior investigations, these features lack generally
accepted explanations. Here, we examine the microscopic origin of these phenomena by extending a recently
developed
ab initio
theory of high-field transport and noise in semiconductors. We find that the drift velocity
anomaly may be attributed to scattering dominated by acoustic phonon emission, leading to an additional regime
of drift velocity saturation at temperatures
40 K for which the acoustic phonon occupation is negligible; while
the nonmonotonic trend in the current noise arises due to the decrease in momentum relaxation time with electric
field. The former conclusion is consistent with the findings of prior work, but the latter distinctly differs from
previous explanations. This work highlights the use of high-field transport and noise phenomena as sensitive
probes of microscopic charge transport phenomena in semiconductors.
DOI:
10.1103/PhysRevB.108.235207
I. INTRODUCTION
The microscopic processes underlying charge transport
in semiconductors are of fundamental and practical interest
[
1
5
]. Numerical approaches to study transport phenomena
have typically employed Monte Carlo simulations [
1
], which
are capable of treating realistic device geometries [
6
,
7
] and
have been extended to full-band simulators [
8
,
9
]. Within the
last decade, the development of the
ab initio
treatment of
the electron-phonon interaction has enabled the calculation
of transport properties in homogeneous systems without any
fitting parameters [
4
,
5
]. Recent findings with these methods
include the importance of two-phonon scattering in GaAs
[
10
,
11
] and Si [
12
] as well as the discovery of simultaneously
high electron and hole mobilities in BAs [
13
15
]. Most
ab
initio
calculations have been restricted to the low-field regime
in which electrons are in thermal equilibrium with the lattice.
However, outside of the low-field regime, phenomena such as
field-dependent mobility, drift velocity saturation, and nega-
tive differential resistance may occur [
16
18
]. Furthermore,
at high fields qualitatively new information can be obtained
from noise quantities, such as the current fluctuation power
spectral density (PSD), due to the breakdown of equilibrium
relationships including the Einstein [
19
] and Nyquist equa-
tions [
20
]. Recent studies have extended the
ab initio
method
to the high-field regime and applied them to both transport and
noise in GaAs [
11
,
21
] and Si [
12
,
22
].
Charge transport in
p
-Si at cryogenic temperatures and
microwave frequencies is of particular interest due to two
anomalous transport characteristics. First, at large electric
fields of tens of kVcm
1
, semiconductors generally exhibit
drift velocity saturation in which the drift velocity no longer
*
aminnich@caltech.edu
increases with field [
1
]. In
p
-Si below 40 K, an additional
regime of drift velocity saturation occurs at considerably
smaller fields, approximately 0.1 kVcm
1
[
17
]. This fea-
ture has been studied numerically by semianalytical [
23
,
24
]
and full-band Monte Carlo methods [
25
]. While an ini-
tial investigation attributed the feature to the nonparabolic
hole dispersion and Bose-Einstein distribution occupation of
phonons [
23
], a later study instead attributed the behavior
solely to population of the spin-orbit band at high fields [
24
].
The “shoulder” feature has more recently been computation-
ally reproduced using full-band Monte Carlo [
25
], but without
analysis. Therefore, the explanation of the anomalous satura-
tion regime remains unresolved.
A second anomalous feature appears in the PSD. It has
been observed experimentally that in
p
-Si at 10 GHz and 77 K,
a nonmonotonic (peak) trend with increasing dc electric field
is observed for current noise both longitudinal and transverse
to the field. This feature has also been observed in
n
-Si and Ge
[
26
,
27
]. The trend was originally ascribed to an initial increase
due to carrier heating followed by a decrease due to a decreas-
ing momentum relaxation time at high fields [
27
]. In Sec. 9.3
of Ref. [
28
], the trend was attributed to a field-dependent
energy relaxation time which influences the convective noise
mechanism (described in Sec. 7.2 of the same reference). The
convective mechanism arises from the coupling of velocity
and energy fluctuations of the charge carriers, leading to a
decrease in the PSD with increasing electric field for materials
with a sublinear current-voltage characteristic. Given these
differing explanations, the origin of the peak remains unclear.
Here, we investigate these phenomena by computing the
drift velocity and PSD versus field of hot holes in silicon
between 6 and 77 K using a modified implementation of the
high-field
ab initio
method which is applicable at cryogenic
temperatures. We show that the additional drift velocity sat-
uration regime arises from scattering dominated by acoustic
2469-9950/2023/108(23)/235207(8)
235207-1
©2023 American Physical Society
DAVID S. CATHERALL AND AUSTIN J. MINNICH
PHYSICAL REVIEW B
108
, 235207 (2023)
phonon emission. For the PSD, we find that the peak occurs
due to the decrease of the momentum relaxation time with
increasing electric field, an explanation not given in prior
literature. This work demonstrates the use of high-field first-
principles transport and noise calculations as useful tools in
investigating the microscopic phenomena in charge transport.
II. THEORY AND NUMERICAL METHODS
A. Overview
The methods used in this work have been described pre-
viously [
11
,
21
,
22
]. Briefly, the transport of charge carriers
due to an applied electric field may be described by the
Boltzmann transport equation (BTE). For a spatially homo-
geneous system, the BTE is given by
q
E
̄
h
·∇
k
f
λ
=−
λ


λλ


f
λ

,
(1)
where
q
is the carrier charge,
E
is the electric field vector,
f
λ
is
the carrier occupation function indexed by
λ
which represents
the combined indices of band
n
and wave vector
k
, and

f
λ
represents the deviation from equilibrium as

f
λ
=
f
λ
f
0
λ
where
f
0
λ
is the Fermi-Dirac distribution at a specified lattice
temperature and chemical potential.

λλ

is the linearized
collision integral as given by Eq. (3) of Ref. [
21
]. The collision
integral depends on the phonon populations which may be
perturbed by carrier scattering, but owing to the small free-
carrier densities in the relevant experiments (

10
14
cm
3
),
this effect may be neglected here [
17
].
Equation (
1
) can be expressed as a linear system of
equations using the linearized collision integral and a finite
difference approximation for the reciprocal space gradient,
as described in Sec.
II A
of Ref. [
21
]. The only modifica-
tion to the governing equations in the present work is a
change of sign for the reciprocal-space gradient to account
for the charge carriers being holes. The BTE then takes
the form
λ


λλ


f
λ

=
γ
e
E
γ
k
B
T
v
λ,γ
f
0
λ
,
(2)
where E
γ
and
v
λ,γ
are the electric field strength and hole
velocity in the
γ
Cartesian axis,
k
B
is the Boltzmann constant,
and
T
is the lattice temperature. The relaxation operator

λλ

is defined as

λλ

=

λλ

γ
e
E
γ
̄
h
D
λλ

,
(3)
where
D
λλ

is the finite -difference matrix (FDM) approxi-
mating the
γ
-axis component of the reciprocal-space gradient.
The solution to the linear system,

f
λ
, can be used to calcu-
late various observables including the mobility and PSD. The
mobility is given as [
29
]
μ
αβ
(
E
)
=
2
e
2
k
B
T
V
λ
v
λ,α
λ


1
λλ

(
v
λ

f
0
λ

)
,
(4)
where
V
is the supercell volume,
α
is the direction along
which the current is measured, and
β
is the direction
along which the electric field is applied. The PSD is
given as
S
j
α
j
β
(
E
)
=
2
(
2
e
V
)
2
×
R
[
λ
v
λ,α
λ

(
i
ω
I
+

)
1
λλ

(
f
s
λ

(
v
λ

V
β
))
]
,
(5)
where
ω
is the angular frequency,
I
is the identity matrix, and
f
s
λ
=
f
0
λ
+

f
λ
is the steady distribution. Here,
V
β
is the drift
velocity, given by
V
β
=
1
N
λ
v
λ,β
f
s
λ
,
(6)
where
N
=
λ
f
λ
is the number of holes in the Brillouin
zone. Further details are available in Sec.
II B
of Ref. [
21
].
B. Finite-difference matrix
Within the Monkhorst-Pack Brillouin zone discretization
scheme used here [
30
], the reciprocal-space gradient may
be approximated by an FDM as defined by Eqs. (8) and
(B2) of Refs. [
31
,
32
], respectively. In prior high-field works
[
11
,
12
,
21
,
22
], the inclusion of only first-nearest neighbors
provided adequate accuracy at
T

77 K. At temperatures
below 77 K, we found the quality of this numerical finite
difference approximation to be poor. This issue could, in prin-
ciple, be addressed by increasing the grid density. However,
for
T
<
77 K the grid density cannot be made sufficiently
large to achieve the needed accuracy due to computational
limitations.
Instead, we utilize the same prior FDM method but use
multiple shells for the finite difference approximation, as orig-
inally proposed in Ref. [
32
]. To determine the shell weights,
an analytical method cannot be used since the combination of
a body-centered cubic Monkhorst-Pack mesh for Si and the
use of multiple shells leaves low-order mixed partial deriva-
tives which cannot be eliminated using a finite number of
shells. To overcome this limitation, we solve for shell weights
for a selected electronic band by the minimization of an ob-
jective function measuring the difference between exact and
approximate derivatives
Error(
T
,
{
w
b
}
)
=
λ
([
f
0
λ
k
γ
]
[
D
γ
(
{
w
b
}
)
f
0
λ
]
)
2
f
0
λ
λ
[
f
0
λ
k
γ
]
2
f
0
λ
,
(7)
where [
f
0
λ
k
γ
] is the analytical derivative of the Fermi-Dirac
function
f
0
λ
and
D
γ
(
{
w
b
}
) is the FDM representation of a
single
γ
-axis component of the reciprocal-space gradient. The
temperature dependence of the objective function arises from
the temperature dependence of
f
0
λ
. The FDM is defined using
{
w
b
}
, which is the set of shell weights and are the minimiza-
tion parameters, subject to the constraint
b
w
b
=
1. The
subscript
b
identifies the nearest-neighbor shell. We define the
total number of shells in a scheme by
B
.
235207-2
HOT-HOLE TRANSPORT AND NOISE PHENOMENA IN ...
PHYSICAL REVIEW B
108
, 235207 (2023)
TABLE I. Calculated FDM shell weights
w
b
using Eq. (
7
)and
T
=
77 K for various numbers of shells
B
.
Shell (
b
)
Total shells (
B
)1 2 3 4 5 6 7
1
1.000
2
0.725
0.275
3
1.390
0.530
0.920
4
1.270
0.600
0.490
0.380
5
1.280
0.590
0.480
0.370
0.020
6
1.310
0.590
0.420
0.680
0.030
0.170
7
1.380
0.730
0.330
0.980
0.450
0.100
0.550
It was found that relative to the one-shell case (
B
=
1) at
T

77 K, the error defined by Eq. (
7
) decreased by 96% with
the inclusion of three shells (
B
=
3), and by a further 77%
with seven shells (
B
=
7). In terms of transport properties,
we compare the high-field dc mobility in different cases via
RMS difference (defined as
||
y
1
y
2
||
2
/
||
y
1
||
2
, where
y
is the
property of interest). The mobility obtained using
B
=
3 and
B
=
7 agreed with the
B
=
1 case to within 2.1% even at 20 K
and up to 10 kVcm
1
. The RMS difference between
B
=
3
and
B
=
7 was 0.3%. This agreement indicates that, despite
being optimized to accurately calculate the gradient of the
Fermi-Dirac function
f
0
λ
, the higher-order FDM still provides
an accurate result when applied to the hot-carrier distribution
function
f
λ
.
For noise properties, however,
B

3 was found to be
necessary to remove numerical artifacts. Therefore, we used
B
=
3 at 77 K and
B
=
7for
T
<
77 K. For points near the
edge of the defined grid for which the default number of shells
(
B
) is not possible, we used the maximum number of shells
with fully defined neighboring points. For uniformity, the
shell weights were calculated using
T
=
77 K, as temperature
had only a negligible effect on
{
w
b
}
. For example, calculating
{
w
b
}
for
B
=
7at
T
=
20 K results in an RMS difference in
the 20 K drift velocities of only 3.7% compared to calculations
using shell weights optimized at 77 K. Shell weights
w
b
for
each value of
b
and
B
are as shown in Table
I
.
C. Numerical details
Band structure and phonon calculations were performed
in Q
UANTUM
ESPRESSO [
33
]asinRef.[
22
]. Fine-grid in-
terpolation and scattering rate calculations were performed
via P
ERTURBO
[
34
]. A 250
3
(250
×
250
×
250) grid with
2.5 meV Gaussian smearing was used for
T
>
6 K, and
a 325
3
grid with 1.25 meV smearing was used for
T
=
6 K. Increasing the Gaussian width to 1.55 meV at 6 K
resulted in a maximum change in drift velocity of 1.7%
at any field. The energy windows for scattering rate calcu-
lations were 80 and 57 meV for
T
>
6 K and
T
=
6K,
respectively. Increasing the energy window to 90 meV at
77 K and a grid of 180
3
resulted in a 3% difference at
the maximum field studied (10 kVcm
1
), and increasing the
window to 70 meV at 6 K and a grid of 200
3
resulted in
a similar 3.5% difference at 0.84 kVcm
1
. The linearized
Boltzmann transport equation [
21
] was then solved via the
GMRES algorithm [
35
] to calculate the steady-state distribu-
tion function, mobility, and noise quantities [
21
,
22
]. Only the
heavy holes were included to reduce computational costs; it
was found that at 77 K omitting the light hole band resulted
in an RMS difference of only 1.5% compared to results with
both the heavy and light holes, and this error decreases with
decreasing temperature owing to the decreasing occupation of
the light and split-off bands.
III. RESULTS
A. Drift velocity
We first investigate the drift velocity characteristics of
p
-Si. The calculated drift velocity-field curves at various cryo-
genic temperatures are shown in Fig.
1
. We observe that
the
ab initio
calculations predict the two regimes of veloc-
ity saturation that are observed experimentally. The good
agreement when only considering electron-phonon scattering
indicates that ionized impurity scattering is negligible in these
high-purity samples (

100 k
cm) even at cryogenic tem-
peratures, as originally concluded in the experimental studies
FIG. 1. Hole drift velocity versus electric field in the [100] direc-
tion. Experimental values from Ref. [
17
] (6 K as red circles, 20 K
as green upward triangles, 30 K as yellow squares, 45 K as blue
downward triangles). Calculated values also shown (dashed lines
with identical color scheme to experiment).
235207-3
DAVID S. CATHERALL AND AUSTIN J. MINNICH
PHYSICAL REVIEW B
108
, 235207 (2023)
FIG. 2. (a) Calculated hole-phonon scattering rate versus energy
from the valence band maximum at temperatures of 6 K (red circles),
20 K (green upward triangles), 30 K (yellow squares), 45 K (blue
downward triangles), and 77 K (magenta crosses). (b) Hole occu-
pation function versus energy from the valence band maximum at
various electric fields applied in the [100] direction. The occupation
is obtained using a kernel density estimate. Colors correspond to the
same lattice temperatures as in (a). Data at 6 K are shown up to the
57 meV energy window at that temperature.
[
17
,
18
]. At

10 Vcm
1
, the drift velocity increases linearly
with field following the Ohmic trend, followed by an in-
crease with lesser slope from 10–200 Vcm
1
corresponding
to the anomalous saturation regime. At
200 Vcm
1
,the
drift velocity increases further with nearly the same slope
as in the Ohmic regime, but eventually enters another satu-
ration regime at
500 Vcm
1
. Near-quantitative agreement
with experiment is achieved to within 8% (RMS) over all
temperatures.
To gain insight into the origin of the low-field saturation
regime, we show the calculated hole-phonon scattering rates
versus energy in Fig.
2(a)
. In the following discussion, we
focus only on the 20 K case for simplicity. At this temper-
ature, we observe that the scattering rates increase rapidly
with increasing energy; from hole energies of 0 to 20 meV
there is an increase in scattering rates by a factor exceeding
200. However, from 20 to 60 meV the rates are relatively
constant, only rising by a factor of 4. Above 60 meV optical
phonon emission becomes possible and the scattering rates
again rapidly rise with energy.
These differing trends below and above
20 meV arise
due to two factors. First, below
40 K, the population of
zone-edge acoustic phonons, which have the largest density of
states among acoustic phonons, rapidly diminishes with tem-
perature. Consequently, the acoustic phonon absorption rate
approaches zero with decreasing temperature while the acous-
tic phonon emission rate approaches a constant, which causes
a large disparity in scattering rates between low and high
energies where absorption and emission are dominant, respec-
tively. Second, the phonon density of states, which increases
quadratically with energy up to a maximum at
20 meV
[
36
], causes a corresponding increase in scattering rates with
energy from 0–20 meV. These factors together result in an
strongly temperature-dependent increase in scattering rates
from 0 to
20 meV and relatively constant rates at higher
energies up to
60 meV.
We next examine the occupation function at various fields
to establish which hole energies and corresponding scatter-
ing rates are relevant when the saturation feature occurs.
We plot the hole occupation function versus energy and
applied electric field in Fig.
2(b)
.At1Vcm
1
, we ob-
serve that nearly all holes are confined to energies less
than 10 meV, and the hole occupation exhibits a clear de-
pendence on temperature. For sufficiently strong electric
fields (

10 Vcm
1
), the weight of the distribution shifts
towards 10–30 meV as shown in Fig.
2(b)
at 100 Vcm
1
,
where the slope of drift velocity with field reaches a lo-
cal minimum. In this energy range, we observe the rapid
increase in scattering rates with increase in energy due to
zone-edge acoustic phonon emission. At 500 Vcm
1
,shown
also in Fig.
2(b)
, the hole distribution has weight in the
30–60 meV region, where scattering rates are relatively
constant. Above 500 Vcm
1
, some holes have energy ex-
ceeding
60 meV and may undergo optical phonon emission
processes.
We may now explain the relationship between acoustic
phonon scattering and the first regime of drift velocity satu-
ration. At low fields (

10 Vcm
1
), the drift velocity is linear
with field corresponding to the Ohmic regime, as holes are
approximately in equilibrium with the lattice and the mean
scattering rate remains independent of field, a trend which
occurs at all temperatures. However, at cryogenic tempera-
tures and high fields (
40–200 Vcm
1
), carrier heating shifts
the occupation to energies where the scattering rates are 1–2
orders of magnitude larger than those relevant at low fields.
The higher scattering rates at these energies more than com-
pensates the higher band velocity, causing a decrease in
mobility and the drift velocity saturation.
At higher fields, above
200 Vcm
1
, the increasing band
velocity with energy compensates the relatively smaller in-
crease in scattering rates. Therefore, the saturation “shoulder”
ends, and drift velocity again increases with increasing field.
Finally, at about
500 Vcm
1
, optical phonon emission be-
gins and another regime of drift velocity saturation occurs.
At higher temperatures

77 K, the scattering rates below
20 meV do not show a pronounced increase, and therefore
the anomalous saturation regime does not occur. From these
235207-4
HOT-HOLE TRANSPORT AND NOISE PHENOMENA IN ...
PHYSICAL REVIEW B
108
, 235207 (2023)
FIG. 3. PSD of current noise versus electric field applied in
the [110] direction at 10 GHz. Experimental values from Ref. [
27
]
(measurements taken parallel to [110] as circles, measurements taken
in the transverse direction [1
10] as upward triangles). Calculated
data also shown [parallel (

) as solid line, perpendicular (
)as
dashed line].
findings, we may conclude that it is the energy-dependence of
the acoustic phonon scattering rates at cryogenic temperatures
which cause the saturation “shoulder” in
p
-Si between 10 and
200 Vcm
1
for
T

40 K, in agreement with a prior work
[
23
].
B. Power spectral density
We now turn to the nonmonotonic trend of the mi-
crowave PSD versus electric field. This feature appears in the
experimental 10 GHz PSD, which first increases and then
decreases with an increasing electric field, leading to a peak
at around 20 Vcm
1
[
26
,
27
]. We have calculated the 10 GHz
PSD for longitudinal and transverse measurement directions
([110] and [1
10], respectively) versus [110] electric field
strength, shown in Fig.
3
. Normalized values are shown to
facilitate the comparison of trends. The absolute PSD values
would be uniformly overestimated by around 25% based on
previous calculations of the mobility and diffusion coefficient
and considering the Nyquist and Einstein relations linking the
PSD, diffusion coefficient, and mobility at equilibrium [
22
].
The calculations do not agree quantitatively with experiment
but show a qualitatively similar nonmonotonic trend. A firm
conclusion regarding the origin of the discrepancy is difficult
to draw owing to discrepancies between experimental reports
(
c.f. Fig. 1(b) of [
26
] and the various other experimental re-
ports of the diffusion coefficient, namely Fig. 2(a) of Ref. [
37
]
and Fig. 3 of Ref. [
38
]) and the fact that only one reference ex-
ists for this particular data set. In the following discussion, we
therefore consider the qualitative experimental trends rather
than focusing on quantitative agreement.
Previous studies have attributed the trend to the com-
petition of carrier heating (increasing PSD) and decreasing
relaxation time (decreasing PSD) [
27
], as well as the
convective noise mechanism which is influenced by the
field-dependent energy relaxation time [
28
]. However, each
explanation has inconsistencies. Regarding the first expla-
nation, we note that the 10 GHz PSD measurements by
Ref. [
26
,
27
] were used to calculate the diffusion coefficient,
which is valid only in the low-frequency limit [
39
](
ωτ
m

1,
where
τ
m
is the momentum relaxation time appearing in the
Drude conductivity expression [
40
]). In this limit, an increase
in transport properties with increasing electric field due to
carrier heating would be expected in not just the PSD but also
in dc quantities such as mobility. However, the hole mobility
in Si is observed to decrease with increasing field. Therefore,
this explanation is inconsistent with the observed trends of
other transport properties. Regarding the second explanation,
the convective noise mechanism only affects the longitudinal
PSD, not the transverse PSD (see Sec. 7.2 of Ref. [
28
]). Be-
cause the PSD peak appears in both measurement directions,
the convective mechanism cannot be responsible for the peak.
Therefore, the underlying cause of the nonmonotonic PSD
feature remains unclear.
To identify the origin of the peak, we calculated the PSD
versus both electric field and frequency. The result is shown
in Fig.
4(a)
. The longitudinal PSD calculation in Fig.
3
cor-
responds to a horizontal slice of this plot at 10 GHz. The
location of the peak, given by the dotted line, is observed to
depend on both electric field and frequency. Because the peak
is observed in our calculations which include only the heavy-
hole band, interband scattering may be ruled out. Additionally,
the feature is also observed when optical phonons are omitted,
eliminating carrier streaming noise as the origin (see Sec. 7.4
of Ref. [
28
]). Therefore, we conclude that the peak must arise
from qualities of the hole and phonon dispersion and acoustic
phonon scattering rates.
Knowing that the convective mechanism cannot be respon-
sible for the PSD peak, and without any other mechanism
by which energy relaxation may produce the feature, we
take as a working hypothesis that the energy dependence of
the momentum relaxation time leads to the trend. However,
within the fully
ab initio
framework, testing this hypothesis is
difficult. Therefore, we introduce a simplified model relating
the PSD to
τ
m
(
E
). Assuming equipartition and neglecting
anisotropy, the PSD takes the form [
41
]
S
j
(
ω,
E
)
=
S
j
(0
,
E
)
1
+
(
ωτ
)
2
,
(8)
where
τ
represents a characteristic relaxation time of the
system (for instance momentum, energy, or intervalley re-
laxation), and
E
represents the electric field strength. A
derivation of this expression, considering the constant re-
laxation time to be
τ
m
, can be found in Secs. 3 and 7
of Ref. [
28
]. A simple although less rigorous derivation
can be performed by combining the Nyquist current noise
(
Refs. [
42
,
43
], Eqs. (
6
) and (
1
), respectively) and Drude ac
conductivity
(
Ref. [
40
], Eq. (1.29)) to yield Eq. (
8
), also
neglecting the field dependence.
For our model, we assume a monoenergetic hole distribu-
tion characterized by a momentum relaxation time depending
on electric field
τ
m
(
E
). For simplicity, we further approxi-
mate
S
j
(0
,
E
)
τ
m
(
E
) as the low-frequency PSD is directly
235207-5
DAVID S. CATHERALL AND AUSTIN J. MINNICH
PHYSICAL REVIEW B
108
, 235207 (2023)
FIG. 4. (a) Contour map of normalized
ab initio
PSD of current fluctuations versus frequency and electric field. Location of peak in the
horizontal slices (PSD vs field curves) also shown (dotted black line). (b) Contour map of normalized model PSD [Eq. (
9
)] versus frequency
and electric field. An analogous peak feature is observed as the
ab initio
PSD. Location of peak in the horizontal slices shown (dotted black
line), as well as
ωτ
m
(
E
)
=
1 (solid black line). (c)
Ab initio
PSD vs electric field in (a) (solid red line) and modeled PSD in (b) (solid blue
line) at 60 GHz. The normalized momentum relaxation time is shown (dotted green line) as well as its inverse (dashed magenta line). The
normalization constant for
τ
m
(
E
) is 3.86 ps.
proportional to the momentum relaxation time through the
fluctuation-diffusion [
39
] and Einstein [
19
] relationships.
These assumptions yield
S
j
(
ω,
E
)
τ
m
(
E
)
1
+
(
ωτ
m
(
E
)
)
2
.
(9)
Now, the Einstein relationship holds only in the low-field
limit, but it remains approximately valid at fields where the
carrier distribution is sufficiently close to the equilibrium one.
In the present case, the mean carrier energy at 20 Vcm
1
,
where nonmonotonic behavior is apparent, is only 2.2% larger
than the equilibrium value, so the error arising from the use
of the Einstein relationship should be on the order of this
difference (

10%). Therefore, it is expected that Eq. (
9
)
provides a qualitatively accurate description of PSD behavior
at the electric fields of interest.
To plot Eq. (
9
), an expression for
τ
m
(
E
) is needed. From
our
ab initio
results, we calculated the momentum relaxation
time at each electric field by fitting a Lorentzian to the fre-
quency dependence of the transverse ac mobility with an
identical functional form as in Eq. (
9
). The resulting function
τ
m
(
E
) is shown in normalized form in Fig.
4(c)
.
Using this function, we compute PSD versus frequency
and electric field strength using Eq. (
9
). The result is given in
Fig.
4(b)
. We observe a qualitatively similar “ridge” trend as
seen in the
ab initio
PSD. Horizontal slices, which yield PSD
versus field strength, taken at high frequencies, display a peak,
the location of which follows the same trend as in Fig.
4(a)
.
Within this model, the location of the peak corresponds to
ωτ
m
(
E
)
=
1, providing evidence that the field dependence of
the momentum relaxation time
τ
m
(
E
) is responsible for the
nonmonotonic PSD trend.
This behavior can be understood by inspecting the rela-
tive magnitude of the numerator
τ
m
(
E
) and the denominator
1
+
(
ωτ
m
(
E
))
2
, which represent carrier scattering and carrier
inertia [
44
], respectively. At sufficiently high temperatures
and low frequencies such that
ωτ
m
(
E
)

1, the PSD is
proportional to
τ
m
(
E
) and decreases with field due to
τ
m
(
E
)
decreasing. However, when
ωτ
m
(
E
)
1 is satisfied, which
occurs at some combinations of sufficiently low temperatures,
high frequencies, and low fields, then, with increasing field
as
τ
m
(
E
) decreases the inertial term 1
+
(
ωτ
m
(
E
))
2
decreases
faster than
τ
m
(
E
) itself. Thus, the PSD rises as PSD(
ω,
E
)
1
m
(
E
). Once
ωτ
m
(
E
)

1 is satisfied again at high enough
fields, which will occur regardless of temperature and fre-
quency, the PSD again decreases with field. These different
behaviors are displayed in Fig.
4(c)
, along with horizontal
slices of Figs.
4(a)
and
4(b)
at 60 GHz. We have chosen this
higher frequency for the plots to make the relevant features
in the model more visible. We observe that at low fields, both
the modeled and
ab initio
PSD exhibit the expected
τ
m
(
E
)
1
behavior, and the
τ
m
(
E
) behavior at high fields. These differ-
ent trends enable a peak in the PSD, driven solely by the field
dependence of the momentum relaxation time.
We now discuss why the
ab initio
PSD results display a
peak even when
ωτ
m
(
E
=
0)

1 is not satisfied, a condition
strictly necessary for a peak to appear in the model. In a real
system, the hole distribution is not characterized by a single
energy, as was assumed in the model. Therefore, although the
momentum relaxation time
τ
m
(
E
) may satisfy
ωτ
m
(
E
)
1,
the lifetime of some states
τ
λ
=

1
λ,λ
may span a range of
values of
ωτ
λ
. Therefore, even at frequencies where
ωτ
m
(
E
=
0)
<
1, some electronic states for which
ωτ
λ
1 may have
non-negligible occupation. The carriers in these states may
then contribute to the occurrence of a nonmonotonic PSD
trend, even at lower frequencies than are required to satisfy
ωτ
m
(
E
=
0)

1.
IV. DISCUSSION
We now examine our findings in the context of prior works.
The role of the field dependence of the momentum relax-
ation time in producing nonmonotonic trends of ac mobility
with electric field has been previously reported [
45
]. Various
235207-6
HOT-HOLE TRANSPORT AND NOISE PHENOMENA IN ...
PHYSICAL REVIEW B
108
, 235207 (2023)
studies have incorporated field dependence into the momen-
tum relaxation time of ac mobility models and observed
different trends depending on the frequency. At low fre-
quencies such that
ωτ
m
(
E
=
0)

1, the ac mobility is both
theoretically predicted [
46
,
47
] and experimentally observed
[
48
] to decrease monotonically with increasing electric field,
as expected for systems with sublinear current-voltage charac-
teristics. However, at frequencies such that
ωτ
m
(
E
=
0)

1,
other studies predicted a nonmonotonic dependence of ac
mobility on electric field [
49
,
50
].
Experimentally, nonmonotonic trends of ac mobility with
electric field have been reported for
n
-InSb [
45
,
51
]. These ex-
periments measured ac mobility at 77 K as a function of field
from 0–165 Vcm
1
and from 0–83
.
7 GHz, which includes the
case
ωτ
m
(
E
=
0)
1. The measurements show nonmono-
tonic behavior when
ωτ
m
(
E
=
0)
>
1, with an initial rise
in the ac mobility at low fields (

50 Vcm
1
), in contrast
to the decreasing trend at lower frequencies. Furthermore,
the data agree qualitatively with a hydrodynamic model [
51
]
which predicts nonmonotonicity in the high-field ac mobility
for frequencies above
50 GHz. Since the ac mobility is
proportional to the PSD at low electric fields through the
Nyquist relationship, these prior findings are consistent with
our conclusions. On the other hand, prior studies of the mi-
crowave PSD concluded the nonmonotonic trend was due to a
competition between carrier heating and decreasing relaxation
time with increasing field [
27
], and convective noise [
28
].
However, our study and prior ac mobility studies support that
the field dependence of the momentum relaxation time alone
as the origin.
V. SUMMARY
We have investigated the origin of anomalous high-field
transport and noise characteristics in
p
-Si at cryogenic tem-
peratures using a modified high-field formalism for
ab initio
charge transport. We find that the additional regime of drift
velocity saturation that occurs for
T
<
40 K is due to acoustic
phonon emission, in agreement with a prior work. We also
find that the peak in the power spectral density of current
fluctuations which occurs at 77 K and 10 GHz is due to the
field dependence of the momentum relaxation time, contrary
to the conclusions of prior studies. This work highlights the
capabilities of the
ab initio
method for providing microscopic
insight into high-field transport and noise phenomena in semi-
conductors.
ACKNOWLEDGMENT
This work was supported by the National Science Founda-
tion under Award No. 1911926.
[1] M. Lundstrom,
Fundamentals of Carrier Transport
, 2nd ed.
(Cambridge University Press, Cambridge, UK, 2000).
[2] D. Ferry,
Semiconductor Transport
(CRC Press, Boca Raton,
FL, 2000).
[3] S. Poncé, W. Li, S. Reichardt, and F. Giustino,
Rep. Prog. Phys.
83
, 036501 (2020)
.
[4] M. Bernardi,
Eur.Phys.J.B
89
, 239 (2016)
.
[5] F. Giustino,
Rev. Mod. Phys.
89
, 015003 (2017)
.
[6] C. Jacoboni and P. Lugli,
The Monte Carlo Method for Semicon-
ductor Device Simulation
, 1st ed. (Springer, New York, 1989).
[7] J. Mateos, T. Gonzalez, D. Pardo, V. Hoel, and A. Cappy,
IEEE
Trans. Electron Devices
47
, 1950 (2000)
.
[8]EditedbyK.Hess,
Monte Carlo Device Simulation: Full Band
and Beyond
, The Springer International Series in Engineering
and Computer Science Vol. 144 (Springer Science + Business
Media, LLC, 1991).
[9] C. Jungemann, S. Keith, M. Bartels, and B. Meinerzhagen,
IEICE Transactions on Electronics
E82-C
, 870 (1999).
[10] N.-E. Lee, J.-J. Zhou, H.-Y. Chen, and M. Bernardi,
Nat.
Commun.
11
, 1607 (2020)
.
[11] P. S. Cheng, J. Sun, S.-N. Sun, A. Y. Choi, and A. J. Minnich,
Phys.Rev.B
106
, 245201 (2022)
.
[12] B. Hatanpää, A. Y. Choi, P. S. Cheng, and A. J. Minnich,
Phys.
Rev. B
107
, L041110 (2023)
.
[13] T.-H. Liu, B. Song, L. Meroueh, Z. Ding, Q. Song, J. Zhou, M.
Li, and G. Chen,
Phys.Rev.B
98
, 081203(R) (2018)
.
[14] J. Shin, G. A. Gamage, Z. Ding, K. Chen, F. Tian, X. Qian, J.
Zhou, H. Lee, J. Zhou, L. Shi, T. Nguyen, F. Han, M. Li, D.
Broido, A. Schmidt, Z. Ren, and G. Chen,
Science
377
, 437
(2022)
.
[15] S. Yue, F. Tian, X. Sui, M. Mohebinia, X. Wu, T. Tong, Z. Wang,
B. Wu, Q. Zhang, Z. Ren, J. Bao, and X. Liu,
Science
377
, 433
(2022)
.
[16] C. Jacoboni and L. Reggiani,
Adv. Phys.
28
, 493 (1979)
.
[17] G. Ottaviani, L. Reggiani, C. Canali, F. Nava, and A.
Alberigi-Quaranta,
Phys.Rev.B
12
, 3318 (1975)
.
[18] C. Canali, C. Jacoboni, F. Nava, G. Ottaviani, and A.
Alberigi-Quaranta,
Phys.Rev.B
12
, 2265 (1975)
.
[19] L. Reggiani, General theory, in
Hot-Electron Transport in
Semiconductors
, edited by L. Reggiani (Springer, Berlin,
Heidelberg, 1985), pp. 7–86.
[20] J. P. Nougier, in
Physics of Nonlinear Transport in Semicon-
ductors
, NATO Advanced Study Institutes Series: Series B,
Physics,Vol.52,editedbyD.K.Ferry,J.R.Barker,andC.
Jacoboni (Plenum Press, New York, 1980), pp. 415–465.
[21] A. Y. Choi, P. S. Cheng, B. Hatanpää, and A. J. Minnich,
Phys.
Rev. Mater.
5
, 044603 (2021)
.
[22] D. S. Catherall and A. J. Minnich,
Phys.Rev.B
107
, 035201
(2023)
.
[23] J. Von Borzeszkowski,
Physica Status Solidi B
73
, 607
(1976)
.
[24] Y. Ohno,
J. Appl. Phys.
64
, 4549 (1988)
.
[25] B. Fischer and K. R. Hofmann,
Appl. Phys. Lett.
76
, 583 (2000)
.
[26] V. Bareikis, R. Barkauskas, A. Galdikas, and R. Katilius, Sov.
Phys. Semicond.
14
, 1046 (1980).
[27] V. Bareikis, A. Galdikas, and R. Miliusyte, in
Hot Electron
Diffusion / Diffuziia goriachikh elektronov
, Electrons in Semi-
conductors / Elektrony v poluprovodnikakh Vol. 3 (Vilnius:
Mokslas, 1981), Chap. 4, pp. 127–161, academy of Sciences
of the Lithuanian SSR / Lietuvos TSR Moksly Akademija.
235207-7
DAVID S. CATHERALL AND AUSTIN J. MINNICH
PHYSICAL REVIEW B
108
, 235207 (2023)
[28] H. L. Hartnagel, R. Katilius, and A. Matulionis,
Microwave
Noise in Semiconductor Devices
(John Wiley & Sons, Inc., New
York, 2001).
[29] W. Li,
Phys. Rev. B
92
, 075405 (2015)
.
[30] H. J. Monkhorst and J. D. Pack,
Phys. Rev. B
13
, 5188
(1976)
.
[31] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt,
and N. Marzari,
Comput. Phys. Commun.
178
, 685 (2008)
.
[32] N. Marzari and D. Vanderbilt,
Phys.Rev.B
56
, 12847 (1997)
.
[33] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C.
Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I.
Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi,
R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M.
Lazzeri, L. Martin-Samos
et al.
,
J. Phys.: Condens. Matter
21
,
395502 (2009)
.
[34] J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M.
Bernardi,
Comput. Phys. Commun.
264
, 107970 (2021)
.
[35] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T.
Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser,
J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman,
N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J.
Carey
et al.
,
Nature Methods
17
, 261 (2020)
.
[36] D. S. Kim, H. L. Smith, J. L. Niedziela, C. W. Li,
D. L. Abernathy, and B. Fultz,
Phys. Rev. B
91
, 014307
(2015)
.
[37] L. Reggiani, J. Phys. Soc. Jpn. Suppl. A
49
, 317 (1980), proc.
15th Int. Conf. Physics of Semiconductors.
[38] D. Gasquet, B. Azais, J. C. Vaissiere, and J. P. Nougier, in
Noise
in Physical Systems and 1/f Noise - 1985
,editedbyA.d’Amico
and P. Mazzetti (Elsevier, Amsterdam, Netherlands, 1986),
pp. 231–233.
[39] S. V. Gantsevich, V. L. Gurevich, and R. Katilius,
Riv. del
Nuovo Cim.
2
, 1 (1979)
.
[40] N. W. Ashcroft and N. D. Mermin,
Solid State Phyiscs
,1sted.
(Saunders College Publishing, Philadelphia, PA, 1976).
[41] Y. K. Pozhela, Transport parameters from microwave conduc-
tivity and noise measurements, in
Hot-Electron Transport in
Semiconductors
, edited by L. Reggiani (Springer, Berlin, Hei-
delberg, 1985), pp. 113–147.
[42] H. Nyquist,
Phys. Rev.
32
, 110 (1928)
.
[43] J. B. Johnson,
Phys. Rev.
32
, 97 (1928)
.
[44] K. Champlin, D. Armstrong, and P. Gunderson,
Proc. IEEE
52
,
677 (1964)
.
[45] E. Bonek,
J. Appl. Phys.
43
, 5101 (1972)
.
[46] A. Gibson, J. Granville, and E. Paige,
J. Phys. Chem. Solids
19
,
198 (1961)
.
[47] B. V. Paranjape,
Phys. Rev.
122
, 1372 (1961)
.
[48] G. H. Glover,
J. Appl. Phys.
44
, 1295 (1973)
.
[49] D. Mukhopadhyay and B. Nag,
Electron. Lett.
5
, 20 (1969)
.
[50] K. Seeger and H. Pötzl, in
The Boltzmann Equation
, edited
by E. G. D. Cohen and W. Thirring (Springer Vienna, Vienna,
1973), pp. 341–378.
[51] E. Bonek, H. Pötzl, and K. Richter,
J. Phys. Chem. Solids
31
,
1151 (1970)
.
235207-8