Microreversibility, fluctuations, and nonlinear transport in transistors

TL;DR

This paper introduces a stochastic model for transistor charge transport that respects thermodynamic laws and microreversibility, verifying fundamental fluctuation relations and nonlinear transport properties through numerical analysis.

Contribution

It develops a comprehensive stochastic framework for transistors that incorporates detailed balance, fluctuation theorems, and nonlinear transport, advancing theoretical understanding.

Findings

Verification of the signal amplification under working conditions.

Confirmation of the fluctuation theorem for joint electric currents.

Numerical evidence of Onsager relations and their nonlinear extensions.

Abstract

We present a stochastic approach for charge transport in transistors. In this approach, the electron and hole densities are governed by diffusion-reaction stochastic differential equations satisfying local detailed balance and the electric field is determined with the Poisson equation. The approach is consistent with the laws of electricity, thermodynamics, and microreversibility. In this way, the signal amplifying effect of transistors is verified under their working conditions. We also perform the full counting statistics of the two electric currents coupled together in transistors and we show that the fluctuation theorem holds for their joint probability distribution. Similar results are obtained including the displacement currents. In addition, the Onsager reciprocal relations and their generalizations to nonlinear transport properties deduced from the fluctuation theorem are numerically shown to be satisfied.

Tables

Table 1. Table 1: The values of dimensionless physical quantities and parameters used in simulating the BJT model in rescaled units.

quantity	value	quantity	value \bigstrut
permittivity	$ϵ=0.01$	length of each cell	$\Delta x=0.1$ \bigstrut
elementary charge	$\left|e\right|=1.0$	width of each cell	$\Delta y=0.2$ \bigstrut
inverse temperature	$\beta =1.0$	number of cells in both $n$-type regions	$L_n=10$ \bigstrut
diffusion coefficient for electrons and holes	$D=0.01$	number of cells in the $p$-type region	$L_p=3$ \bigstrut
generation and recombination rate constants	$k_{+}=k_{-}=0.01$	number of cells in contact with the Base	$L_B=1$ \bigstrut

Table 2. Table 2: The set of parameter values used in Sec.  III .

parameter	value	parameter	value \bigstrut
volume of each cell	$\Omega ={10}^9$	section area	$\Sigma ={10}^{10}$, ${\Sigma }_B=5\times {10}^9$ \bigstrut
number of electrons for the Collector	${\overline{N}}_C={10}^{13}$	number of holes for the Collector	${\overline{P}}_C={10}^5$ \bigstrut
number of electrons for the Base	${\overline{N}}_B={10}^8$	number of holes for the Base	${\overline{P}}_B={10}^{10}$ \bigstrut
number of electrons for the Emitter	${\overline{N}}_E={10}^{13}$	number of holes for the Emitter	${\overline{P}}_E={10}^5$ \bigstrut

Table 3. Table 3: The set of parameter values used in Sec.  IV and Sec.  V .

parameter	value	parameter	value \bigstrut
volume of each cell	$\Omega =1000$	section areas	$\Sigma =10000$, ${\Sigma }_B=5000$ \bigstrut
number of electrons for the Collector	${\overline{N}}_C=10000$	number of holes for the Collector	${\overline{P}}_C=100$ \bigstrut
number of electrons for the Base	${\overline{N}}_B=100$	number of holes for the Base	${\overline{P}}_B=10000$ \bigstrut
number of electrons for the Emitter	${\overline{N}}_E=10000$	number of holes for the Emitter	${\overline{P}}_E=100$ \bigstrut

Table 4. Table 4: The comparison between the numerical affinities and their theoretical expectations. The statistics used to evaluate the numerical affinities is obtained by simulations with the time step d ​ t = 0.05 𝑑 𝑡 0.05 dt=0.05 , the total time t = 2.5 × 10 3 𝑡 2.5 superscript 10 3 t=2.5\times 10^{3} , and 5 × 10 5 5 superscript 10 5 5\times 10^{5} trajectories for every case.

case	$A_C^{\left(\mathrm{th}\right)}$	$A_C^{\left(\mathrm{num}\right)}$	$A_B^{\left(\mathrm{th}\right)}$	$A_B^{\left(\mathrm{num}\right)}$ \bigstrut
1	$1.0$	$0.9914\pm 0.0034$	$0.7$	$0.6942\pm 0.0027$ \bigstrut
2	$0.8$	$0.7919\pm 0.0025$	$0.4$	$0.3952\pm 0.0019$ \bigstrut
3	$0.5$	$0.5018\pm 0.0033$	$1.2$	$1.2007\pm 0.0041$ \bigstrut
4	$0.0$	$0.0000\pm 0.0000$	$0.0$	$0.0000\pm 0.0000$ \bigstrut
5	$-0.4$	$-0.4002\pm 0.0018$	$0.6$	$0.5975\pm 0.0020$ \bigstrut
6	$-0.5$	$-0.4864\pm 0.0029$	$-0.7$	$-0.6864\pm 0.0029$ \bigstrut
7	$-1.0$	$-1.0058\pm 0.0039$	$0.4$	$0.4022\pm 0.0028$ \bigstrut
8	$-1.2$	$-1.3924\pm 0.0084$	$1.3$	$1.4118\pm 0.0084$ \bigstrut

Table 5. Table 5: The numerical values of the quantities used in the fluctuation-dissipation and the Onsager reciprocal relations.

$L_{\alpha ,\beta }$	${D_{\alpha \beta }|}_{(0,0)}$	$L_{\alpha ,\beta }-{D_{\alpha \beta }|}_{(0,0)}$ \bigstrut
${\frac{\partial J_C}{\partial A_C}|}_{(0,0)}=93.106\pm 0.019$	${D_{CC}|}_{(0,0)}=92.991\pm 1.039$	$0.115$ \bigstrut
${\frac{\partial J_C}{\partial A_B}|}_{(0,0)}=-56.288\pm 0.019$	${D_{CB}|}_{(0,0)}=-56.343\pm 0.488$	$0.055$ \bigstrut
${\frac{\partial J_B}{\partial A_C}|}_{(0,0)}=-56.303\pm 0.020$	${D_{BC}|}_{(0,0)}=-56.343\pm 0.488$	$0.040$ \bigstrut
${\frac{\partial J_B}{\partial A_B}|}_{(0,0)}=112.603\pm 0.020$	${D_{BB}|}_{(0,0)}=113.158\pm 0.487$	$-0.555$ \bigstrut

Table 6. Table 6: The numerical values of the quantities used in the nonlinear transport relations ( 90 ).

$M_{\alpha ,\beta \gamma }$	$R_{\alpha \beta ,\gamma }$	$R_{\alpha \gamma ,\beta }$	$M_{\alpha ,\beta \gamma }-R_{\alpha \beta ,\gamma }-R_{\alpha \gamma ,\beta }$ \bigstrut
${\frac{{\partial }^2{J}_C}{\partial A_C^2}|}_{(0,0)}=-67.388\pm 0.620$	${\frac{\partial D_{CC}}{\partial A_C}|}_{(0,0)}=-33.642\pm 9.897$	${\frac{\partial D_{CC}}{\partial A_C}|}_{(0,0)}=-33.642\pm 9.897$	$-0.104$ \bigstrut
${\frac{{\partial }^2{J}_C}{\partial A_B^2}|}_{(0,0)}=-45.325\pm 0.620$	${\frac{\partial D_{CB}}{\partial A_B}|}_{(0,0)}=-22.474\pm 4.639$	${\frac{\partial D_{CB}}{\partial A_B}|}_{(0,0)}=-22.474\pm 4.639$	$-0.377$ \bigstrut
${\frac{{\partial }^2{J}_C}{\partial A_C\partial A_B}|}_{(0,0)}=68.747\pm 0.097$	${\frac{\partial D_{CC}}{\partial A_B}|}_{(0,0)}=47.409\pm 9.900$	${\frac{\partial D_{CB}}{\partial A_C}|}_{(0,0)}=20.992\pm 4.642$	$0.346$ \bigstrut
${\frac{{\partial }^2{J}_B}{\partial A_C^2}|}_{(0,0)}=42.064\pm 0.667$	${\frac{\partial D_{CB}}{\partial A_C}|}_{(0,0)}=20.992\pm 4.642$	${\frac{\partial D_{CB}}{\partial A_C}|}_{(0,0)}=20.992\pm 4.642$	$0.080$ \bigstrut
${\frac{{\partial }^2{J}_B}{\partial A_B^2}|}_{(0,0)}=90.066\pm 0.665$	${\frac{\partial D_{BB}}{\partial A_B}|}_{(0,0)}=45.068\pm 4.644$	${\frac{\partial D_{BB}}{\partial A_B}|}_{(0,0)}=45.068\pm 4.644$	$-0.070$ \bigstrut
${\frac{{\partial }^2{J}_B}{\partial A_C\partial A_B}|}_{(0,0)}=-44.777\pm 0.107$	${\frac{\partial D_{CB}}{\partial A_B}|}_{(0,0)}=-22.474\pm 4.639$	${\frac{\partial D_{BB}}{\partial A_C}|}_{(0,0)}=-22.330\pm 4.630$	$0.027$ \bigstrut

Equations

$$d(\mathbf{r})=d\theta \left(-x-l_p/2\right)+d\theta \left(x-l_p/2\right)\text{,}$$

$$a(\mathbf{r})=a\theta \left(x+l_p/2\right)\theta \left(-x+l_p/2\right)\text{,}$$

$$\rho =e(p-n+d-a),$$

$$n_{\mathrm{eq}}(\mathbf{r})\sim {\mathrm{e}}^{+\beta {\varphi }_{\mathrm{eq}}(\mathbf{r})}\text{and}{p}_{\mathrm{eq}}(\mathbf{r})\sim {\mathrm{e}}^{-\beta e{\varphi }_{\mathrm{eq}}(\mathbf{r})}$$

$$({\varphi }_C-{\varphi }_E{)}_{\mathrm{eq}}=\frac{1}{\beta e}\ln \frac{n_C}{n_E}=\frac{1}{\beta e}\ln \frac{p_E}{p_C}$$

$$({\varphi }_B-{\varphi }_E{)}_{\mathrm{eq}}=\frac{1}{\beta e}\ln \frac{n_B}{n_E}=\frac{1}{\beta e}\ln \frac{p_E}{p_B}\text{,}$$

$$V_C={\varphi }_C-{\varphi }_E-\frac{1}{\beta e}\ln \frac{n_C}{n_E}\text{,}$$

$$V_B={\varphi }_B-{\varphi }_E-\frac{1}{\beta e}\ln \frac{n_B}{n_E}\text{,}$$

$$A_C\equiv \beta e{V}_C\text{mbox}and{A}_B\equiv \beta e{V}_B\text{,}$$

$${\mu }_n=\beta e{D}_n\text{and}{\mu }_p=\beta e{D}_p\text{.}$$

$${\partial }_{t}n+\nabla \cdot {\mathbf{j}}_n={\sigma }_n\text{,}$$

$${\partial }_{t}p+\nabla \cdot {\mathbf{j}}_p={\sigma }_p\text{,}$$

$${\nabla }^2\varphi =-\frac{\rho }{ϵ}\text{,}$$

$${\sigma }_n={\sigma }_p=k_{+}-k_{-}np+\delta \sigma \text{,}$$

$${\mathbf{j}}_n=-{\mu }_{n}n\mathcal{E}-D_n\nabla n+\delta {\mathbf{j}}_n\text{,}$$

$${\mathbf{j}}_p=+{\mu }_{p}p\mathcal{E}-D_p\nabla p+\delta {\mathbf{j}}_p\text{,}$$

$$\mathcal{E}=-\nabla \varphi \text{,}$$

$$⟨\delta {\mathbf{j}}_n(\mathbf{r},t)⟩=⟨\delta {\mathbf{j}}_p(\mathbf{r},t)⟩=0\text{,}$$

$$⟨\delta {\mathbf{j}}_n(\mathbf{r},t)\otimes \delta {\mathbf{j}}_n({\mathbf{r}}^{\prime },t^{\prime })⟩={\Gamma }_{nn}(\mathbf{r},t){\delta }^3(\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-t^{\prime })1\text{,}$$

$$⟨\delta {\mathbf{j}}_p(\mathbf{r},t)\otimes \delta {\mathbf{j}}_p({\mathbf{r}}^{\prime },t^{\prime })⟩={\Gamma }_{pp}(\mathbf{r},t){\delta }^3(\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-t^{\prime })1\text{,}$$

$$⟨\delta {\mathbf{j}}_n(\mathbf{r},t)\otimes \delta {\mathbf{j}}_p({\mathbf{r}}^{\prime },t^{\prime })⟩=0\text{,}$$

$$⟨\delta \sigma (\mathbf{r},t)⟩=0\text{,}$$

$$⟨\delta \sigma (\mathbf{r},t)\delta \sigma ({\mathbf{r}}^{\prime },t^{\prime })⟩={\Gamma }_{\sigma \sigma }(\mathbf{r},t){\delta }^3(\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-t^{\prime })\text{,}$$

$$⟨\delta \sigma (\mathbf{r},t)\delta {\mathbf{j}}_n({\mathbf{r}}^{\prime },t^{\prime })⟩=⟨\delta \sigma (\mathbf{r},t)\delta {\mathbf{j}}_p({\mathbf{r}}^{\prime },t^{\prime })⟩=0\text{,}$$

$${\Gamma }_{nn}(\mathbf{r},t)\equiv 2D_{n}n(\mathbf{r},t)\text{,}$$

$${\Gamma }_{pp}(\mathbf{r},t)\equiv 2D_{p}p(\mathbf{r},t)\text{,}$$

$${\Gamma }_{\sigma \sigma }(\mathbf{r},t)\equiv k_{+}+k_{-}n(\mathbf{r},t)p(\mathbf{r},t)$$

$${\bar{N}}_C=n_C\Omega \text{,}{\bar{P}}_C=p_C\Omega \text{,}$$

$${\bar{N}}_B=n_B\Omega \text{,}{\bar{P}}_B=p_B\Omega \text{,}$$

$${\bar{N}}_E=n_E\Omega \text{,}{\bar{P}}_E=p_E\Omega \text{.}$$

$$n_C{p}_C=n_B{p}_B=n_E{p}_E={\nu }^2\text{.}$$

$$n_C=n_E\text{,}{p}_C=p_E\text{,}$$

$$\alpha ={\left(\frac{\partial J_C}{\partial J_B}\right)}_{A_C}$$

$$\alpha (A_C=20,A_B=6)\simeq 4.278\text{,}$$

$${\mathcal{I}}_C(t)=\sum _{n=-\infty }^{+\infty }{q}_n^{(C)}\delta (t-t_n^{(C)})\text{,}$$

$${\mathcal{I}}_B(t)=\sum _{n=-\infty }^{+\infty }{q}_n^{(B)}\delta (t-t_n^{(B)})\text{,}$$

$$Z_C(t)=\frac{1}{e}{\int }_0^t{\mathcal{I}}_C(t^{\prime })d{t}^{\prime }\text{,}{Z}_B(t)=\frac{1}{e}{\int }_0^t{\mathcal{I}}_B(t^{\prime })d{t}^{\prime }.$$

$${\tilde{\mathcal{I}}}_C(t)={\mathcal{I}}_C(t)-ϵ{\partial }_t{\partial }_x\varphi {\Sigma }_C\text{,}$$

$${\tilde{\mathcal{I}}}_B(t)={\mathcal{I}}_B(t)-ϵ{\partial }_t{\partial }_y\varphi {\Sigma }_B\text{,}$$

$$J_C\equiv \underset{t\to \infty }{\lim }\frac{1}{t}⟨Z_C(t)⟩=\underset{t\to \infty }{\lim }\frac{1}{t}⟨{\tilde{Z}}_C(t)⟩\text{,}$$

Full text

I Introduction

Transistors are the main compounds of semiconductor electronic technology. The core of transistors is composed of three semiconducting materials concatenated in series, thus forming double junctions. The middle semiconductor is doped with charged impurities different from those in the two other semiconductors. Since transistors have three ports and currents flow between pairs of ports, two electric currents are coupled together inside transistors, enabling the amplification of signals SST51 ; EM54 ; SS04 ; CC05 ; B05 ; SN07 .

The fundamental issue is that the coupling between the electric currents is ruled by microreversibility, as in any type of device or process. In linear regimes close to thermodynamic equilibrium, microreversibility implies the Onsager-Casimir reciprocal relations O31a ; O31b ; C45 . However, transistors are functioning in highly nonlinear regimes beyond the domain of application of the Onsager-Casimir reciprocal relations. Remarkably, the generalizations of these relations beyond the linear regime are known today S92 ; AG04 ; AG07JSM ; HPPG11 ; BG18 . They can be deduced from the fluctuation theorem for currents, which is based on the time-reversal symmetry of the microscopic dynamics of electrons and ions AG07JSP ; AGMT09 ; AG09 ; EHM09 ; CHT11 ; S12 ; G13 . The fluctuation theorem is valid not only in the linear regimes, but also in the nonlinear regimes, and can thus be used to investigate the nonlinear transport properties of transistors.

In our previous paper GG18 , the fluctuation theorem was considered for diodes that are also nonlinear electronic devices. Here, our purpose is to extend these considerations to transistors. The novel aspect is that two currents are flowing in transistors, instead of only one in diodes. As a consequence of the nonlinear coupling between the two currents, the generalizations of Onsager-Casimir reciprocal relations to nonlinear transport can be tested in transistors.

For this purpose, the stochastic approach of Ref. GG18 is extended from the single junction of diodes to the double junction of $n$-$p$-$n$ transistors. The approach is based on diffusion-reaction stochastic partial differential equations for electrons and holes, including their Coulomb interaction described by the Poisson equation. This scheme satisfies local detailed balance in consistency with microreversibility. The stochastic description is presented in Sec. II. The functionality of transistors is studied in Sec. III. Section IV is devoted to the fluctuation theorem for the two currents of the transistor. Section V shows that the linear response coefficients obey the Onsager-Casimir reciprocal relation and the fluctuation-dissipation theorem, and that the next-order nonlinear response coefficients satisfy higher-order generalizations. Section VI gives concluding remarks.

II Stochastic description of transistors

II.1 The bipolar $n$-$p$-$n$ junction transistor

There exists many types of transistors SS04 ; CC05 ; B05 ; SN07 . The bipolar $n$-$p$-$n$ junction transistor (BJT) is one of the most common of them. BJTs consist of three small doped regions of a piece of silicon, respectively typed as $n$, $p$, and $n$, thus forming two junctions, as shown in Fig. 1. The electrons ${\rm e}^{-}$ and holes ${\rm h}^{+}$ are the two mobile charge carriers across the bipolar $n$-$p$-$n$ junction, with electrons being the majority ones in $n$-type semiconductor, and holes the majority ones in $p$-type semiconductor. The positively-charged donors and negatively-charged acceptors are respectively anchored in $n$-type semiconductors and $p$-type semiconductors. Each doped region has a port and the three ports are in contact with some charge carrier reservoir. They are respectively called Collector, Base, and Emitter (see Fig. 1).

In order to model the transistor, a Cartesian coordinate system is associated with the system. As shown in Fig. 1(b), the semiconducting material extends from $x=-l/2$ to $x=+l/2$ and is divided in three parts. The part from $x=-l/2$ to $x=-l_{p}/2$ is of $n$-type, the one from $x=-l_{p}/2$ to $x=+l_{p}/2$ of $p$-type, and the one from $x=+l_{p}/2$ to $x=+l/2$ of $n$-type. The three parts are respectively of lengths $l_{n}=(l-l_{p})/2$, $l_{p}$, and $l_{n}=(l-l_{p})/2$. The Collector is in contact at $x=-l/2$, the Emitter at $x=+l/2$, and the Base along a length $l_{B}$ symmetrically located around the origin $x=0$. The length of the contact with the Base is smaller than the one of the $p$-type part: $l_{B}<l_{p}$. The geometry is chosen to be symmetric with respect to $x=0$ for simplicity.

In addition, the bipolar $n$-$p$-$n$ double junction has the section area $\Sigma$ in the transverse $y$- and $z$-directions. The section areas of the contacts with the Collector and Emitter are assumed to be equal: $\Sigma_{C}=\Sigma_{E}=\Sigma$. Accordingly, the semiconducting material extends over a domain of volume $V=l\Sigma$. Moreover, we denote $\Sigma_{B}$ the section area of the contact with the Base.

The donor density $d({\bf r})$ and acceptor density $a({\bf r})$ are supposed to be uniform in the different types of semiconductor. Therefore, they can be expressed as

[TABLE]

in terms of two constant values $a$ and $d$, combined with Heaviside’s step function $\theta(x)$ defined such that $\theta(x)=1$ if $x>0$ and $\theta(x)=0$ otherwise. The charge density is thus given by

[TABLE]

with the elementary electric charge $e=|e|$, and the densities of holes $p$, electrons $n$, donors $d$, and acceptors $a$. Here, we have assumed that every donor gives one electron and every acceptor one hole. Because of the electrostatic interaction between the charges, these densities are coupled to the electric potential $\phi({\bf r})$.

The electron and hole densities as well as the electric potential have fixed boundary values at the contacts with the three reservoirs. They are respectively given by $n_{C}$, $p_{C}$, $\phi_{C}$ at the Collector; $n_{B}$, $p_{B}$, $\phi_{B}$ at the Base; and $n_{E}$, $p_{E}$, $\phi_{E}$ at the Emitter.

If the transistor is at equilibrium without flow of charge carriers, detailed balance between the generation and recombination of electron-hole pairs requires that $n_{\rm eq}p_{\rm eq}=\nu^{2}$, where $\nu$ is called the intrinsic carrier density. Moreover, the electron and hole densities are given at equilibrium by

[TABLE]

in terms of the electric potential determined across the whole system by the Poisson equation and the boundary conditions at the contacts with the three reservoirs. If the BJT is at equilibrium, the inhomogeneous distributions of the charge carriers thus produce the Nernst potentials

[TABLE]

and

[TABLE]

where $\beta\equiv(k_{\rm B}T)^{-1}$ is the inverse temperature.

The transistor is driven out of equilibrium by applying voltage differences with respect to the Nernst potentials

[TABLE]

which induce currents across the BJT. In the following, we use the associated affinities or thermodynamic forces

[TABLE]

which are dimensionless. The equilibrium state is recovered if they vanish, i.e., if the applied voltages are equal to zero $V_{C}=V_{B}=0$.

II.2 Stochastic diffusion-reaction equations

The thermal agitation inside the BJT generates incessant erratic motion for the electrons and holes, in turn causing local fluctuations in the currents and reaction rates. These fluctuations can be described within the stochastic approach by introducing Gaussian white noise fields in the diffusion-reaction equations for the electron and hole densities. The advantage of this approach is that the usual phenomenological parameters suffice for the stochastic description.

The mobilities of electrons and holes are related with their diffusion coefficients through Einstein’s relations

[TABLE]

Besides, the electron-hole pairs are randomly generated and recombined according to the reactions

where $k_{+}$ and $k_{-}$ are respectively the generation and recombination rate constants. In general, the quantities $D_{n}$, $D_{p}$, and $k_{\pm}$ are spatially dependent in an inhomogeneous medium. However, for simplicity, we assume that they are uniform across the whole BJT.

Considering the diffusion and generation-recombination processes as well as the electrostatic interaction between the charges, we have the following stochastic partial differential equations for the charge carrier densities coupled to the Poisson equation for the electric potential,

[TABLE]

where

[TABLE]

are the reaction rates, the current densities, and the electric field, while $\rho$ is the charge density given by Eq. (3) and $\epsilon$ the dielectric constant of the material GG18 . The fluctuations $\delta{\bf j}_{n}$, $\delta{\bf j}_{p}$, and $\delta\sigma$ are Gaussian white noise fields characterized by

[TABLE]

where ${\boldsymbol{\mathsf{1}}}$ is the $3\times 3$ identity matrix and

[TABLE]

are the noise spectral densities associated with the electron and hole diffusions, and the reaction.

Because of Eqs. (18) and (22), we recover the mean-field equations of the macroscopic description by averaging the stochastic partial differential equations over the noises.

II.3 Numerical method for simulating the transistor

For the numerical simulation of the transistor, a Markov jump process is associated with the stochastic partial differential equations (11)-(17), as described in detail in Appendix A. Space is discretized into $L$ cells of length $\Delta x=l/L$, section area $\Sigma$, and volume $\Omega=\Sigma\Delta x$, located at the coordinates $x_{i}=(i-0.5)\Delta x-l/2$ ($i=1,2,\dots,L$). Consistently with Fig. 1(b), there are $L_{n}=l_{n}/\Delta x$ cells in both parts of $n$-type, $L_{p}=l_{p}/\Delta x$ cells for the part of $p$-type, and $L_{B}=l_{B}/\Delta x$ cells in contact with the Base. The numbers of electrons, holes, acceptors, and donors in each cell of the BJT are related to the corresponding densities by $N_{i}=n(x_{i})\Omega$, $P_{i}=p(x_{i})\Omega$, $A_{i}=a(x_{i})\Omega$, and $D_{i}=d(x_{i})\Omega$. The state of the discretized BJT is fully characterized by the electron numbers ${\bf N}=(N_{i})_{i=1}^{L}$ and the hole numbers ${\bf P}=(P_{i})_{i=1}^{L}$ in the cells. The master equation ruling the time evolution of their probability distribution ${\cal P}({\bf N},{\bf P},t)$ is given in Appendix A.1. Moreover, the Poisson equation (13) is also discretized along the chain of $L$ cells forming the system, taking into account the electric potentials of the Collector, the Base, and the Emitter, as explained in Appendix A.2. The resulting stochastic process can be simulated numerically by Gillespie’s algorithm G76 , which is an exact method for generating random trajectories in this case.

In order to speed up the simulation, the Markov jump process is approximated by a Langevin stochastic process under the assumption that the numbers of electrons and holes are large enough in every cell, $N_{i}\gg 1$ and $P_{i}\gg 1$. Accordingly, these numbers obey stochastic differential equations expressed in terms of the fluxes of particles between the cells, the reaction rates, and Gaussian white noises for their fluctuations, as shown in Appendix B.

At the contacts with the three reservoirs, the boundary conditions on the charge carrier densities determine the boundary values for the corresponding particle numbers

[TABLE]

Furthermore, the three parts of the transistor are supposed to be doped from a semiconducting material of uniform intrinsic density $\nu$, so that the boundary values of the electron and hole densities should satisfy the conditions

[TABLE]

We further set

[TABLE]

to have a system that is symmetric with respect to $x=0$, as depicted in Fig. 1(b).

In numerical simulations, the statistical averages of any observable quantity $X$ can be evaluated by the time average $\langle X\rangle=\lim_{T\to\infty}(1/T)\int_{0}^{T}X(t)\,dt$, which is equivalent by ergodicity to the ensemble average $\langle X\rangle=\sum_{{\bf N},{\bf P}}X\,{\cal P}_{\rm st}({\bf N},{\bf P})$ over the stationary probability distribution ${\cal P}_{\rm st}$. In the continuum limit, the volume of the cells is supposed to vanish together with the particle numbers, so that the electron and hole densities can be recovered as $n(x_{i})=N_{i}/\Omega$ and $p(x_{i})=P_{i}/\Omega$.

\bigstrutjot

=2pt

We assume for simplicity that the electron and hole diffusion coefficients are equal $D_{n}=D_{p}\equiv D$. As done in our previous paper GG18 , the quantities of interest may be rescaled using the intrinsic carrier density $\nu$, the intrinsic carrier lifetime $\tau=1/(k_{-}\nu)$, the intrinsic carrier diffusion length before recombination $\ell=\sqrt{D\tau}$, the inverse temperature $\beta$, and the elementary electric charge. After this rescaling, the quantities of interest become dimensionless. Table 1 gives the values of the so-rescaled quantities used in the following numerical simulations of the BJT model.

III The Functionality of transistors

The purpose of this section is to show that the properties characterizing the functionality of transistors can be described within the stochastic approach.

In electronic technology, transistors are primarily used to amplify signals in electric circuits. This amplification results from the coupling between the two electric currents, $J_{C}$ and $J_{B}$. By this coupling, one current can serve as input and the other as output. The amplification factor is defined as the ratio of these two currents, $J_{C}/J_{B}$. We may also introduce the differential amplification factor as follows. When the affinity $A_{C}$ is fixed, the variation of the other affinity $A_{B}$ leads to variations of $J_{C}$ and $J_{B}$. The amplification factor is defined as the ratio between these two variations

[TABLE]

under specific working conditions. To achieve the functionality of signal amplification, the transistors should satisfy the following requirements:

• •

The concentration of the majority charge carriers in the Collector region should be overwhelmingly larger than the concentration of minority charge carriers in the Base region.

• •

The concentration of the majority charge carriers in the Emitter region should be overwhelmingly larger than the concentration of minority charge carriers in the Base region.

• •

The Collector-Base junction should be reverse biased.

• •

The Emitter-Base junction should be forward biased.

• •

The Base region should be very thin so that the majority charge carriers in the Emitter region can easily get swept to the Collector region.

• •

The contacting section areas $\Sigma_{C}$ and $\Sigma_{E}$ should be larger than $\Sigma_{B}$.

Table 2 gives a set of parameter values approaching these requirements in order to show that the present stochastic model can describe transistors in such regimes. The first two conditions are satisfied since $\bar{N}_{C}=\bar{N}_{E}\gg\bar{N}_{B}$, and the last one because $\Sigma=\Sigma_{C}=\Sigma_{E}>\Sigma_{B}$.

If the transistor was at equilibrium without applied voltage ($A_{C}=A_{B}=0$), the Nernst potentials (5) and (6) would take the values $(\phi_{C}-\phi_{E})_{\rm eq}=0$ and $(\phi_{B}-\phi_{E})_{\rm eq}=-11.5$ with the parameter set of Table 2. At equilibrium, the electric field would have a symmetric profile around $x=0$ with $(\phi_{C}-\phi_{B})_{\rm eq}=(\phi_{E}-\phi_{B})_{\rm eq}=11.5$.

Figure 2 shows the profiles of charge carrier densities and current densities together with the electric potential under nonequilibrium conditions with applied voltages corresponding to $A_{C}=20$ and $A_{B}=6$. In Fig. 2(a), we see that the Base region is thin in the model, so that the fifth condition is satisfied. As observed in Fig. 2(b), the current densities are non-vanishing because the transistor is out of equilibrium. According to Eqs. (7)-(8), we here have that $\phi_{C}-\phi_{E}=20$ and $\phi_{B}-\phi_{E}=-5.5$, so that $\phi_{C}-\phi_{B}=25.5$ and $\phi_{E}-\phi_{B}=5.5$, in agreement with the electric field plotted in Fig. 2(c). Since $\phi_{C}-\phi_{B}=25.5$ is larger than $(\phi_{C}-\phi_{B})_{\rm eq}=11.5$, the Collector-Base junction is reverse biased, as it should by the third condition. Moreover, $\phi_{E}-\phi_{B}=5.5$ is smaller than $(\phi_{E}-\phi_{B})_{\rm eq}=11.5$, so that the Emitter-Base junction is forward biased and the fourth condition is also satisfied. Under these conditions, the transistor can indeed achieve signal amplification, as demonstrated in Fig. 3. The currents $J_{C}$ and $J_{B}$ are shown in Fig. 3(a) as functions of $A_{B}$, with $A_{C}$ fixed. Since the current $J_{C}$ is greater than $J_{B}$, the amplification factor $J_{C}/J_{B}$ is larger than unity, as expected. Furthermore, Fig. 3(b) depicts how the current $J_{C}$ increases with the other current $J_{B}$ and the associated affinity $A_{B}$. For $A_{B}=6$, the differential amplification factor (44) is evaluated to be

[TABLE]

which is also larger than unity, as required. It should be noticed that the amplification factors can take different values for different working conditions of the transistor.

These results show that the stochastic approach is relevant to study transistors in their regimes of signal amplification. We proceed in the next Sec. IV and Sec. V with the study of their fluctuation properties.

IV Fluctuation Theorem for Currents

IV.1 Generalities

We consider the fluctuating electric currents flowing respectively across the contact with the Collector and the contact with the Base. These electric currents are due to the random motion of electrons and holes crossing the contact sections between the transistor and the corresponding reservoirs. The instantaneous electric currents are thus defined as

[TABLE]

where $t_{n}^{(C)}$ (resp. $t_{n}^{(B)}$) are the random times of the crossing events and $q_{n}^{(C)}$ (resp. $q_{n}^{(B)}$) are the transferred charges equal to $\pm e$ depending on whether the carrier is an electron or a hole and if its motion is inward or outward the transistor. The corresponding random numbers of charges accumulated over the time interval $[0,\,t]$ are defined as

[TABLE]

We also define the instantaneous total electric currents including the contribution of displacement currents as

[TABLE]

which are the experimentally measured electric currents GG18 ; AG09 ; BB00 ; S38 ; R39 , as well the corresponding accumulated charge numbers $\tilde{Z}_{C}(t)$ and $\tilde{Z}_{B}(t)$ with definitions as in Eq. (48).

The mean values of the charge currents are given by

[TABLE]

and the corresponding electric currents by $I_{C}=eJ_{C}$ and $I_{B}=eJ_{B}$. The equality between the mean values without and with the displacement currents comes from the fact that the displacement currents are given by a time derivative.

The diffusivities of the currents are defined as

[TABLE]

in terms of the variances and the covariances between the accumulated random charge numbers

[TABLE]

The diffusivities also take the same value whether the displacement currents are included or not. Since the covariance between two random variables is symmetric under their exchange, we have the symmetry $D_{CB}=D_{BC}$.

We suppose that the voltages (7) and (8) are applied at the boundaries of the transistor. Consequently, the transistor is driven out of equilibrium and the stochastic process of charge transfers between the reservoirs eventually reaches a nonequilibrium steady state. This latter is expected to depend on the applied voltages, or equivalently on the affinities

[TABLE]

which are determined by the differences of electrochemical potentials between the corresponding reservoirs. The dependences of the mean values of the currents on the affinities define the characteristic functions of the transistor: $J_{C}(A_{C},A_{B})$ and $J_{B}(A_{C},A_{B})$. At equilibrium, the affinities are vanishing together with the applied voltages and the mean values of the currents, so that $J_{C}(0,0)=J_{B}(0,0)=0$. However, the diffusivities do not necessarily vanish at equilibrium.

Beyond the mean values of the currents and the diffusivities, the process can be characterized by higher cumulants or the full probability distribution $P_{A_{C},A_{B}}(Z_{C},Z_{B},t)$ that $Z_{C}(t)$ and $Z_{B}(t)$ charges are crossing the Collector and the Base during the time interval $[0,t]$, while the transistor is in a nonequilibrium steady state of affinities $A_{C}$ and $A_{B}$. This steady state is given by the stationary solution of the master equation of the Markov jump process described in Appendix A. Using the network representation of this Markov jump process and its decomposition into cyclic paths S76 , the process can be shown to obey a fluctuation theorem for all the currents as a consequence of local detailed balance AG07JSP ; AG09 . This theorem states that the joint distribution of random variables $Z_{C}$ and $Z_{B}$ at time $t$ satisfies the following fluctuation relation

[TABLE]

A similar fluctuation relation holds if the displacement currents are included in the accumulated charge numbers AG09 .

As a consequence of the fluctuation theorem, the thermodynamic entropy production is always non-negative in accord with the second law of thermodynamics. The entropy production can indeed be expressed as the Kullback-Leibler divergence between the probability distributions of opposite fluctuations of the currents G13 , giving the dissipated power divided by the thermal energy

[TABLE]

as expected.

We notice that the fluctuation relation (61) holds in the long-time limit. The convergence time is determined by diffusion GGHK18 and it can be estimated to range between the time of diffusion across the middle part, $t_{\rm diff}\sim l_{p}^{2}/D\sim 9$, and the one before recombination, $t_{\rm diff}\sim\ell^{2}/D\sim 100$.

IV.2 Numerical results

The direct test of the fluctuation relation (61) requires the availability of an overlap between the probability distributions $P(Z_{C},Z_{B},t)$ and $P(-Z_{C},-Z_{B},t)$. Since the maxima of these distributions move apart under nonequilibrium conditions, the overlap rapidly decreases as time increases. Therefore, the direct test of the fluctuation relation is restricted to short times. Nevertheless, the test is possible as shown in Fig. 4 for the joint probability distributions of the accumulated charge numbers without and with the displacement currents using the set of parameter values given in Table 3. For the bare charge numbers, Fig. 4(a) depicts the joint distribution itself at time $t=20$, which is roughly Gaussian and shifted with respect to the origin because of the elapsed time. There is a significant overlap with the opposite distribution $P(-Z_{C},-Z_{B},t)$ and Fig. 4(b) shows several contours of the two-dimensional function $\ln\left[P(Z_{C},Z_{B},t)/P(-Z_{C},-Z_{B},t)\right]$ in the plane of the variables $Z_{C}$ and $Z_{B}$. These contours appear straight given the presence of statistical errors, in agreement with the prediction of the fluctuation theorem that the function should be linear. The function $\ln\left[P(Z_{C},Z_{B},t)/P(-Z_{C},-Z_{B},t)\right]$ can thus be fitted to a linear function $A_{C}(t)\,Z_{C}+A_{B}(t)\,Z_{B}$, defining the finite-time affinities $A_{C}(t)$ and $A_{B}(t)$. However, their values remain smaller than the applied affinities $A_{C}=A_{B}=0.1$ because convergence is expected for $t\gg t_{\rm diff}$ and has not yet been reached in Fig. 4.

As shown in Fig. 4(c) and Fig. 4(d), similar results hold for the joint probability distribution $P(\tilde{Z}_{C},\tilde{Z}_{B},t)$ of the charge numbers with the displacement currents. As seen in Fig. 4(c), the displacement currents have for effect that the distribution $P(\tilde{Z}_{C},\tilde{Z}_{B},t)$ is narrower than $P(Z_{C},Z_{B},t)$ depicted in Fig. 4(a). Consequently, the finite-time affinities $\tilde{A}_{C}(t)$ and $\tilde{A}_{B}(t)$ are larger than $A_{C}(t)$ and $A_{B}(t)$ and the convergence in time towards the asymptotic values of the affinities should be faster for the statistics of the transferred total charges $\tilde{Z}_{C}$ and $\tilde{Z}_{B}$ including the displacement currents, than for the statistics of the transferred charges $Z_{C}$ and $Z_{B}$. Figure 5 confirms that the finite-time affinities $\tilde{A}_{C}(t)$ and $\tilde{A}_{B}(t)$ approach their asymptotic value $A_{C}=A_{B}=0.1$, as time increases. Since the overlap between the opposite distributions rapidly decreases, statistical errors increase for $t>20$. The exponential fits of the finite-time affinities provide estimations of the convergence times in the range of values expected by charge carrier diffusion.

In order to test the convergence of the finite-time affinities towards their asymptotic values over longer time scales, we develop a method using the following coarse-grained model,

[TABLE]

where the charges are supposed to be transferred between the three reservoirs with the transition rates $\{W_{kl}\}_{k,l=C,B,E}$, as formulated in Appendix C. This constitutes the minimal model in the sense that the values of its rates can be fully determined from the knowledge of the mean currents and diffusivities, if the conditions of local detailed balance are satisfied. This simple model is related to the Ebers-Moll transport model of bipolar junction transistors EM54 ; SS04 . Given the values $J_{C}$, $J_{B}$, $D_{CC}$, $D_{BB}$, and $D_{CB}$ of the mean currents and the diffusivities, the six rates $W_{kl}$ can be determined, giving the values of the affinities according to $A_{kl}=\ln(W_{kl}/W_{lk})$ with $k,l=C,B,E$. Since this model results from the coarse graining of the complete description, it has a domain of validity limited to moderate values of the applied voltages. In this domain, the parameter values of the model can thus be fitted to the numerical values of the mean currents (51)-(52) and the diffusivities (53)-(55) of the full model in order to obtain the affinities.

Table 4 shows the comparison between the numerical affinities and the theoretical predictions for several cases. Accurate agreement is found if the affinities remain moderate, confirming the convergence of the finite-time affinities $A_{C}(t)$ and $A_{B}(t)$ towards their expected asymptotic values (59) and (60) within the domain of validity of the model (80). Despite the limited scope of application of this method, the agreement between the numerical and theoretical values of the affinities brings further numerical support to the fluctuation relation for the currents. In the next section, the consequences of the fluctuation theorem on the linear and nonlinear transport properties will be tested.

V Linear and nonlinear response properties

V.1 Deduction of the properties from the fluctuation theorem

The fluctuation theorem provides a unified framework for deducing the Onsager reciprocal relations and their generalizations to the nonlinear transport properties S92 ; AG04 ; AG07JSM ; HPPG11 ; BG18 . For this purpose, it is convenient to introduce the cumulant generating function

[TABLE]

where ${\boldsymbol{\lambda}}=(\lambda_{C},\lambda_{B})$ are the so-called counting parameters and the macroscopic affinities are written in vectorial notation ${\bf A}=(A_{C},A_{B})$. As a consequence of the fluctuation theorem (61), the cumulant generating function obeys the following symmetry relation

[TABLE]

Now, the mean currents and the diffusivities can be obtained by taking the successive derivatives of the generating function (81) with respect to the counting parameters:

[TABLE]

for $\alpha,\beta=C,B$. Besides, we may expand the mean currents in power series of the affinities as

[TABLE]

in terms of the response coefficients defined by

[TABLE]

The coefficients $L_{\alpha,\beta}$ characterize the linear response properties and the coefficients $M_{\alpha,\beta\gamma}$ the nonlinear response properties of the currents at second order in the affinities. The coefficients of higher orders can also be introduced AG07JSM ; BG18 .

If we take the derivatives of the symmetry relation Eq. (82) with respect to $\lambda_{\alpha}$ and $A_{\beta}$, and set ${\boldsymbol{\lambda}}={\bf 0}$ and ${\bf A}={\bf 0}$, we obtain the fluctuation-dissipation relations

[TABLE]

and the Onsager reciprocal relations

[TABLE]

as a consequence of the symmetry $D_{\alpha\beta}=D_{\beta\alpha}$ resulting from the definition (84) of the diffusivities.

If we take a further derivative of the symmetry relation (82) with respect to $A_{\gamma}$ before setting ${\boldsymbol{\lambda}}={\bf 0}$ and ${\bf A}={\bf 0}$, we find that

[TABLE]

giving the nonlinear response coefficient $M_{\alpha,\beta\gamma}$ in terms of the first responses of the diffusivities around equilibrium. The relations (90) as well as the Onsager reciprocal relations (89) find their origin in the microreversibility underlying the fluctuation theorem for currents S92 ; AGMT09 ; EHM09 ; CHT11 ; S12 ; G13 .

V.2 Numerical test of the linear transport properties

In this subsection, we focus on the numerical test of the fluctuation-dissipation relations (88) and the Onsager reciprocal relation (89) for $\alpha,\beta=C,B$. Here, we use the methods given in Appendix D for the numerical evaluation of derivatives and their error analysis.

The evaluation of the linear response coefficients relies on the determination of the mean currents as a function of the affinities. To achieve this evaluation, we have computed the mean currents for several values of the affinities, as shown in Fig. 6. We have used the Lagrange interpolation method to obtain one-variable polynomials approximating $J_{C}(A_{C},A_{B}=0)$, $J_{C}(A_{C}=0,A_{B})$, $J_{B}(A_{C},A_{B}=0)$, and $J_{B}(A_{C}=0,A_{B})$ based on the numerical data plotted in Fig. 6. Subsequently, the linear response coefficients can be computed by taking the first partial derivatives of the Lagrange polynomials at the equilibrium point $(A_{C}=0,A_{B}=0)$. Their numerical values are given in the first column of Table 5. This computation already confirms that the Onsager reciprocal relation $L_{C,B}=L_{B,C}$ is satisfied within the numerical accuracy.

Furthermore, the equilibrium values of the diffusivities are computed using Eqs. (53)-(55), giving the values in the second column of Table 5. The difference between the linear response coefficients and the diffusivities are reported in the third column of Table 5, showing that the fluctuation-dissipation relations (88) are also satisfied within the numerical accuracy.

V.3 Numerical test of the nonlinear transport properties

The numerical values of the charge currents $J_{C}$ and $J_{B}$ are computed for different values of the affinities $A_{C}$ and $A_{B}$ in order to construct the two-variable functions $J_{C}(A_{C},A_{B})$ and $J_{B}(A_{C},A_{B})$ using two-dimensional Lagrange interpolations, as shown in Fig. 7. The values of second derivatives at the equilibrium point $(A_{C}=0,A_{B}=0)$,

[TABLE]

are thus numerically evaluated in order to determine the nonlinear response coefficients $M_{\alpha,\beta\gamma}$, using the numerical method explained in Appendix D. On the other hand, the diffusivities $D_{\alpha\beta}$ are again computed using Eqs. (53)-(55), but for the transistor driven away from equilibrium. They are plotted in Fig. 8 as functions of the affinities. Therefore, the derivatives of the diffusivities with respect to the affinities

[TABLE]

can also be evaluated numerically at the equilibrium point $(A_{C}=0,A_{B}=0)$. The results for the quantities $M_{\alpha,\beta\gamma}$ and $R_{\alpha\beta,\gamma}$ are given in Table 6 where we calculate the differences, $M_{\alpha,\beta\gamma}-R_{\alpha\beta,\gamma}-R_{\alpha\gamma,\beta}$, testing the validity of the prediction (90) of the fluctuation theorem beyond the linear transport properties. We see that these differences are smaller than the numerical errors in agreement with the predictions.

VI Conclusion and Perspectives

Using a spatially extended stochastic description of charge transport in bipolar $n$-$p$-$n$ junction transistors, we have shown in this paper that a fluctuation theorem holds for the two electric currents that are coupled together in the double junction of the transistor. We have also shown that, as a corollary of the fluctuation theorem for the currents, nonlinear transport generalizations of the fluctuation-dissipation and Onsager reciprocal relations are satisfied in the transistor. In particular, we have verified in detail that the second-order nonlinear response coefficients of the currents are related to the first-order responses of the diffusivities, as predicted by theory AG04 ; AG07JSM ; BG18 .

These results are based on stochastic partial differential equations describing the diffusion of electrons and holes, as well as their generation and recombination. These stochastic diffusion-reaction equations are coupled to the Poisson equation for the electric potential and they obey local detailed balance. The scheme is consistent with the laws of electricity, thermodynamics, and microreversibility. The stochastic process is driven out of equilibrium by boundary conditions due to the voltages applied to the reservoirs in contact with the three ports of the transistor. In this case, the transistor is the stage of a nonequilibrium steady state, manifesting highly nonlinear transport properties. The key point raised in this paper is that, besides their amazing technological importance, transistors can be used to address the fundamental issue of microreversibility in nonequilibrium statistical physics.

The one-variable fluctuation theorem has already been experimentally investigated in linear $RC$ electric circuits GC05 ; JGC08 . Our previous paper GG18 has shown that the one-variable fluctuation theorem can be studied in nonlinear devices such as diodes. In transistors, the experimental test of the two-variable fluctuation theorem can also be envisaged, either by the direct measurement of current fluctuations, or by testing its consequences, namely, the time-reversal symmetry relations generalizing the fluctuation-dissipation and Onsager reciprocal relations to the nonlinear transport properties. Such tests would require accurate noise measurements with large enough statistics. In this way, these symmetry relations, finding their origins in the fundamental law of microreversibility, could be tested experimentally in common devices of modern technology.

Acknowledgments

The authors thank Sergio Ciliberto for stimulating discussions. Financial support from the China Scholarship Council under the Grant No. 201606950037, the Université libre de Bruxelles (ULB), and the Fonds de la Recherche Scientifique - FNRS under the Grant PDR T.0094.16 for the project “SYMSTATPHYS” is acknowledged.

Appendix A Discretized Markov jump process

To describe the BJT by a Markov jump process, the system is spatially discretized into cells of volume $\Omega$, each containing some numbers $N_{i}$ and $P_{i}$ of electrons and holes, respectively. These numbers are supposed to change in time because of random transitions at rates to be specified here below. The Markov jump process is fully defined by these transition rates and the master equation ruling the time evolution of the probability that the cells contain given numbers of particles. In the continuum limit, the Markov jump process leads to the stochastic reaction-diffusion equations (11)-(27), as shown in Appendix B. This method is similar to the one used in Refs. AG09 ; GG18 ; G05 .

A.1 Master equation of the process

At any time, the state of the discretized BJT is fully characterized by the electron numbers ${\bf N}=(N_{i})_{i=1}^{L}$ and hole numbers ${\bf P}=(P_{i})_{i=1}^{L}$ in all the cells. The time evolution of these numbers is ruled by a Markov jump process corresponding to the following network:

[TABLE]

On the left-hand side, the Collector C is a reservoir of electron and holes where their numbers $\bar{N}_{C}$ and $\bar{P}_{C}$ take fixed values. On the right-hand side, it is the Emitter E that fixes the values of $\bar{N}_{E}$ and $\bar{P}_{E}$. In the middle, similar transitions happen with the Base B, fixing the values of $\bar{N}_{B}$ and $\bar{P}_{B}$. The transitions with the rates $W_{i}^{(\pm N)}$ describe the diffusive transfers of electrons between the cells, and those with the rates $W_{i}^{(\pm P)}$ the diffusive transfers of holes. The transitions with the rates $W_{i}^{(\pm)}$ describe the generation and recombination of electron-hole pairs, respectively.

The probability ${\cal P}(N_{1},\dots,N_{L},P_{1},\dots,P_{L},t)$ to find the system in a certain state is thus governed by the master equation

[TABLE]

with the transition rates given by

[TABLE]

For electron, the transition rates at the boundaries are given by

[TABLE]

and similar expressions for holes. We note that, in the network shown above, the cell $i=m$ is the only one in contact with the Base, in which case the sum $\sum_{iB}$ in Eq. (113) has the sole term $i=m$.

$U$ is the total electrostatic energy stored in the BJT and $\Delta U$ is the energy difference associated with the change of the BJT state. $\psi(\Delta U)$ is a function defined by

[TABLE]

which satisfies the local detailed balance condition

[TABLE]

A.2 Discretized Poisson equation

The Poisson equation is replaced by its discretized version

[TABLE]

with the boundary conditions $\phi_{0}=\phi_{C}$ and $\phi_{L+1}=\phi_{E}$ at two ends of BJT, and the symbol $\chi_{iB}=1$ if the $i^{\rm th}$ cell is in contact with the Base and $\chi_{iB}=0$ otherwise. This linear system should be solved after every electron or hole transfer between cells. We suppose that the electric potential $\phi_{B}$ of the Base is set on both sides of the chain in the transverse $y$-direction, in order to get a symmetric geometry.

The electrostatic energy is given by

[TABLE]

where the electric potential

[TABLE]

obeys the discretized Poisson equation

[TABLE]

with the symmetric matrix

[TABLE]

expressed in terms of the Kronecker symbol such that $\delta_{i,j}=1$ if $i=j$ and $\delta_{i,j}=0$ otherwise, the coefficients

[TABLE]

and

[TABLE]

The change of electrostatic energy during the transfer of an electron of charge $-e$ from the $i^{\rm th}$ to the $(i+1)^{\rm th}$ cell is given by

[TABLE]

where

[TABLE]

so that

[TABLE]

A similar expression holds for hole transfers since they have the charge $+e$,

[TABLE]

For electron transfers at the boundary, we have

[TABLE]

and similar expressions for holes.

Appendix B Langevin stochastic process

In the limit where $N_{i}\gg 1$ and $P_{i}\gg 1$, the Markov jump process described here above can be replaced by a Langevin stochastic process GG18 ; G05 , which is ruled by another master equation obtained by expanding the operators $\exp(\pm\partial_{X})$ up to second order in the partial derivatives $\partial_{X}$ in Eq. (113). In this way, we find that the corresponding probability density $\mathscr{P}$ obeys the following Fokker-Planck equation:

[TABLE]

This shows that the variables $N_{i}$ and $P_{i}$ obey stochastic differential equations of Langevin type:

[TABLE]

with the following fluxes and reaction rates:

[TABLE]

expressed in terms of the Gaussian white noises:

[TABLE]

These Langevin stochastic equations are numerically implemented by discretizing time into equal intervals $\Delta t$ and replacing the white noises by independent identically distributed Gaussian random variables. The stochastic partial differential equations (11)-(27) are recovered in the continuum limit GG18 .

Appendix C Coarse-grained Markov jump process

For the simple coarse-grained model (80), the joint probability distribution $P(Z_{C},Z_{B},t)$ to observe the charge transfers $Z_{C}$ and $Z_{B}$ during the time interval $[0,\,t]$ evolves according to the following master equation

[TABLE]

According to the central limit theorem, the joint probability distribution $P(Z_{C},Z_{B},t)$ after a long enough time interval $[0,\,t]$ becomes Gaussian of the following form,

[TABLE]

with the vectorial and matricial notations

[TABLE]

and T denoting the transpose. The mean charge currents ${\bf J}$ and the diffusivities ${\boldsymbol{\mathsf{D}}}$ can be numerically evaluated through

[TABLE]

where $\langle\cdot\rangle$ denotes the statistical average over the data sample.

For this model, the mean currents and the diffusivities can be expressed in terms of the transition rates of the master equation (161) according to the following relations:

[TABLE]

By local detailed balance, the affinities are given by

[TABLE]

The natural condition

[TABLE]

leads to

[TABLE]

Eqs. (165)-(169) and (174) form a set of six nonlinear equations that can be solved numerically with the Newton-Raphson method to find the six transition rates $\{W_{kl}\}_{k,l=C,B,E}$. Thereafter, the affinities are readily evaluated by Eqs. (170)-(172). Taking the Emitter as the reference reservoir, we may more shortly write $A_{CE}$ as $A_{C}$, and $A_{BE}$ as $A_{B}$.

We note that these considerations lead to the Ebers-Moll transport model of bipolar junction transistors EM54 ; SS04 if we assume that $W_{CB}=J_{s}/\beta_{R}$, $W_{EB}=J_{s}/\beta_{F}$, and $W_{EC}=J_{s}\exp(\beta eV_{BC})$, where $J_{s}$ is the reverse saturation current, $\beta_{R}$ the reverse common emitter current gain, and $\beta_{F}$ the forward common emitter current gain, in addition to the local detailed balance conditions $W_{BC}=W_{CB}\exp(\beta eV_{BC})$ and $W_{BE}=W_{EB}\exp(\beta eV_{BE})$ given by Eqs. (171) and (172). The well-known expressions for the mean currents of this model (e.g., given Ref. SS04 pp. 387-389) are thus recovered from Eqs. (165) and (166) by using Eq. (174).

Appendix D Numerical differentiation and error analysis

The differentiation can be approximated by numerical differences using several points AS72 . Given the values of the one-variable function $f(x)$ at the five equispaced points $-2h$, $-h$, [math], $h$, $2h$, we have the following centered-difference formulae

[TABLE]

respectively giving the first- and second-order derivatives up to numerical errors of $O(h^{4})$. These two difference formulae can be obtained using the Lagrange polynomial

[TABLE]

that interpolates the five points at $x_{i}=-2h,\,-h,\,0,\,h,\,2h$. Here, it is easy to obtain Lagrange polynomial corresponding to the two-variable function $f(x,y)$ using points distributed on a grid

[TABLE]

The mixed second derivative of $f(x,y)$ at the point $(0,0)$ can be approximated by the midpoint formula

[TABLE]

which is accurate up to $O(h_{1}^{2}h_{2}^{2})$.

Apart from the numerical error itself, another source of errors comes from the statistical evaluation of the function at the different points. Suppose that the variances of the numerical values of the function are denoted as $\sigma^{2}\left[f(x_{i})\right]$ and $\sigma^{2}\left[f(x_{i},y_{j})\right]$, then the mean square errors on the derivative (175) can be evaluated as

[TABLE]

up to a correction of $O(h^{8})$ coming from the error in the numerical differentiation. Similar expressions hold for the mean square errors on the other derivatives (176) and (179).

Given the random sample $\{X_{1},\dots,X_{n}\}$ of size $n$ from a Gaussian distribution of mean value $\mu$ and variance $\sigma^{2}$, the sample average is defined as $\langle X\rangle=(1/n)\sum_{i=1}^{n}X_{i}$, having the expected value equal to $\mu$. The sample average $\langle X\rangle$ has the mean square error ${\rm MSE}\left(\langle X\rangle\right)=\sigma^{2}/n$. The unbiased sample variance $S_{n-1}^{2}=\sum_{i=1}^{n}(X-\langle X\rangle)^{2}/(n-1)$ has the expected value $\sigma^{2}$ and its mean square error is equal to ${\rm MSE}\left(S_{n-1}^{2}\right)=2\sigma^{4}/(n-1)$. If we define the average current $J=\langle X\rangle/t$ and diffusivity $D=S_{n-1}^{2}/(2t)$, their mean square errors can thus be estimated as

[TABLE]

The procedure used to estimate the error on the numerical computation of the affinities $A_{C}=\ln\left(W_{CE}/W_{EC}\right)$ and $A_{B}=\ln\left(W_{BE}/W_{EB}\right)$ by the method of Appendix C is the following. The expressions of the affinities are differentiated with respect to the mean values of the currents and diffusivities to obtain linear approximations such as

[TABLE]

in terms of some coefficients $a$, $b$, $c$, $d$, and $e$, which are related to the rates. Accordingly, the mean square error is estimated as

[TABLE]

and similarly for the error on $A_{B}$.