A drift-diffusion model based on the fractional exclusion statistics

We propose a drift-diffusion model for systems which obey fractional exclusion statistics (FES), in a framework where the species include classical degrees of freedom such as positions. The transition rates are calculated and the relation between the step and acceptance probabilities on one hand and the diffusion and drift processes on the other hand are established. A Monte Carlo scheme is implemented on a prototypical double-junction system of particles with screened Coulomb interactions. In our approach the properties of interacting quantum gases are locally included using the FES methodology. The model is suitable to describe transient as well as stationary regimes.


Introduction
During the past decades the drift-diffusion models have been invaluable tools for describing the transport in semiconductor devices and the physics of cold plasmas and gas discharges. Efficient solvers based on Boltzmann transport equation implementing Monte Carlo schemes have been devised. However, in recent years, the quantum effects in modern semiconductor devices became increasingly important.
Solving microscopic quantum models-although they are most accurate-is becoming an increasingly difficult task as the size of the analyzed system is getting larger in the attempt to reach the experimentally attainable length scales. Therefore hybrid microscopic-macroscopic models are usually employed, which combine a microscopic description for short length scales with macroscopic fluid-type models, e.g. quantum drift diffusion model (QDD) [1], quantum corrected drift diffusion model (QCDD) [2], Schrödinger-Poisson-Drift-Diffusion model (SPDD) [3].
We propose here a method to locally incorporate in a transport model the behavior of interacting quantum gases, using fractional exclusion statistics (FES) [4][5][6]. The method was shown to be a tool for the study of the thermodynamic properties of interacting particle systems, regarded as ideal FES gases. The transition rates for homogeneous FES gases were established [7][8][9], which opened the possibility of analyzing the non-equilibrium dynamics of interacting bosonic and fermionic systems via Monte Carlo simulations [9]. The extension of this kinetic model to non-homogeneous FES systems with mutual exclusion parameters was realized in Refs. [10][11][12]. We make here a step further and introduce a transport model based on FES, using Monte Carlo simulations.

Model and formalism
The system under investigation is a generic device made up by an active region and contacts. The transport properties are dominated by hopping between states localized at the sites r i . The contacts have the role of particle reservoirs which fix the chemical potentials according to the applied bias. The key elements of the active region, which determine properties of the system, are the local density of states σ(r, ) and the interacting potentials which, in general, may be position dependent. The Hamiltonian of the system is [11,12] where r I is the energy and n r I is the occupation number of the site I; the total particle number is N = I n r I . The quasiparticle energies are defined as˜ In order to apply the FES formalism we work in the quasicontinuous limit and we define the species as in Refs. [11,12] by dividing the real space and the quasiparticle energy axis into elementary volumes and segments δr and δ . If the space is s-dimensional (sD), then the species will be the (s + 1)D volumes δr ξ × δ i , where by ξ (Greek) and i (Latin) we index the volumes and energy intervals, respectively. Each species (ξ, i) contains N ξi particles and G ξi (for bosons) or T ξi (for fermions) single-particle states. Then the number of micro-configurations corresponding to a distribution {G ξi , N ξi } (bosons) or {T ξi , N ξi } (fermions) of states and particles is The FES parameters are defined by the change of the number of states in one species (ξi) due to the change of the number of particles in any (other or the same) species (ηj), which is δG ξi = −α ξi,ηj δN ηj for bosons and δT ξi = −α ξi,ηj δN ηj for fermions. For the Hamiltonian (1) and our choice of quasiparticle energies the FES parameters are [11,12]. Note that if the σ ξ ( i ) does not depend on the energy, the FES parameters simplify to α ξi;ηj = δ ij σ ξ ( i )V (|r ξ − r η |).
In FES the description of bosons and fermions may be unified by the redefinition of the FES parameters and of the number of states obtaining what we may call bosonic or fermionic descriptions. In the bosonic description-the one that we shall use in the following-the number of states is G ξi for bosons and T ξi − N ξi + 1 for fermions. Similarly, the FES parameters are α ξi,ηj for bosons and α ξi,ηj + δ ξη δ ij for fermions. If we make these transformations then the number of microconfigurations is given by W B ξi for both, bosons and fermions. If we assume that the system is in thermal equilibrium, then the probability to find the system in the state Υ, defined by the set {N ξi , G ξi } ξi , is where Z is the partition function and E Υ = ξ,i N ξi ξi is the total energy-we assume that all the particles that belong to one species (ξi) have the same energy ξi . Two states of the FES system, Υ and Φ, are considered neighbors if they differ by one particle jump from one species (ξ, i) to another (η, j), (where (ζ, k) = (ξ, i), (η, j)). Upon extracting one particle from species (ξ, i) and inserting it into species (η, j), the number of available states may change in any species (ζk) by G ζk = G ζk + α ζk,ξi − α ζk,ηj . Introducing a discrete-time reversible Markov chain one can employ the detailed balance equation, where Γ ΦΥ , Γ ΥΦ are the transition rates from Υ to Φ and vice-versa, respectively. From the detailed balance equation we obtain the ratio between the transition rates where, in last step, we used a first order expansion for log(p eq Φ ) around Υ. In the thermodynamic limit, i.e. for large N ξi and G ξi , we apply the Stirling formula to calculate log(p eq Υ ) and by using the identity α ξi,ηj = −∂G ηj /∂N ξi we find where (ζ, k) = (ξ, i), (η, j) and a similar expression for ∂(log p eq Υ )/∂N ηj by interchanging (ξ, i) and (η, j). Defining the number of available states in the absence of particles [5] by G 0 ξi = G ξi + η,j α ξi,ηj (N ηj − δ ξi,ηj ) and the occupation numbers as n ξi = N ξi /G 0 ξi , from Eq. (7) we find the transition rate for the Υ → Φ process where A ξi = η,j α ξi,ηj n ηj . One should note that in the case of a homogeneous system with diagonal statistical parameters, α ξi,ηj = αδ ξη δ ij , one recovers the transition rates reported in Ref. [9]: The transition rate in equation (9) represents the step probability for one particle in species (ξ, i) to select (η, j) as target species. Once the destination species is obtained, the move is accepted with a probability proportional with the Boltzmann factor, i.e. ∼ min[1, exp(−(E Φ − E Υ )/k B T )]. The step probability depends on the occupancies of both source and target species and can be associated with a diffusion process. The acceptance probability depends on the energy diference between the two species and, in the case of external fields, is influenced and may be associated with the drift process.

Application -Transport in a double junction system with screened Coulomb interactions
We apply the drift-diffusion model on a one-dimensional test-case system. We consider a two terminal device partitioned into three regions: the left contact, the active region and the right contact. The particles are interacting via a screened Coulomb interaction, which is parametrized by: where r is the distance between the particles, R 0 is a cut-off radius which removes the singularity at the origin, V 0 ≡ V (R 0 , λ), and λ is the screening length. We assume the interaction inside the contacts is strongly screened, while in the active region we consider both limiting situations which correspond to short and large screening lengths.
Our prototypical system consists of N x ×N species, with N x = 16 and N = 50, and constant density of states. The species are equally spaced in the interval [−L x /2, L x /2], where L x = 1. The length of each of the two contact regions is L x /4 and it is assumed they have a short screening length λ C = L x /N x , while for the active region, of length L x /2, we compare two screening lengths, namely, λ AR = λ C and λ AR → ∞. In the following, all energies are scaled with respect to E F , the Fermi energy of the non-interacting system. The chemical potentials in the left and the right contacts are set to μ L = 2.5 and μ R = 2, respectively, their difference reflecting the applied bias. The cutoff radius is set to R 0 = 0.05 and V 0 = 8. The temperature of the system is k B T = 0.25.
The species ξ = 0 and ξ = N x − 1 have fixed populations, as imposed by the chemical potentials in the contacts. Whenever a particle enters / leaves these species, the particle is removed / added back, respectively.
We first performed a number of 10 10 Monte Carlo steps in order to equilibrate the system, and then using an equal number of steps the thermal averages are performed. The particle densities, ρ ξ = i n ξi , are presented in Fig. 1 (a) and (b) for the two cases -short and long screening lengths in the active region. For the homogeneous system, λ AR = λ C , one obtains a linear dependence of the particle density, as particles are moving from the left to the right contact. By contrast, for pure Coulombian interaction one observes a decrease in the particle density in the active region. The depletion region is skewed by the finite potential difference between the two contacts and a spike-like feature appears in the particle density in the vecinity of the right contact interface, in the species ξ = 11. However, the regions of the non-interacting contacts embedded in the active region still retain a linear dependence in the particle distribution. Correspondingly, Figs. 1 (c) and (d) show the populations in the two contacts, as well as for the species at three other sites. For ξ = 10 one obtains a minimum in the particle density for λ AR → ∞ and the populations of the quasiparticle energy levels are comparatively smaller. This charge depletion effect is due to the unscreened Coulombian repulsion and it may be found, e.g., in a n++ / intrinsic / n++ semiconductor structure.
Within the proposed method one may extract the current of particles flowing between terminals, from the number of particles inserted or extracted from the boundary species, ξ = 0 and N x − 1, respectively. Although in stationary regime the two currents are equal up to some finite system fluctuations, generally, in non-stationary regime, they may differ and one can describe in real time transient processes, e.g. charge accumulation or depletion.

Conclusions
We proposed a drift-diffusion model for systems which obey fractional exclusion statistics. We determined the transition rates, pointing out the relation between the step and acceptance probabilities and the diffusion and drift processes.
The model is employed on a prototypical double-junction system of particles with screened Coulomb interactions.
These types of systems are relevant for device applications in semiconductor physics, e.g. field effect transistors. Depletion effects are analyzed with respect to the screening length in the device active region.
Our approach takes advantage of a local description of the system in terms of ideal FES gases. This enables us to handle larger systems, in the range where exact quantum mechanical calculations are difficult to be performed and yet incorporating locally the quantum behavior of the interacting particle systems by means of the FES methodology.