Silicon–germanium heterojunction bipolar transistors (HBTs) are of interest as low-noise microwave amplifiers due to their competitive noise performance and low cost relative to III–V devices. The fundamental noise performance limits of HBTs are thus of interest, and several studies report that quasiballistic electron transport across the base is a mechanism leading to cryogenic non-ideal IV characteristics that affect these limits. However, this conclusion has not been rigorously tested against theoretical predictions because prior studies modeled electron transport with empirical approaches or approximate solutions of the Boltzmann equation. Here, we study non-diffusive transport in narrow-base SiGe HBTs using an exact, semi-analytic solution of the Boltzmann equation based on an asymptotic expansion approach. We find that the computed transport characteristics are inconsistent with experiments, implying that quasiballistic electron transport is unlikely to be the origin of cryogenic non-ideal IV characteristics. Our work helps to identify the mechanisms governing the lower limits of the microwave noise figure of cryogenic HBT amplifiers.

## I. INTRODUCTION

Silicon–germanium heterojunction bipolar transistors (HBTs) are widely used in microwave applications such as radar and communication systems^{1,2} and show potential in space science^{3} and imaging applications^{4–6} owing to their competitive microwave performance combined with ease of integration with the CMOS process, high yield, and low cost relative to III–V technologies.^{7,8} SiGe HBTs are now approaching the THz domain due to fabrication process advancements that enable the continued scaling of key parameters such as base width and resistance.^{2,9} Moreover, at cryogenic temperatures, they have demonstrated DC gains exceeding 45 dB, noise temperatures $\u22725$ K at 3 GHz with bandwidths of 2–3 GHz,^{10,11} lower $1/f$ noise relative to III–V devices,^{12} and power consumption in the hundreds of microwatts.^{13} These competitive figures of merit make SiGe HBTs of interest for low-noise amplifiers and cryo-CMOS circuits^{14,15} in radio astronomy and quantum computing applications^{13,16} and in oscillators with low phase noise.^{12,17,18}

Following the first reports of SiGe HBTs grown by molecular beam epitaxy,^{19,20} studies of the cryogenic performance of HBTs were performed in the early 1990s.^{19,21} Subsequent work focused on understanding and optimizing base doping and Ge grading profiles for cryogenic performance.^{22} While these works demonstrated the enhanced collector current and transconductance expected at cryogenic temperatures, later studies reported marked deviations from the DC current–voltage characteristics predicted from drift–diffusion theory.^{23,24} Specifically, it was observed that below $\u224877$ K, the collector $(IC)$ and base current $(IB)$ exceeded the predicted values at a given temperature and bias, and were independent of temperature, with the transconductance $gm$ plateauing instead of increasing as $T\u22121$. These non-ideal $IV$ characteristics limit the gain, cutoff frequency, and ultimately the microwave noise figure.

The origin of this behavior has been attributed to several mechanisms, including trap-assisted tunneling at base–emitter voltages ($VBE$) well below the built-in potential;^{23} and non-equilibrium carrier transport, in which electrons quasiballistically traverse the base, at biases comparable to the built-in potential.^{24} The latter effect has been phenomenologically described using an elevated electron temperature that is taken as the effective temperature for thermionic emission at the base–emitter junction.^{10,11,24} As base widths were scaled further down to $\u224820$ nm in recent years, direct tunneling of electrons was reported to contribute to the collector current at biases approaching the built-in potential.^{25–27}

There are several discrepancies with these explanations, however. First, non-ideal collector currents have been reported in first-generation devices with base widths on the order of $100$ nm,^{27} for which direct tunneling is unlikely. Furthermore, at base doping levels above $\u22483\xd71018$ cm$\u22123$ common in modern devices,^{25,27} ionized impurity scattering is expected to dominate. Therefore, quasiballistic transport is not expected to be markedly more pronounced at cryogenic temperatures relative to room temperature as evidenced by the weak temperature dependence of the minority carrier mobility.^{11,28–30} Finally, while electrons may be heated as they traverse the base by the built-in field, the collector current and transconductance depend on the electron temperature at the base–emitter junction prior to transport across the base and thus cannot be affected by the built-in field.

A quantitative, non-empirical description of the quasiballistic transport process would allow a more thorough test of whether quasiballistic transport is a possible origin of the cryogenic DC non-idealities. Steady particle transport across a slab of thickness comparable to the MFPs of the particles is described by the Boltzmann equation.^{31,32} The slab problem has been extensively investigated for radiative,^{32–34} neutron,^{35} phonon,^{32,36} and electron transport^{37–40} as well as for rarefied gases.^{41–43} Although analytical solutions under the “one-speed” or “gray” approximations are possible,^{34,37} solutions considering energy-dependent relaxation times are less tractable analytically. Alternative numerical approaches are computationally expensive owing to the need to discretize both the spatial and reciprocal space coordinates.^{37,38,44} Several works have reported asymptotic series expansion methods, based on the original theory by Grad,^{45} for rarefied gas dynamics^{46} and phonon transport.^{47} These methods could enable the efficient solution of the present problem of electronic transport in a narrow base, but they have yet to be adapted for electronic transport.

Here, we report a study of quasiballistic electron transport across a narrow base using an exact, semi-analytic asymptotic expansion approach to solve the Boltzmann equation. We show that quasiballistic transport leads to decreased collector current for a given base–emitter voltage and does not alter the transconductance. Both of these findings contradict experimental observations, indicating that quasiballistic transport is unlikely to be the origin of cryogenic non-idealities in SiGe HBTs. Alternatively, we hypothesize that the origin of these non-idealities could arise from local inhomogeneities in built-in potential, which are readily observed in Schottky diodes^{48,49} and have been postulated to exist in SiGe HBTs.^{50} Our work advances the understanding of cryogenic transport processes in HBTs that limit gain, transconductance, and thus ultimately microwave noise performance.

## II. THEORY

We begin by describing the semi-analytical asymptotic approach used to solve the one-dimensional Boltzmann equation describing electron transport across the base. Focusing on the DC characteristics, we assume steady transport across a base of width $L$ with prescribed forward and reverse going electron distribution functions at the left and right boundaries, respectively. Then, the Boltzmann equation for the electronic distribution function $f\lambda (x)$ is given by

where $C\lambda \lambda \u2032$ is the collision matrix, $vx\lambda $ is the group velocity for electronic state $\lambda $, $M\lambda \lambda \u2032$ is a dimensionless matrix defined as above, $x$ is the spatial coordinate, $x\xaf\u2261x/L$, and $f\lambda $ is the desired distribution function normalized so that $V\u22121\u2211\lambda f\lambda (x)=n(x)$, where $V$ is a normalizing volume. Macroscopic quantities such as electric current are computed using the appropriate Brillouin zone sum.^{32}

In this work, we will solve this equation using an exact, semi-analytic asymptotic expansion method originally reported for rarefied gases^{46} and recently adapted for phonons.^{47} The theory for phonons is given in detail in Ref. 47 for an isotropic crystal under the relaxation time approximation. Here, we generalize this theory to incorporate the full electron collision matrix and arbitrary crystals and extend its applicability beyond linear deviations from a reference distribution.

In the asymptotic expansion approach, $f\lambda $ and associated observables are written as series expansions with the average Knudsen number $\u03f5\u2261\u27e8\Lambda \u27e9/L$ as the expansion parameter, where the average mean free path $\u27e8\Lambda \u27e9\u2261\u2211\lambda kvx\lambda \tau \lambda f\lambda eq/neq$. Here, $f\lambda eq$ and $neq$ are the reference equilibrium distribution function and electron density at zero base–emitter voltage, and $\tau \lambda $ is the wave vector dependent relaxation time. As an example, $f\lambda $ is written as

where $f\lambda ,l$ is the distribution function for order $l$. Substituting this expansion into Eq. (1), we get

Since $M\lambda \lambda \u2032$ scales with $\u03f5\u22121$, the left and right side of Eq. (3) are offset by one order of $\u03f5$. Equating terms with equal order of *ɛ* for $l=0$, we find that the zeroth-order solution satisfies

where we have split the solution into a pure $x$-dependent term $g0(x)$ and wave vector dependent term $e0,\lambda $; i.e., $f\lambda ,0(x)=e0,\lambda g0(x)$. We note here that the null eigenvector $e0,\lambda $ can be any Boltzmann distribution function, but for this work, we choose $e0,\lambda =f\lambda eq$ corresponding to the Boltzmann distribution for the case where no bias is applied.

We now consider the first-order terms. From Eq. (3), the solution can be written using the zeroth-order solution,

This solution is also determined only up to an $x$-dependent multiple of the null eigenvector, denoted $g1(x)$, which will be obtained after deriving the governing equation for $x$-dependence. The higher-order solutions proceed similarly.

The $x$-dependence of the zeroth-order solution can be derived by enforcing current continuity for the electric current density $Jx$ in the $x$-direction for a 1D steady slab,

which must be satisfied at each order of the expansion. Substituting from Eq. (5) for the first-order term, we find that

The second term in the parentheses vanishes since $vx\lambda $ is odd in $\lambda $, whereas $e0,\lambda $ is even. Therefore, we find

Applying $V\u22121\u2211\lambda f\lambda ,0(x)=n0(x)$ to the above equation, we see that the zeroth-order electron density $n0(x)$ satisfies the diffusion equation. Similarly, it can be shown that higher-order terms of $f\lambda ,l(x)$ also satisfy the diffusion equation (see Appendix A of Ref. 47).

To solve these equations, the boundary conditions at each order must be specified. The boundary conditions at zeroth order are simply the prescribed isotropic, hemispherical electron distribution functions at the edges of the slab specified by $nL$ and $nR$. Note that in general, the actual carrier density at the edges of the slab will differ from these values owing to quasiballistic transport. Thus, the solution at zeroth order is just the diffusion equation solution. Note that the boundary conditions of the original problem are completely satisfied at zeroth order.

We now discuss the boundary conditions at higher orders. Because the zeroth-order solution completely satisfies the boundary conditions of the original problem, to enforce the boundary conditions at higher orders, we must introduce a boundary correction term $\Phi \lambda (x)$ that satisfies the Boltzmann equation in the boundary region. This boundary correction must exactly cancel the contribution from $f\lambda (x)$ at the boundaries and vanish in the bulk.

The Boltzmann equation for the boundary solution at first-order is

where we introduce the stretched boundary-region coordinate $\eta \u2261x/\u03f5$. Focusing on the left boundary at $x\xaf=0$, the condition enforced on this first-order boundary term is that it must cancel the first-order bulk term: $\Phi \lambda ,1|0=\u2212f\lambda ,1|0$. Using Eq. (5), the condition for the left boundary is

where, as in Ref. 47, we anticipate the boundary term to scale as the gradient of the previous order solution,

An analogous boundary condition applies to the right boundary. The unknown constant $c1$ captures the jump-type boundary condition for non-diffusive corrections at first order. Solving for $c1$ involves obtaining the eigenvalues and eigenvectors of $M$ and is described in detail in Ref. 47. In brief, the general solution to Eq. (9) can be written as a linear combination of the eigenvectors corresponding to the negative eigenvalues of $M$,

where $Ai$ are unknown coefficients and $\rho i$ and $h\lambda ,i$ are the negative eigenvalues and the corresponding eigenvectors, respectively. Only the negative eigenvalues are used so that the boundary solution tends to zero in the bulk. Then, Eq. (10) gives a linear system of equations for $Ai$ and $c1$. We note that in the formulation in Ref. 47, $c1$ represents a deviation in a linearized quantity, such as temperature, relative to an equilibrium value. However, in the present formulation, $c1$ multiplies the absolute distribution function rather than a deviation and is not restricted to linear deviations from that distribution.

After calculating $c1$, the first-order solution is completely specified using the diffusion equation for $g1(x)$ and the boundary condition from Eq. (11). The analysis progresses similarly for higher-order solutions. For the specific slab problem here, it can be shown that the jump coefficients associated with second-order derivatives vanish, as in Appendix C of Ref. 47. Therefore, $c1$ is the only required coefficient, allowing the asymptotic expansion to be summed over all higher-order terms. After summing $f\lambda $ over the Brillouin zone to obtain the carrier density, we obtain

where $n0(x)$ is the zeroth-order carrier density and the higher-order terms correct for non-diffusive transport. Because the macroscopic constitutive relation between the electric current density and the number density gradient applies in the bulk,^{47} the current density including the effect of quasiballistic transport is given by

Here, $De,Si$ is the bulk diffusivity. These results have been extensively validated by comparison to Monte Carlo simulations in Ref. 47. We also observe from Eq. (13) that the dependence of carrier density and thus electric current on base–emitter voltage $VBE$ is only through the left boundary condition $nL\u2248eqVBE/kT$, allowing the current–voltage characteristics for arbitrary $VBE$ to be obtained once $c1$ is computed.

We briefly discuss the features of the present asymptotic approach compared to those reported previously for electronic transport. An early study of the slab electron transport problem employed a direct scattering matrix solution to the Boltzmann equation^{38} that required discretization in the wave vector and spatial coordinates and was thus numerically expensive. Approximate methods such as the McKelvey flux method^{39,51} and the Landauer approach^{44,52} were next applied for electron and phonon transport. These flux-based methods provide approximate solutions to the Boltzmann equation with minimal computational expense due to simplifications to the collision term. The advantage of the present approach is that it provides the exact solution to the Boltzmann equation with the full collision matrix while requiring only discretization in wave vector space rather than in both real space and wave vector space. The resulting calculation is thus orders of magnitude faster than the purely numerical solution of the Boltzmann equation. Furthermore, compared to Ref. 47, this work does not require linearization around an equilibrium distribution and can accommodate arbitrary crystals, rather than isotropic ones, and the full collision matrix, rather than the relaxation time approximation. Finally, with a standard numerical treatment, the current density would need to be recomputed when any variable is changed. With the present method, the expensive calculation for $c1$ only needs to be performed once at each temperature because the average Knudsen number and $VBE$ are independent of $c1$ and included directly in Eq. (13).

## III. RESULTS

We apply this approach to study 1D steady-state electron transport in slabs of non-degenerate silicon ($NA=3\xd71018$ cm$\u22123$) and of width $L\u224840$ nm–$100$ $\mu $m. The hemispherical distribution function at the left and right boundaries is a Boltzmann distribution with carrier density $nL=nL,eqexp\u2061(qVBE/kT)$ and $nR=0$, respectively, where $nL,eq$ is the equilibrium minority carrier concentration. Note that at high biases, the law of the junction used for $nL$ may not be strictly satisfied due to degenerate carrier statistics, but this approximation will not affect our conclusions as discussed later. The junction is assumed to be a step junction with a uniform Ge fraction of 0.2, with bandgap values for SiGe obtained from Ref. 53. Owing to lack of precise knowledge of the electronic structure and scattering mechanisms in the strained SiGe:C films used in devices, in this work, we assume parabolic bands and employ the relaxation time approximation, with energy-dependent impurity and acoustic phonon scattering rates based on the forms given in Chap. 3 of Ref. 31. Several works report that the temperature dependence of minority carrier mobility in p-doped SiGe does not obey the $\u2248T1.5$ scaling typical of impurity scattering at cryogenic temperatures.^{28–30} Therefore, we extract coefficients and exponents for relaxation time $\tau =\mu 0TpEs$ by fitting the computed mobility to minority carrier mobility data reported in Fig. 3 of Ref. 30 and combining the relaxation times for acoustic and impurity scattering using Matthiessen’s rule. The fitting coefficients and exponents for acoustic (impurity) scattering rates are $\mu 0=350(700)$ cm$2$ V$\u22121$ s$\u22121$, $p=\u22121.5(0.25)$, and $s=\u22120.5(\u22120.5)$, respectively. Our conclusions are robust to these assumptions. We solve for $c1$ at the temperatures corresponding to Ref. 25, assuming an isotropic crystal using a $2200\xd72200$ grid in energy and an angular coordinate space. The collision matrix used in this work corresponds to Eq. (B.5) of Ref. 47. For example, at 300 K, we find $c1=0.7341$.

Figure 1(a) shows the carrier density vs normalized distance at 300 K for $L=40$ nm and $VBE=0.8$ V calculated at different orders of the expansion. The zeroth-order solution is by definition the drift–diffusion solution and exhibits a linear profile with a carrier density at the left edge of the slab of $nL=nL,eqexp\u2061(qVBE/kT)=0.87\xd71017$ cm$\u22123$. The higher-order non-diffusive corrections predict a reduced carrier concentration gradient across the base relative to the zeroth-order case. This feature occurs because at $\u03f5=0.38$, transport is quasiballistic, with electron mean free paths being on the order of the base width. Under these conditions, the ballistic propagation of some electrons across the base lowers the carrier density at the left boundary and increases it at the right boundary, resulting in jump-type boundary conditions at each edge and a shallower carrier density across the base. The solution at the first-order over-corrects for non-diffusive effects, but the infinite series accounts for higher-order terms and lies in between the zeroth and first-order solutions. This non-monotonicity occurs because the closed form solution switches signs at each order in Eq. (13).

To study how the electric current density is impacted by quasiballistic transport, we calculate the current density vs Knudsen number in Fig. 1(b) for two temperatures values reported in Ref. 25. The electric current density from Eq. (14) is proportional to the carrier concentration gradient in Fig. 1(a) for a given voltage, temperature, and Knudsen number. At 300 K, the current density is observed to decrease from the diffusive limit as $\u03f5\u21921$ (quasiballistic limit). As seen from Fig. 1(a), quasiballistic corrections lower the concentration gradient, and the contribution of the correction terms is multiplied by Knudsen number, resulting in greater deviation with increasing Knudsen number. Therefore, a shallower concentration gradient and consequently a decreased current density are expected with increasing Knudsen number. At 36 K corresponding to the temperature reported in Ref. 25, the calculation for $c1$ is repeated to give a value of $0.733$. At this temperature, we see that the deviation from the diffusive limit is comparable to that at 300 K. The difference between the two curves is due to changes in bulk diffusivity at 36 K that are independent of $c1$ and predict a slightly smaller current density than at 300 K. This calculation shows that even at cryogenic temperatures, quasiballistic effects yield a lower collector current density relative to the drift–diffusion prediction, opposite to that observed experimentally.

We now compare the predictions of our calculation with experimentally reported DC IV characteristics. The data from Ref. 25 are selected due to the availability of Gummel curves for a state-of-the-art device at a reasonably dense concentration of temperatures. Figure 2(a) plots calculated and measured collector current density vs base–emitter voltage $VBE$. We observe that at all temperatures, quasiballistic transport predicts a weakly reduced current density compared to the drift–diffusion prediction, as expected from the discussion of Fig. 1(b). The two calculations nearly coincide in Fig. 2(a) because the deviation of the quasiballistic current density from the diffusive current density is small on the logarithmic scale. At 300 K, 200 K, and 77 K, both predictions agree well with measured data. However, at 36 K, quasiballistic and drift–diffusion current density predictions are similar in magnitude but are both orders of magnitude lower than the measured data for a given $VBE$. Furthermore, the predicted slope of the $J$–$V$ curve at 36 K is unchanged from the drift–diffusion prediction and does not saturate to a value similar to that at 77 K as observed in the measured data.

To further examine this result, we present transconductance vs temperature at a fixed value of collector current density in Fig. 2(b). The transconductance per unit area $gm$ is calculated as the slope of the current–voltage curve at a fixed bias, and for an ideal SiGe HBT, it is given by $gm=qJC/kT$. The measured transconductance exhibits a plateau to a lower value than the ideal value as the temperature decreases. However, the transconductance values calculated from the slope of the calculated results in Fig. 2(a) are identical to the ideal values. Therefore, quasiballistic transport does not change the temperature dependence of the transconductance. The reason can be seen in Eq. (13) in which $VBE$ affects the electron density only through the left boundary condition $nL$; therefore, a change in voltage affects current density in a way that is independent of the mechanism of carrier transport and thus preserves the diffusive temperature dependence regardless of the Knudsen number.

## IV. DISCUSSION

Our analysis has shown that the predicted influence of quasiballistic transport on the IV characteristics of HBTs is inconsistent with the experimental observations. We now address several mechanisms that we have not included in our analysis. First, we have not incorporated the built-in field in the base region due to Ge grading. However, as in Ref. 54, the only effect of a field term is to modify $M$ and thus change the value of $c1$, and our conclusions are robust to such changes. Second, modern devices may have base doping levels exceeding $1019$ $cm\u22123$ to minimize base resistance and prevent carrier freeze-out.^{25,27,55} This high level of doping leads to degenerate carrier statistics, while we have assumed non-degenerate statistics. However, non-idealities in the IV characteristics are observed at temperatures up to $\u224880$ K for which the electrons are non-degenerate, and the present analysis thus applies under relevant conditions. Furthermore, in degenerate conditions, the law of the junction overpredicts the minority carrier density,^{56} implying that accounting for degenerate statistics would predict a further decreased collector current that again contradicts experiment. Last, the asymptotic solution is restricted to values of the average Knudsen number such that $2\u03f5c1<1$, which limits our analysis to a minimum base width of 40 nm for the chosen scattering rates. However, DC non-idealities were reported in devices with base widths on this order,^{25–27} and again, our analysis applies. Given these considerations, we conclude that quasiballistic transport is not responsible for non-ideal cryogenic DC characteristics of SiGe HBTs.

Finally, we discuss alternate explanations for the cryogenic non-idealities. Prior works have suggested that direct tunneling from the emitter to the collector is a possible origin of non-idealities in devices with narrow-base widths $\u224810$ nm.^{25,26} However, non-ideal cryogenic transconductance has been observed in first-generation SiGe HBTs with base widths on the order of $100$ nm^{24,27} for which the direct tunneling current is expected to be negligible. Other reported evidence for a transition from trap-assisted transport to a tunneling or quasiballistic mechanism is the change in the slope of the collector I–V characteristics.^{26,27} However, similar trends are observed in the forward regime of Schottky diodes^{49,57} for which direct tunneling or quasiballistic transport is not relevant, and the possibility of a similar mechanism occurring in SiGe HBTs has not been excluded.

We instead offer a hypothesis for the origin of cryogenic non-idealities as originating from spatial inhomogeneities in the base–emitter potential barrier. The earliest works on heavily doped p–n junctions reported temperature-independent slopes of cryogenic forward-bias I–V characteristics for voltages where band-to-band tunneling is unlikely to occur.^{58–60} Extensive studies of Schottky diodes have reported a variety of cryogenic $I$–$V$ non-idealities including the so-called $T0$ anomaly^{61} and temperature-dependent ideality factors.^{62–64} These non-idealities have been explained by inhomogeneities in barrier height^{49} that exist even in high-quality epitaxial junctions such as NiSi$2$/Si junctions and have been linked to the atomic structure of the interface.^{48,65} In SiGe HBT junctions, non-ideal base currents and the resulting random telegraph noise have been attributed to voltage barrier height fluctuations arising from trap states in the base–emitter space charge region^{50} such as those of electrically active carbon defects.^{66} Future work will examine whether such barrier inhomogeneities are able to explain the collector current non-idealities of SiGe HBTs.

## V. CONCLUSION

We have reported a study of quasiballistic transport in SiGe HBTs using an exact, semi-analytic solution to the Boltzmann equation based on an asymptotic expansion method. We find that the predicted IV characteristics including quasiballistic transport are inconsistent with experiment. Specifically, our calculations including quasiballistic transport predict collector currents that are orders of magnitude smaller than the measured currents for a given base–emitter voltage, and an unaltered temperature dependence of transconductance relative to the ideal value, both of which contradict experimental observations. We suggest that local fluctuations in the base–emitter barrier height could account for the non-ideal collector current as has been observed in Schottky diodes. Our work advances the understanding of electronic transport in cryogenic SiGe HBTs that is required to further improve their microwave noise performance.

## ACKNOWLEDGMENTS

The authors thank Mark Lundstrom, J. P. Peraud, and Nicolas Hadjiconstantinou for useful discussions. This work was supported by the National Science Foundation (NSF) (Award No. 1911926).

## DATA AVAILABILITY

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## REFERENCES

*2012 IEEE 12th Topical Meeting on Silicon Monolithic Integrated Circuits in RF Systems*(IEEE, 2012), pp. 133–136.

*2013 IEEE Bipolar/BiCMOS Circuits and Technology Meeting (BCTM)*(IEEE, 2013), pp. 215–218.

*2021 IEEE International Reliability Physics Symposium (IRPS)*(IEEE, 2021), pp. 1–6.

*2020 IEEE Radio Frequency Integrated Circuits Symposium (RFIC)*(IEEE, 2020), pp. 227–230.

*Technical Digest—International Electron Devices Meeting (IEDM)*(IEEE, 2017), Vol. 2, p. 3.1.1.

*IEEE MTT-S International Microwave Symposium Digest, 4–6 August 2020*(IEEE, 2020), p. 181.

*1997 21st International Conference on Microelectronics. Proceedings*(1997), Vol. 1, pp. 215–222; available at http://ir.lib.ncu.edu.tw:88/thesis/view_etd.asp?URN=91521052&fileName=GC91521052.pdf

*International Technical Digest on Electron Devices*(IEEE, 1990), Vol. 1, p. 171.

*Proceedings of the IEEE Bipolar/BiCMOS Circuits and Technology Meeting*(IEEE, 2017), p. 17.

*Technical Digest—International Electron Devices Meeting*(IEEE, 1988), p. 298.

*Monte Carlo Principles and Neutron Transport Problems*, Addison-Wesley Series in Computer Science and Information Processing (Addison-Wesley Publishing Company, 1969).

*2018 IEEE BiCMOS and Compound Semiconductor Integrated Circuits and Technology Symposium (BCICTS 2018)*(IEEE, 2018), p. 174.

*Rarefied Gas Dynamics*(Academic Press, 1969), p. 243.