Two spinorial drift-diffusion models for quantum electron transport in graphene

TL;DR

This paper derives and analyzes two quantum drift-diffusion models for electron transport in graphene that incorporate spin effects, providing theoretical properties and numerical simulations to understand their behavior.

Contribution

It introduces two new spinorial drift-diffusion models for graphene based on a spinorial Wigner equation, including their derivation, mathematical analysis, and numerical validation.

Findings

Proved global existence of weak solutions for one model

Established entropy-dissipation and decay properties

Performed numerical simulations demonstrating model behavior

Abstract

Two drift-diffusion models for the quantum transport of electrons in graphene, which account for the spin degree of freedom, are derived from a spinorial Wigner equation with relaxation-time or mass- and spin-conserving matrix collision operators using a Chapman-Enskog expansion around the thermal equilibrium. Explicit models are computed by assuming that both the semiclassical parameter and the scaled Fermi energy are sufficiently small. For one of the models, the global existence of weak solutions, entropy-dissipation properties, and the exponential long-time decay of the spin vector are proved. Finally, numerical simulations of a one-dimensional ballistic diode using both models are presented, showing the temporal behavior of the particle density and the components of the spin vector.

Tables

Table 1. Table 1. Units used for the numerical simulations.

space unit	${10}^{-7}$ m
time unit	$0.5\times {10}^{-13}$ s
voltage unit	$1.25\times {10}^{-2}$ V
particle density unit	${10}^{17}$ m-2
current density unit	$2\times {10}^{23}$ m-1 s-1

Equations

$$H_0=-i\hslash v_F\left({\sigma }_1\frac{\partial }{\partial x_1}+\sigma \frac{\partial }{\partial x_2}\right),$$

$$H=H_0-\left(\frac{{\hslash }^2}{2m}\Delta \right){\sigma }_0,$$

$${\partial }_t{n}_0-\text{div}{J}_0=0,J_0=\nabla n_0+n_0\nabla V,$$

$${\partial }_t{n}_j-\text{div}{J}_j=F_j,J_j=A_0(\vec{n}/n_0)\nabla \vec{n}+\vec{n}\otimes \nabla V+B_0(n_0,\vec{n}),j=1,2,3,$$

$${\partial }_t{n}_0-\text{div}{J}_0=0,J_0=A_1(\nabla n_0+n_0\nabla V)+B_1\cdot (\nabla \vec{n}+\vec{n}\otimes \nabla V)+C_1(n_0,\vec{n}),$$

$${\partial }_t{n}_j-\text{div}{J}_j=G_j,J_j=A_2(\nabla n_0+n_0\nabla V)+B_2\cdot (\nabla \vec{n}+\vec{n}\otimes \nabla V)+C_2(n_0,\vec{n}),$$

$${\partial }_{t}u-\Delta u-\nabla (\log n_0+V)\cdot \nabla u=2G[\vec{n}/n_0],$$

$${\sigma }_0=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),{\sigma }_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\sigma }_2=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),{\sigma }_3=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)$$

$${\partial }_t{w}_0+\left(\frac{\vec{p}}{m}\cdot \vec{\nabla }\right)w_0+v_F\vec{\nabla }\cdot \vec{w}+{\theta }_{\hslash }[V]w_0$$

$${\partial }_t\vec{w}+\left(\frac{\vec{p}}{m}\cdot \vec{\nabla }\right)\vec{w}+v_F\left(\vec{\nabla }{w}_0+\frac{2}{\hslash }\vec{w}\wedge \vec{p}\right)+{\theta }_{\hslash }[V]\vec{w}$$

$$({\theta }_{\hslash }[V]w)(x,p)=\frac{i}{\hslash }\frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}{\int }_{{\mathbb{R}}^2}\delta V(x,\xi )w(x,p^{\prime })e^{-i(p-p^{\prime })\cdot \xi }d\xi d{p}^{\prime },$$

$$\delta V(x,\xi )=V\left(x+\frac{\hslash }{2}\xi \right)-V\left(x-\frac{\hslash }{2}\xi \right).$$

$$x\to \hat{x},t\to \hat{t}t,p\to \hat{p}p,V\to \hat{V}V,$$

$$\frac{2v_F\hat{p}}{\hslash }=\frac{\hat{V}}{\hat{x}\hat{p}},\frac{2\hat{p}{v}_F{\tau }_c}{\hslash }=\frac{\hslash }{2\hat{p}{v}_F\hat{t}},\hat{p}=\sqrt{m{k}_{B}T}.$$

$$\frac{\hat{x}}{L}=\frac{1}{2}\frac{\hat{V}(\hslash /{\tau }_c)}{({\hat{p}}^2/m)(m{v}_F^2)}.$$

$$\frac{\hat{t}}{{\tau }_c}=\frac{1}{4}\frac{(\hslash /{\tau }_c{)}^2}{({\hat{p}}^2/m)(m{v}_F^2)}$$

$$\epsilon =\frac{\hslash }{\hat{x}\hat{p}},\tau =\frac{2\hat{p}{v}_F{\tau }_c}{\hslash },c=\sqrt{\frac{m{v}_F^2}{k_{B}T}}.$$

$$\gamma :=\frac{c}{\epsilon }=O(1)\text{mbox}as\epsilon \to 0.$$

$$\tau {\partial }_t{w}_0+\frac{1}{2\gamma }(\vec{p}\cdot \vec{\nabla })w_0+\frac{\epsilon }{2}\vec{\nabla }\cdot \vec{w}+{\theta }_{\epsilon }[V]w_0$$

$$\tau {\partial }_t\vec{w}+\frac{1}{2\gamma }(\vec{p}\cdot \vec{\nabla })\vec{w}+\frac{\epsilon }{2}\vec{\nabla }{w}_0+\vec{w}\wedge \vec{p}+{\theta }_{\epsilon }[V]\vec{w}$$

$$A[w]={\int }_{{\mathbb{R}}^2}{\int }_{{\mathbb{R}}^2}\text{mbox}tr\left(w\left(\text{mbox}Log(w)-1+\frac{h(p)}{k_{B}T}\right)\right)dxdp,$$

$$(\text{mbox}{Op}_{\epsilon }(\gamma )\psi )(x)=\frac{1}{(2\pi \hslash {)}^2}{\int }_{{\mathbb{R}}^2}{\int }_{{\mathbb{R}}^2}\gamma \left(\frac{x+y}{2},p\right)\psi (y)e^{i(x-y)\cdot p/\hslash }dydp,$$

$$h(p)=\frac{|p{|}^2}{2m}{\sigma }_0+v_F(p_1{\sigma }_1+p_2{\sigma }_2)$$

$$A\big{[}g[n_{0},\vec{n}]\big{]}=\min_{w}\left\{A[w]:\int_{\mathbb{R}^{2}}w_{0}dx=n_{0},\quad\int_{\mathbb{R}^{2}}\vec{w}dx=\vec{n}\right\},$$

$$g[n_0,\vec{n}]=\text{mbox}Exp(-h_{A,B}),h_{A,B}=\left(\frac{|p{|}^2}{2}+A\right){\sigma }_0+(c\vec{p}+\vec{B})\cdot \vec{\sigma },$$

$${\int }_{{\mathbb{R}}^2}g[n_0,\vec{n}](x,p,t)dp=n_0(x,t),{\int }_{{\mathbb{R}}^2}\vec{g}[n_0,\vec{n}](x,p,t)dp=\vec{n}(x,t),$$

$$g[n_0,\vec{n}]=\text{mbox}Exp(a+\epsilon b),a=-\left(\frac{|p{|}^2}{2}+A\right){\sigma }_0-\vec{B}\cdot \vec{\sigma },b=-\gamma \vec{p}\cdot \vec{\sigma }.$$

$$g^{(0)}$$

$$g^{(1)}$$

$$\displaystyle\phantom{xx}{}-\Bigg{[}\left(\left(\cosh|\vec{B}|-\frac{\sinh|\vec{B}|}{|\vec{B}|}\right)\frac{\vec{B}\otimes\vec{B}}{|\vec{B}|^{2}}+\frac{\sinh|\vec{B}|}{|\vec{B}|}I\right)\vec{p}$$

$$\displaystyle\phantom{xx}{}+\left(\cosh|\vec{B}|-\frac{\sinh|\vec{B}|}{|\vec{B}|}\right)\frac{((\vec{p}\cdot\vec{\nabla}_{x})\vec{B})\wedge\vec{B}}{2\gamma|\vec{B}|^{2}}\Bigg{]}\cdot\vec{\sigma}\Bigg{\}},$$

$$e^{-A}=\frac{1}{2\pi }\sqrt{n_0^2-|\vec{n}{|}^2}+O({\epsilon }^2),\vec{B}=-\frac{\vec{n}}{|\vec{n}|}\log \sqrt{\frac{n_0+|\vec{n}|}{n_0-|\vec{n}|}}+O({\epsilon }^2),$$

$$\lambda :=\frac{c}{\tau }=\frac{\epsilon \gamma }{\tau }=O(1)\text{mbox}as\tau \to 0.$$

$$f_0=-\frac{1}{2\gamma }(\vec{p}\cdot \vec{\nabla })g_0+\vec{\nabla }V\cdot {\vec{\nabla }}_p{g}_0,\vec{f}=-\frac{1}{2\gamma }(\vec{p}\cdot \vec{\nabla })\vec{g}+\vec{\nabla }V\cdot {\vec{\nabla }}_p\vec{g}-\vec{g}\wedge \vec{p}.$$

Full text

The second author acknowledges partial support from the Austrian Science Fund (FWF), grants P20214, P22108, and I395, and from the Austrian-French Project of the Austrian Exchange Service (ÖAD).

1. Introduction

Graphene is a new semiconductor material, which is a subject of great interest for nanoscale electronic applications. The reason for this interest is due to the very remarkable properties of graphene, such as the high electron mobility and long coherence length. Therefore, graphene is a promising candidate for the construction of a new generation of electronic devices with far better performances than current silicon devices [8]. Potential applications include, for instance, spin field-effect transistors [12, 24], extremely sensitive gas sensors [23], one-electron graphene transistors [19], and graphene spin transistors [3].

Physically, graphene is a two-dimensional semiconductor with a zero-width band gap, consisting of a single layer of carbon atoms arranged in a honeycomb lattice. In the energy spectrum, the valence band intersects the conduction band at some isolated points, called the Dirac points. Around these points, quasiparticles in graphene exhibit the linear dispersion relation $E=v_{F}|p|$, where $p$ denotes the crystal momentum and $v_{F}\approx 10^{6}$ m/s is the Fermi velocity [16]. This energy spectrum resembles the Dirac spectrum for massless relativistic particles, $E=c|p|$, where $c$ is the speed of light. Hence, the Fermi velocity $v_{F}\approx c/300$ takes the role of the speed of light. The system Hamiltonian can be approximated near a Dirac point, for low energies and in the absence of a potential, by the Dirac-like operator

[TABLE]

where $\hbar$ is the reduced Planck constant, and $\sigma_{1}$ and $\sigma_{2}$ are the Pauli matrices (see (3)).

In order to understand and predict the charge carrier transport in graphene, transport models, which incorporate the spin degree of freedom, have to be devised. Theoretical models for spin-polarized transport involve fluid-type drift-diffusion equations, kinetic transport equations, and Monte-Carlo simulation schemes; see the references in [18]. A hierarchy of fluiddynamic spin models was derived from a spinor Boltzmann transport equation in [2]. Suitable matrix collision operators were suggested and analyzed in [20]. Drift-diffusion models for spin transport were considered in several works; see, e.g., [1, 6, 22]. A mathematical analysis of spin drift-diffusion systems for the band densities is given in [10].

Fluiddynamic equations provide a compromise between physical accuracy and numerical efficiency. Another advantage is that they contain already the physically interesting quantities, such as the particle density, momentum, and spin densities, whereas other models usually involve variables which do not have an immediate physical interpretation, like wavefunctions, density operators, and Wigner distributions. In the latter case, further computations have to be made to obtain the quantities of physical interest.

In this work, we address the quantum kinetic and diffusion level of spin-polarized transport in graphene. More precisely, starting from a spinorial Wigner equation, we aim to derive via a moment method and a Chapman-Enskog expansion macroscopic drift-diffusion models for the particle density and spin vector. Furthermore, we prove the global existence of weak solutions to one of these models and we illustrate the behavior of the solutions in a ballistic diode by numerical experiments..

We note that there are only very few articles concerned with kinetic or macroscopic transport models for graphene. In the physics literature, the focus is on transport properties such as the carrier mobility [11], charged impurity and phonon scattering [4], and Klein tunneling [17]. Wigner models were investigated in [15]. Starting from a Wigner equation, hydrodynamic spin models were derived in [27], and the work [26] is concerned with the derivation of drift-diffusion models for the band densities. In contrast, we will work in the present paper with all components of the spin vector. Furthermore, we provide a mathematical analysis of one of the models and numerical simulations of both models.

In the following, we describe our approach and the main results. ¿From the unique features of graphene follow that Fermi-Dirac statistics would be more suitable to describe quantum transport in the material than Maxwell-Boltzmann statistics, since the energy spectrum of the Hamiltonian (1) is not bounded from below. We overcome this problem by modifying the Hamiltonian $H_{0}$. In fact, we assume that the system Hamiltonian is approximated by the following operator which is bounded from below:

[TABLE]

where $m>0$ is a parameter with the dimension of a mass and $\sigma_{0}$ is the unit matrix in $\mathbb{R}^{2\times 2}$ (see (3)). This is not a very restrictive assumption since the operator (1) is itself only an approximation of the correct system Hamiltonian, valid for small values of the momentum $|p|$.

Starting from the Hamiltonian $H+V\sigma_{0}$, where $V$ is the electric potential, Wigner equations were derived from the Von-Neumann equation in [27]. In order to derive diffusion models, we consider two types of collision operators in the Wigner model.

First, we employ relaxation-time terms of BGK-type, $Q(w)=(g-w)/\tau_{c}$, where $w$ is the Wigner distribution, $\tau_{c}>0$ is the mean free path, and $g$ is the thermal equilibrium distribution derived from the quantum minimum entropy principle [5]; see Sections 2.1 and 2.2. We assume that the wave energy is much smaller than the typical kinetic energy (semiclassical hypothesis) and that the the scaled Planck constant is of the same order as the scaled Fermi energy (low scaled Fermi speed hypothesis). Performing a diffusive limit and a Chapman-Enskog expansion around the equilibrium distribution formally yields the first quantum spin diffusion model (QSDE1) for the particle density $n_{0}$ and the spin vector $\vec{n}=(n_{1},n_{2},n_{3})$, which are the zeroth-order moments of the Wigner distribution $w$:

[TABLE]

where $A_{0}$ and $B_{0}$ are some functions, and $F_{j}$ depends on $n_{0}$, $\vec{n}$, $\nabla\vec{n}$, and $\nabla V$. We refer to Section 2.3 for details. Note that the equations for the particle density and spin vector decouple; only the spin vector equation depends nonlinearly on $n_{0}$ and $\vec{n}$. The functions $A_{0}$ and $B_{0}$ are well defined only if $0\leq|\vec{n}|/n_{0}<1$. Hence, the main difficulty in the analysis of this model is the proof of lower and upper bounds for $|\vec{n}|/n_{0}$.

Second, we employ a mass- and spin-conserving matrix collision operator suggested in [20] for a semiconductor subject to a magnetic field. Performing a diffusive limit and a Chapman-Enskog expansion similarly as for the first model, we derive the second quantum spin diffusion model (QSDE2) in which the equations for the particle density $n_{0}$ and the spin vector $\vec{n}$ are fully coupled:

[TABLE]

where $A_{j}$, $B_{j}$, and $C_{j}$ are some functions depending also on the (given) pseudo-spin polarization and the direction of the local pseudo-magnetization, and $G_{j}$ depends on $\vec{n}$ and $J_{j}$. We refer to Section 2.4 for details. Because of the cross-diffusion structure, the analysis of this model is not immediate, and we solve this model only numerically (see Section 4).

Thanks to the decoupled structure of the model QSDE1, we are able to perform an analytical study. More precisely, we show in Section 3 the global existence and uniqueness of weak solutions, some entropy-dissipation properties, and the exponential long-time decay of the spin vector. As mentioned above, the main challenge is the proof of $|\vec{n}|/n_{0}<1$. By the maximum principle, it is not difficult to prove that $n_{0}$ is strictly positive. However, an application of the maximum principle to the equation for the spin vector is less obvious. Our idea is to show that $u=1-|\vec{n}|^{2}/n_{0}^{2}$ satisfies the equation

[TABLE]

where $G[\vec{n}/n_{0}]$ is some nonnegative function. This simple structure comes from the fact that certain antisymmetric terms in $A_{0}$ and $B_{0}$ cancel in this situation. By Stampacchia’s truncation method, we conclude that there exists a positive lower bound for $u$ which proves that $|\vec{n}|/n_{0}<1$.

Finally, we present in Section 4 some numerical results for the models QSDE1 and QSDE2, applied to a simple ballistic diode in one space dimension. The equations are discretized by a Crank-Nicolson finite-difference method. We illustrate the behavior of the particle density $n_{0}$ and the spin components $n_{j}$ for various instants of time and the exponential convergence of the particle density to the steady state.

The paper is organized as follows. Section 2 is concerned with the derivation of the models QSDE1 and QSDE2. The model QSDE1 is analyzed in Section 3. Finally, numerical experiments are presented in Section 4.

2. Modeling

2.1. A kinetic model for graphene

We describe the kinetic model for the quantum transport in graphene associated to the Hamiltonian $H+V\sigma_{0}$, where $H$ is given by (2), and $V$ is the electric potential. Let $w(x,p,t)$ denote the system Wigner distribution, depending on the position $x\in\mathbb{R}^{2}$, momentum $p\in\mathbb{R}^{2}$, and time $t\geq 0$. The Wigner function takes values in the space of complex Hermitian $2\times 2$ matrices, which is an Hilbert space with respect to the scalar product $(A,B)=\frac{1}{2}\mbox{tr}(AB)$, where $\mbox{tr}(A)$ denotes the trace of the matrix $A$. The set of Pauli matrices

[TABLE]

is a complete orthonormal system on that space. Therefore, we can develop the Wigner function $w$ in terms of the Pauli matrices, $w=\sum_{j=0}^{3}w_{j}\sigma_{j}$, where $w_{j}(x,p,t)$ are real-valued scalar functions. We set $\vec{w}=(w_{1},w_{2},w_{3})$, $\vec{p}=(p_{1},p_{2},0)$, $p=(p_{1},p_{2})$, $\vec{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$, and we abbreviate $\partial_{t}=\partial/\partial t$ and $\vec{\nabla}=(\partial/\partial x_{1},\partial/\partial x_{2},0)$. With this notation, we can write $w=w_{0}\sigma_{0}+\vec{w}\cdot\vec{\sigma}$. By applying the Wigner transform to the Von-Neumann equation, associated to the Hamilonian $H+V$, the following Wigner equations for the quantum transport in graphene have been derived in [27]:

[TABLE]

where $\hbar$ is the reduced Planck constant. The parameter $m$, which has the dimension of a mass, appears in the Hamiltonian $H$; see (2). The expressions $((\vec{p}/m)\cdot\vec{\nabla})w_{j}$ are originating from the quadratic term in the Hamiltonian $H$. Compared to Formula (12) in [27], we have allowed for BGK-type collision operators on the right-hand sides of (4) with the relaxation time $\tau_{c}$ and the thermal equilibrium distribution $g=g_{0}\sigma_{0}+\vec{g}\cdot\vec{\sigma}$ which is defined in Section 2.2. The pseudo-differential operator $\theta_{\hbar}[V]w$ is given by

[TABLE]

with its symbol

[TABLE]

In order to derive macroscopic diffusive models, we perform a diffusion scaling. We introduce a typical spatial scale $\widehat{x}$, time scale $\widehat{t}$, momentum scale $\widehat{p}$, and potential scale $\widehat{V}$:

[TABLE]

where the scales are related to

[TABLE]

Here, $T$ is the (constant) system temperature and $k_{B}$ the Boltzmann constant. The third relation means that the typical value of the momentum equals to the thermal momentum. Let $L$ denote the average distance which a particle travels with the Fermi velocity $v_{F}$ between two consecutive collisions, i.e. $L=\tau_{c}v_{F}$. Then the first relation can be written as

[TABLE]

Thus, the ratio of the typical length scale and the “Fermi mean free path” is assumed to be of the same order as the quotient of the electric/wave energies and the kinetic/Fermi energies. The second relation

[TABLE]

means that the ratio of the typical time scale and the relaxation time is of the same order as the quotient of the square of the wave energy and the kinetic/Fermi energies.

We introduce the semiclassical parameter $\varepsilon$, the diffusion parameter $\tau$, and the scaled Fermi speed $c$, given by

[TABLE]

We suppose the semiclassical hypothesis $\varepsilon\ll 1$ and the so-called low scaled Fermi speed hypothesis,

[TABLE]

With the above scaling, equations (4) become

[TABLE]

The “drift” terms $(\varepsilon/2)\vec{\nabla}\cdot\vec{w}$ and $(\varepsilon/2)\vec{\nabla}w_{0}$ are of order $O(\varepsilon)$, whereas the “precession” term $\vec{w}\wedge\vec{p}$ is of order one. This means that we have chosen a time scale which is of the same order as the magnitude of the precession period of the spin around the current, which is smaller than the typical time scale of the drift process.

2.2. Thermal equilibrium distribution

We define now the thermal equilibrium distribution $g=g_{0}\sigma+\vec{g}\cdot\vec{\sigma}$ using the minimum entropy principle. We introduce the (unscaled) quantum entropy by

[TABLE]

where $\mbox{Log}(w)=\mbox{Op}_{\varepsilon}^{-1}\log\mbox{Op}_{\varepsilon}(w)$ is the so-called quantum logarithm introduced by Degond and Ringhofer [5], $\mbox{Op}_{\varepsilon}$ is the Weyl quantization, defined for any symbol $\gamma(x,p)$ and any test function $\psi$ by [7, Chapter 2]

[TABLE]

and

[TABLE]

is the symbol of the Hamiltonian $H$, i.e. $H=\mbox{Op}_{\hbar}(h)$.

According to the theory of Degond and Ringhofer [5], we define the Wigner distribution at local thermal equilibrium related to the given functions $n_{0}$ and $\vec{n}$ as the formal solution $g=g[n_{0},\vec{n}]$ (if it exists) to the problem

[TABLE]

where the minimum is taken over all Wigner functions with complex Hermitian values and $w$ is decomposed according to $w=w_{0}\sigma_{0}+\vec{w}\cdot\vec{\sigma}$. This problem can be solved formally by means of Lagrange multipliers; see [26, Section 3.2]. For scalar-valued Wigner functions, such problems are studied analytically in [14]. Formally, the (scaled) solution is given by

[TABLE]

where $A=A(x,t)$ and $\vec{B}=\vec{B}(x,t)=(B_{1},B_{2},B_{3})(x,t)$ are the Lagrange multipliers determined by

[TABLE]

and $\mbox{Exp}(w)=\mbox{Op}_{\varepsilon}^{-1}\exp\mbox{Op}_{\varepsilon}(w)$ is the quantum exponential [5].

We wish to find an approximate but explicit expression for $g$. To this end, we expand the quantum exponential in terms of powers of $\varepsilon$, using the semiclassical and the low scaled Fermi speed hypotheses. The expansion follows the lines of Section 3.4 in [26]. We obtain from (7):

[TABLE]

Note that $a$ and $b$ are of order one, in view of (5). Employing formulas (29), (37), and (38) of [26], we deduce that $g=g^{(0)}+\varepsilon g^{(1)}+O(\varepsilon^{2})$, where

[TABLE]

where $I$ denotes the unit matrix in $\mathbb{R}^{3\times 3}$. Finally, it remains to express the Lagrange multipliers $A$ and $\vec{B}$ in terms of $n_{0}$ and $\vec{n}$ by means of the contraints (8). We find after tedious but straightforward computations that

[TABLE]

where $I$ denotes the unit matrix. Equations (9) and (10) provide an explicit approximation of the thermal equilibrium distribution. The functions $g_{0}$ and $\vec{g}$ are defined by $g^{(0)}+\varepsilon g^{(1)}=g_{0}\sigma_{0}+\vec{g}\cdot\vec{\sigma}$.

2.3. Derivation of the first model

We derive our first spinorial drift-diffusion model. We assume that both the semiclassical parameter $\varepsilon$ and the diffusion parameter $\tau$ in (6) are small and of the same order. We will perform the limit $\tau\to 0$ and $\varepsilon\to 0$, setting

[TABLE]

From (6) follows that the lowest-order approximations of $w_{0}$ and $\vec{w}$ are $g_{0}$ and $\vec{g}$, respectively. In order to compute the first-order approximation, we employ a Chapman-Enskog expansion of the Wigner function $w=w_{0}+\vec{w}\cdot\vec{\sigma}$ around the equilibrium distribution $g$. Inserting the expansions $w_{0}=g_{0}\sigma_{0}+\tau f_{0}$, $\vec{w}=\vec{g}+\tau\vec{f}$ into (6) and performing the formal limit $\tau\to 0$, we infer that

[TABLE]

Here, we have used the expansion $\theta_{\varepsilon}[V]=-\vec{\nabla}V\cdot\vec{\nabla}_{p}+O(\varepsilon)$.

The moment equations of (6) read as

[TABLE]

since $\int_{\mathbb{R}^{2}}\theta_{\varepsilon}[V]w_{j}dp=0$ for $j=0,1,2,3$. We need to compute the first-order moments

[TABLE]

in order to close the moment equations (11). For this, we insert in these integrals the expansions for $w_{0}$ and $\vec{w}$ as well as the expansions (9)-(10). After long but straightforward computations and rescaling $x\to x/(2\gamma)$ and $V\to V/(2\gamma)$ in order to get rid of the factor $1/(2\gamma)$, we arrive at the expressions, up to terms of order $O(\tau^{2})$,

[TABLE]

where we have set $\partial_{1}=\partial/\partial x_{1}$, $\partial_{2}=\partial/\partial x_{2}$, $\partial_{3}=0$,

[TABLE]

using Einstein’s summation convention, and $(\eta_{jk\ell})$ is the only antisymmetric 3-tensor which is invariant under cyclic index permutations such that $\eta_{123}=1$. In other words, $\eta_{jk\ell}a_{j}b_{k}=(\vec{a}\wedge\vec{b})_{\ell}$ for $\vec{a}$, $\vec{b}\in\mathbb{R}^{3}$. Inserting these expressions into (11), we find the model QSDE1:

[TABLE]

where $j$, $s=1,2,3$. The functions

[TABLE]

satisfy $0<b_{k}[\vec{v}]<\lambda$ for all $0<|\vec{v}|<1$, $\lim_{|\vec{v}|\to 0}b_{k}[\vec{v}]=0$, and the function $|\vec{v}|\mapsto\Phi(|\vec{v}|)=b_{k}[\vec{v}]/(\lambda v_{k})$ is increasing. This allows us to set $b_{k}[0]=0$ such that $b_{k}$ is defined for all $0\leq|\vec{v}|<1$.

The above equations are complemented by the Poisson equation

[TABLE]

for the electric potential, where $\lambda_{D}>0$ is the scaled Debye length.

2.4. Derivation of the second model

In the model (12)-(14), the particle density $n_{0}$ evolves independently from the spin vector $\vec{n}$. We will modify this model in order to derive a fully coupled system by adding a “pseudo-magnetic” field which is supposed to be able to interact with the charge carrier pseudo-spin.

Possanner and Negulescu [20] consider a semiconductor subject to a magnetic field which interacts with the electron spin and they build a purely semiclassical diffusive model for the particle density $n_{0}$ and the spin vector $\vec{n}$ by a Chapman-Enskog expansion around the equilibrium distribution. Instead of the relaxation-time model used in Section 2.3, we employ here the mass and spin conserving collision operator (49) of [20],

[TABLE]

where $g$ is the equilibrium distribution, defined in Section 2.2, and $P=\sigma_{0}+\zeta\vec{\omega}\cdot\vec{\sigma}$ is the polarization matrix with the pseudo-spin polarization $\zeta(x,t)$ of the scattering rate and $\vec{\omega}(x,t)$ is the direction of the pseudo-magnetization (see [20, Section 4.1]). The quantity $\zeta(x,t)\in(0,1)$ satisfies

[TABLE]

where $s_{\uparrow\downarrow}$ are the scattering rates of electrons in the upper and lower band, respectively, and it holds $|\vec{\omega}(x,t)|=1$. We introduce the collision operators $Q_{0}$ and $\vec{Q}$ by

[TABLE]

We start from the scaled Wigner equations

[TABLE]

Compared with the Wigner system (6) in Section 2, the second equation in (16) contains the (heuristic) term $\tau\vec{\omega}\wedge\vec{w}$, which describes the “precession” of $\vec{w}$ around the local pseudo-magnetization. We assume again that $\lambda:=\varepsilon\gamma/\tau$ is of order one (as $\tau\to 0$) and we perform a Chapman-Enskog expansion of the Wigner distribution $w=w_{0}\sigma_{0}+\vec{w}\cdot\vec{\sigma}$. The result reads as follows:

[TABLE]

where

[TABLE]

To compute $P^{-1/2}T[g]P^{-1/2}$, we employ the following lemma whose proof is an elementary computation.

Lemma 1**.**

For all Hermitian matrices $a=a_{0}\sigma_{0}+\vec{a}\cdot\vec{\sigma}$, it holds

[TABLE]

Since the collision operator $Q(w)$ is conserving mass and spin, the moment equations of (16) become

[TABLE]

The first-order moments $\int_{\mathbb{R}^{2}}p_{k}w_{0}dp$ and $\int_{\mathbb{R}^{2}}p_{k}w_{s}dp$ can be calculated up to order $O(\varepsilon^{2})=O(\tau^{2})$ by using (17) and Lemma 1, leading to the spin diffusion model QSDE2:

[TABLE]

where

[TABLE]

In contrast to the model QSDE1 (12)-(14), this system is fully coupled. It is possible to show that the system is uniformly parabolic if $\|\zeta\|_{L^{\infty}(0,T;L^{\infty}(\Omega))}<1$ but the presence of the cross-diffusion terms makes it hard to prove any $L^{\infty}$ bounds, in particular the bound $|\vec{n}|/n_{0}<1$.

3. Analysis for the first model

In this section, we consider the model QSDE1 (12)-(15) in a bounded domain $\Omega\subset\mathbb{R}^{2}$ with $\partial\Omega\in C^{1,1}$. For convenience, we recall the equations:

[TABLE]

and the functions

[TABLE]

We impose the following boundary and initial conditions:

[TABLE]

Finally, we abbreviate $\Omega_{T}=\Omega\times(0,T)$.

3.1. Existence of solutions

We impose the following conditions on the data:

[TABLE]

for some $p>2$. Under these assumptions, we are able to prove the existence of strong solutions $(n_{0},V)$ to the drift-diffusion model (20) and (23).

Theorem 2**.**

Let $T>0$ and assume (26)-(28). Then there exists a unique solution $(n_{0},V)$ to (20) and (23) subject to the initial and boundary conditions in (24)-(25) satisfying

[TABLE]

where $\mu=\lambda_{D}^{-2}$ and

[TABLE]

Proof.

The existence and uniqueness of a weak solution $(n_{0},V)$ to (20), (23), and (24)-(25) satisfying

[TABLE]

is shown in [9], also see Section 3.9 in [13]. It remains to prove the regularity assertions.

First, we show that $n_{0}$ is bounded from above and below. Employing $(n_{0}-M)^{+}=\max\{0,n_{0}-M\}\in L^{2}(0,T;H^{1}_{0}(\Omega))$ as a test function in the weak formulation of (20) and using (23), we find that

[TABLE]

since $n_{0}-C(x)\geq 0$ on $\{n_{0}>M\}$, by the definition of $M$. This implies that $(n_{0}-M)^{+}=0$ and hence, $n_{0}\leq M$ in $\Omega$, $t>0$.

Next, we employ the test function $(n_{0}-me^{-\mu t})^{-}=\min\{0,n_{0}-me^{-\mu t}\}$ in the weak formulation of (20):

[TABLE]

since we integrate over $\{n_{0}<me^{-\mu t}\}$. By the definition of $\mu$, the last integral vanishes. Then, the Gronwall lemma implies that $(n_{0}-me^{-\mu t})^{-}=0$ and hence, $n_{0}\geq me^{-\mu t}$.

The above bounds show that the right-hand side of the Poisson equation is an element of $L^{\infty}(0,T;L^{\infty}(\Omega))$. Then, by elliptic regularity, $V\in L^{\infty}(0,T;W^{2,p}(\Omega))$, where $p>2$ is given in (28). Since $W^{2,p}(\Omega)\hookrightarrow W^{1,\infty}(\Omega)$ (we recall that $\Omega\subset\mathbb{R}^{2}$), it follows that $\nabla V\in L^{\infty}(0,T;L^{\infty}(\Omega))$. Consider

[TABLE]

The right-hand side of this equation satisfies $\partial_{t}n_{0}\in L^{2}(0,T;(H^{1}(\Omega))^{\prime})$. Hence, $\partial_{t}V\in L^{2}(0,T;H^{1}(\Omega))$.

Finally, we prove the higher regularity for $n_{0}$. For this, we consider the equation satisfied by $\rho=n_{0}-n_{D}$:

[TABLE]

where $f=-\partial_{t}n_{D}+\textnormal{div}(\nabla n_{D}+n_{D}\nabla V)$. We use the following result: If $f$, $\partial_{t}f\in L^{2}(\Omega_{T})$ and $\rho(\cdot,0)\in H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$ then $n_{0}\in C^{0}([0,T];H^{2}(\Omega))$, $\partial_{t}n_{0}\in L^{2}(0,T;H^{1}(\Omega))$, $\partial_{t}^{2}n_{0}\in L^{2}(0,T;(H^{1}(\Omega))^{\prime})$ (see, e.g., [28, Theorem 1.3.1]).

The boundedness of $\nabla V$ implies that $f\in L^{2}(\Omega_{T})$. This gives the regularity $\partial_{t}n_{0}\in L^{2}(\Omega_{T})$. As a consequence, $V\in H^{1}(0,T;H^{2}(\Omega))$. Our assumptions on the data show that $-\partial_{t}^{2}n_{D}+\Delta\partial_{t}n_{D}\in L^{2}(\Omega_{T})$, and it remains to prove that $\partial_{t}\textnormal{div}(n_{D}\nabla V)\in L^{2}(\Omega_{T})$. Now,

[TABLE]

The first term on the right-hand side lies in $L^{2}(\Omega_{T})$ since $\nabla V\in L^{\infty}(0,T;L^{\infty}(\Omega))$. Furthermore, $\nabla n_{D}\in L^{\infty}(0,T;H^{1}(\Omega))$ and $\nabla\partial_{t}V\in L^{2}(0,T;H^{1}(\Omega))$ from we conclude that the second term is in $L^{2}(\Omega_{T})$. In a similar way, this property can be verified for the third and fourth terms. This shows the claim and the regularity statements for $n_{0}$. ∎

Next, given $(n_{0},V)$ as the solution to (20) and (23), we prove the existence of a solution $\vec{n}$ to (21) with the corresponding boundary and initial conditions in (24)-(25), satisfying the bound $|\vec{n}|/n_{0}<1$.

Theorem 3**.**

Let $(n_{0},V)$ be the solution to (20), (23), and (24)-(25), according to Theorem 2, and let $\vec{n}_{I}\in H_{0}^{1}(\Omega)^{3}$ satisfy

[TABLE]

Then there exists a weak solution $\vec{n}\in L^{2}(0,T;H^{1}(\Omega))^{3}\cap H^{1}(0,T;H^{-1}(\Omega))^{3}$ to (21)-(22) and (24)-(25) satisfying

[TABLE]

Furthermore, there exists at most one weak solution satisfying (29) and $\vec{n}\in L^{\infty}(0,T;$ $W^{1,4}(\Omega))^{3}$.

Proof.

We prove first the existence of solutions to a truncated problem by applying the Leray-Schauder fixed-point theorem. Let $0<\chi<1$ be a fixed parameter and let $\phi_{\chi}\in C^{0}(\mathbb{R})$ be a nonincreasing function satisfying $\phi_{\chi}(y)=1$ for $y\leq 1-\chi$ and $\phi_{\chi}(y)=0$ for $y\geq 1$. Then define

[TABLE]

Step 1: Application of the fixed-point theorem. In order to define the fixed-point operator, let $\vec{\rho}\in L^{2}(\Omega_{T})^{3}$ and $\sigma\in[0,1]$. We wish to solve the linear problem

[TABLE]

where

[TABLE]

for $\vec{z}\in H_{0}^{1}(\Omega)^{3}$. The bilinear form $a:H_{0}^{1}(\Omega)^{3}\times H_{0}^{1}(\Omega)^{3}\to\mathbb{R}$ is continuous since $|b_{k}^{\chi}[\vec{\rho}/n_{0}]|\leq\lambda$ and $|\nabla V|\in L^{\infty}(0,T;L^{\infty}(\Omega))$. Furthermore, using the antisymmetry of $\eta_{jk\ell}$,

[TABLE]

where $\|\nabla\vec{n}\|^{2}=\sum_{j,k=1}^{3}(\partial_{j}n_{k})^{2}$. All the terms on the right-hand side can be written as a product of $n_{j}$, $\partial_{k}n_{\ell}$, and possibly an $L^{\infty}$ function. Note that the only term, which does not have this structure, $b_{k}^{\chi}\eta_{jk\ell}\partial_{s}n_{\ell}\partial_{s}n_{j}$, vanishes because of the antisymmetry of $\eta_{jk\ell}$. Therefore, the Hölder and Cauchy-Schwarz inequalities yield

[TABLE]

for some constant $c>0$ which depends on the $L^{\infty}$ norm of $\nabla V$. Hence, there exists a unique weak solution $\vec{n}\in L^{2}(0,T;H^{1}_{0}(\Omega))^{3}\cap H^{1}(0,T;H^{-1}(\Omega))^{3}$ to the linear problem (30) [25, Corollary 23.26]. Moreover, there exists a constant $c>0$ independent of $\rho$ and $\sigma$ such that

[TABLE]

This defines the fixed-point operator $F:L^{2}(\Omega_{T})^{3}\times[0,1]\to L^{2}(\Omega_{T})^{3}$, $F(\vec{\rho},\sigma)=\vec{n}$. We note that $F(\vec{\rho},0)=0$.

Next, we show that $F$ is continuous. Let $(\vec{\rho}^{(k)})\subset L^{2}(\Omega_{T})^{3}$ and $(\sigma^{(k)})\subset\mathbb{R}$ such that $\vec{\rho}^{(k)}\to\vec{\rho}$ in $L^{2}(\Omega_{T})^{3}$ and $\sigma^{(k)}\to\sigma$ as $k\to\infty$. Since $b_{k}^{\chi}$ is bounded, it follows that $b_{j}^{\chi}[\rho^{(k)}/n_{0}]\to b_{j}^{\chi}[\rho/n_{0}]$ in $L^{r}(\Omega_{T})$ for all $r<\infty$. Let $n^{(k)}=F(\rho^{(k)},\sigma^{(k)})$. The uniform estimate (31) shows that, up to a subsequence,

[TABLE]

since the embedding $L^{2}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;H^{-1}(\Omega))\hookrightarrow L^{2}(\Omega_{T})$ is compact, by the Aubin lemma. These convergence results are sufficient to perform the limit $k\to\infty$ in the weak formulation of (30) with $\vec{n}^{(k)}$ instead of $\vec{n}$ and $\sigma^{(k)}$ instead of $\sigma$. The limit equation shows that $\vec{n}=F(\vec{\rho},\sigma)$. As the solution to the linear problem is unique, the convergence $\vec{n}^{(k)}\to\vec{n}$ in $L^{2}(\Omega_{T})^{3}$ holds for the whole sequence, and hence, $S$ is continuous. Furthermore, by the Aubin lemma, $F$ is compact. Finally, let $\vec{n}$ be a fixed point of $F(\cdot,\sigma)$. By the boundedness of $b_{j}^{\chi}$ and estimate (31), we find uniform estimates for $\vec{n}$ in $L^{2}(\Omega_{T})^{3}$. Therefore, we can apply the fixed-point theorem of Leray-Schauder which yields a weak solution to the truncated problem (30) with $\vec{\rho}$ replaced by $\vec{n}$.

Step 2: $L^{\infty}$ bounds for $\vec{n}$. We prove that the solution to (30) is bounded in $\Omega_{T}$. To this end, we define the function $\psi=\sqrt{1+|\vec{n}|^{2}}$. Then, in the sense of distributions,

[TABLE]

where

[TABLE]

Inserting the identities

[TABLE]

in the above equation for $\psi$, we deduce that

[TABLE]

Since $|\textrm{curl}\,\vec{v}|^{2}\leq 2\|\nabla\vec{v}\|^{2}$ for all $\vec{v}\in H_{0}^{1}(\Omega)^{3}$, Young’s inequality gives

[TABLE]

Elementary computations show that

[TABLE]

Hence, we can estimate (33) by

[TABLE]

By the maximum principle, $\psi\leq c(T)$ in $\Omega_{T}$, where $c(T)>0$ depends on the end time $T>0$. Taking into account that $\psi\geq 1$ by definition, we conclude that $\psi\in L^{\infty}(0,T;L^{\infty}(\Omega))$ and hence, $|\vec{n}|\in L^{\infty}(0,T;L^{\infty}(\Omega))$.

Step 3: Proof of $|\vec{n}|/n_{0}<1$. We show that there exists $0<\kappa<1$ such that $|\vec{n}|/n_{0}\leq\kappa<1$. Then, choosing $\chi>0$ sufficiently small, we can remove the truncation obtaining a solution to the original problem.

Let $u=1-|\vec{n}|^{2}/n_{0}^{2}$. Note that this function is well defined since $n_{0}$ is strictly positive, see Theorem 2. A tedious computation shows that $u$ solves

[TABLE]

in the sense of distributions, where $G$ is defined in (32). We prove a lower bound for $u$ by testing the weak formulation of this equation by $U:=(u-k)^{-}=\min\{0,u-k\}\in L^{2}(0,T;H_{0}^{1}(\Omega))$, where $k=\min\{\inf_{\partial\Omega\times(0,T)}u,\inf_{\Omega\times\{0\}}u\}>0$.

[TABLE]

Since $G(\vec{n}/n_{0}]\geq 0$, the last integral is nonpositive. The first integral is estimated by employing the lower bound for $n_{0}$, obtained in Theorem 2, and applying Young’s and Hölder’s inequalities:

[TABLE]

The Gagliardo-Nirenberg, Hölder, and Poincaré inequalities imply that

[TABLE]

where $\varepsilon>0$. By Theorem 2, $\nabla n_{0}$ is an element of $L^{\infty}(0,T;W^{1,4}(\Omega))$, in view of the embedding $H^{2}(\Omega)\hookrightarrow W^{1,4}(\Omega)$. Collecting the above estimates, we infer that

[TABLE]

Inserting this inequality into (35) and choosing $\varepsilon<2$, it follows that

[TABLE]

Then Gronwall’s lemma and $U(\cdot,0)=0$ gives $U=0$ in $\Omega_{T}$ and hence, $u\geq k>0$ in $\Omega_{T}$.

Step 4: Uniqueness of solutions. Let $\vec{u}$ and $\vec{v}$ be two solutions to (21) and (24)-(25) satisfying (29) and $\vec{u}\in L^{\infty}(0,T;W^{1,4}(\Omega))^{3}$. Set $\vec{w}=\vec{u}-\vec{v}$. Taking the difference of the equations satisfied by $\vec{u}$ and $\vec{v}$, respectively, and employing $\vec{w}$ as a test function, we find that

[TABLE]

Thanks to the $L^{\infty}$ bounds on $\nabla V$, $\vec{u}$, and $\vec{v}$, we can estimate as follows:

[TABLE]

where $c>0$ is some generic constant. Since $\vec{w}(0)=0$, the $W^{1,4}$ regularity for $\vec{u}$ and Gronwall’s lemma imply the assertion. This finishes the proof. ∎

3.2. Entropy dissipation

Let $(n_{0},\vec{n},V)$ be a solution to (20)-(23), (24)-(25) according to Theorems 2 and 3. We assume that the boundary data is in thermal equilibrium, i.e.

[TABLE]

where $V_{D}=V_{D}(x)$ is time-independent. In this subsection, we will show that the macroscopic entropy

[TABLE]

is nonincreasing in time. Note that $n_{0}<|\vec{n}|$ by Theorem 3 such that $\log(n_{0}-|\vec{n}|)$ is well defined.

The functional $S(t)$ can be derived as follows. Inserting the thermal equilibrium distribution $g[n_{0},\vec{n}]$ in the quantum entropy $A[w]$, defined in Section 2.2, and taking into account the electric energy contribution, it follows that the total macroscopic free energy reads as

[TABLE]

Then the expansion of $g[n_{0},\vec{n}]$ (see (7), (9), and (10)) yields the above formula for $\widetilde{S}(t)=S(t)+O(\varepsilon^{2})$.

Proposition 4**.**

The entropy dissipation $-dS/dt$ can be written as

[TABLE]

where $G$ is defined in (32).

Note that in the drift-diffusion model without spin contribution, i.e. $|\vec{n}|=0$, we recover the standard entropy dissipation term $\int_{\Omega}n_{0}|\nabla(\log n_{0}+V)|^{2}dx$.

Proof.

Taking the time derivative of $S$, we find after some computations that

[TABLE]

To be precise, the second term in the integral has to be understood in the sense of $L^{2}(0,T;H^{-1}(\Omega))$. Since $V_{D}$ does not depend on time, we have $\partial_{t}V=0$ on $\partial\Omega$, $t>0$. Hence,

[TABLE]

and thus, the last two terms in $dS/dt$ cancel.

In order to compute the second term in $dS/dt$, we observe that

[TABLE]

where $G$ is defined in (32). Then, inserting this expression and (12) into (37) and integrating by parts, we infer that

[TABLE]

Since

[TABLE]

we can write

[TABLE]

Straightforward computations show that

[TABLE]

which allows us to reformulate the second integral. Together with some manipulations in the first integral, we obtain

[TABLE]

We observe that the expression

[TABLE]

is integrable because $\inf_{\Omega_{T}}n_{0}>0$, $\sup_{\Omega_{T}}|\vec{n}|/n_{0}<1$, the map

[TABLE]

is bounded for all $\varepsilon>0$ and $|\vec{n}|^{2}G[\vec{n}/|\vec{n}|]\in L^{1}(\Omega_{T})$. Since

[TABLE]

this finishes the proof. ∎

3.3. Long-time decay of the solutions

Let $(n_{0},\vec{n},V)$ be a solution to (20)-(23), (24)-(25) according to Theorems 2 and 3. We will show that, under suitable assumptions on the electric potential, the spin vector converges to zero as $t\to\infty$.

Theorem 5**.**

Let $2<p<\infty$. Then there exists a constant $\varepsilon_{p}>0$ such that if the condition $\|\nabla V\|_{L^{\infty}(0,T;L^{\infty}(\Omega))}\leq\varepsilon_{p}$ holds then

[TABLE]

for some constant $\kappa_{p}>0$ which depends on $p$, $\Omega$, and the $L^{\infty}$ norm of $\nabla V$. Furthermore, there exists $\varepsilon_{2}>0$ such that if $\|\Delta V\|_{L^{\infty}(0,T;L^{\infty}(\Omega))}\leq\varepsilon_{2}$ then

[TABLE]

for some constant $\kappa_{2}>0$ which depends on $\Omega$ and the $L^{\infty}$ norm of $\Delta V$.

Note that we may set $T=\infty$ yielding the desired convergence result. For the proof of the second part of the theorem, we need the following lemma.

Lemma 6**.**

There exists a constant $c_{G}>0$, depending only on $\Omega$, such that for all $\vec{u}\in H^{1}_{0}(\Omega)^{3}$,

[TABLE]

where $G$ is defined in (32).

Proof.

Let $\mu>0$ and consider the following bilinear form on $H_{0}^{1}(\Omega)^{3}$:

[TABLE]

where $\nabla\vec{u}:\nabla\vec{v}=\sum_{j,k=1}^{3}\partial_{j}u_{k}\partial_{j}v_{k}$. The bilinear form $B$ is symmetric, continuous, and coercive on $H_{0}^{1}(\Omega)^{3}$, since, using (34) and the Poincaré inequality,

[TABLE]

Hence, by the Lax-Milgram lemma, for all $\vec{f}\in L^{2}(\Omega)^{3}$, there exists a unique solution $\vec{u}\in H_{0}^{1}(\Omega)^{3}$ to

[TABLE]

This defines the linear operator $L:L^{2}(\Omega)^{3}\to L^{2}(\Omega)^{3}$, $L(\vec{f})=\vec{u}$. Since the range of $L$, $R(L)$, is a subset of $H_{0}^{1}(\Omega)^{3}$, which embeddes compactly into $L^{2}(\Omega)^{3}$, $L$ is compact. Furthermore, $L$ is symmetric (since $B$ is symmetric) and positive in the sense of $\int_{\Omega}L(\vec{f})\cdot\vec{f}dx>0$ for all $\vec{f}\neq 0$ (since $B$ is coercive). By the Hilbert-Schmidt theorem for symmetric compact operators (see, e.g., [21, Theorem VI.16]), there exists a complete orthonormal system $(\vec{u}^{(k)})$ of $L^{2}(\Omega)$ of eigenfunctions of $L$,

[TABLE]

Note that $(\lambda_{k})$ depends on $\mu$ since $L$ and $B$ do so. In particular, $\vec{u}^{(k)}\in R(L)\subset H_{0}^{1}(\Omega)^{3}$ and

[TABLE]

We claim that $\lambda_{1}^{-1}>\mu$. Otherwise, if $\lambda_{1}^{-1}\leq\mu$, the definition of $B$ and (38) yield

[TABLE]

The terms containing $\mu$ cancel, which gives

[TABLE]

However, in view of (34), the integral is nonnegative and hence, it must vanish. Therefore,

[TABLE]

On the other hand, using $|\textrm{curl}\,\vec{u}^{(1)}|^{2}\leq 2\|\nabla\vec{u}^{(1)}\|^{2}$,

[TABLE]

from which we infer that $\textrm{curl}\,\vec{u}^{(1)}=2\vec{u}^{(1)}$ and, by (39), $\|\nabla\vec{u}^{(1)}\|^{2}+6|\vec{u}^{(1)}|^{2}=0$. This implies that $\vec{u}^{(1)}=0$ which is absurd. Hence, $\lambda_{1}^{-1}>\mu$.

Now let $\vec{u}\in H_{0}^{1}(\Omega)^{3}\cap H^{2}(\Omega)^{3}$. We can decompose $\vec{u}$ in the orthonormal set $(\vec{u}^{(k)})$, $\vec{u}=\sum_{k\in\mathbb{N}}c_{k}\vec{u}^{(k)}$ for some $c_{k}\in\mathbb{R}$. It follows from (38) and the orthogonality of $(\vec{u}^{(k)})$ on $L^{2}(\Omega)^{3}$ that

[TABLE]

By a density argument, this inequality also holds for all $\vec{u}\in H_{0}^{1}(\Omega)^{3}$. Therefore,

[TABLE]

and the lemma follows with $c_{G}=\lambda_{1}^{-1}-\mu>0$. ∎

Proof of Theorem 5..

Let $p\geq 2$. Using the test function $p|\vec{n}|^{p-2}n_{j}\in L^{2}(0,T;H_{0}^{1}(\Omega))$ in (21) and summing over $j=1,2,3$, we find after elementray computations that

[TABLE]

Now, we distinguish two cases. First, let $p>2$. Employing Young’s inequality with $\alpha>0$ and $G[\vec{n}]\geq 0$, we find that

[TABLE]

Then, choosing $\alpha=4(p-2)/p$ and employing the Poincaré inequality

[TABLE]

we infer that

[TABLE]

This proves the first part after setting $\varepsilon_{p}<\sqrt{8}(p-2)C_{P}/(p(p-1))$ and applying the Gronwall lemma.

For the second part, let $p=2$ in (40). By Lemma 6 and integration by parts in the term containing the potential, we obtain

[TABLE]

With the choice $\varepsilon_{2}<2c_{G}$ and the Gronwall lemma, the theorem follows. ∎

If the total space charge $n_{0}-C(x)=-\lambda_{D}^{2}\Delta V$ is positive, we are able to prove that $\vec{n}(\cdot,t)$ converges to zero in the $L^{\infty}$ norm.

Proposition 7**.**

Let $0<T\leq\infty$. The following $L^{\infty}$ estimate holds:

[TABLE]

Proof.

Let $2<p<\infty$ be arbitrary. From (40) we deduce that, by integrating by parts,

[TABLE]

Therefore

[TABLE]

Passing to the limit $p\to\infty$ in this inequality yields the claim. ∎

4. Numerical simulations

In this section, we present some numerical results for the models (12)-(15) and (18)-(19), (15) in one space dimension, $\Omega=(0,1)$. We choose the boundary conditions

[TABLE]

where $U(x)=V_{A}x$ and $V_{A}=80$ is the scaled applied potential, and the initial conditions

[TABLE]

where $V_{\rm eq}$ is the equilibrium potential, defined by

[TABLE]

We choose $\lambda_{D}^{2}=10^{-3}$. The doping profile corresponds to that of a ballistic diode:

[TABLE]

where $C_{\rm min}=0.025$ and $\bar{x}=0.2$. The pseudo-spin polarization and the direction of the local magnetization are defined by

[TABLE]

Table 1 shows the values of the units which allows for the computation of the physical values from the scaled ones.

Models QSDE1 (12)-(15) and QSDE2 (18)-(19), (15) with the corresponding initial and boundary conditions (24)-(25) are discretized with the Crank-Nicolson finite-difference scheme and the space step $\triangle x=10^{-2}$. The resulting nonlinear discrete ODE system is solved by using the Matlab routine ode23s.

Since the initial spin vector is assumed to vanish, the particle density $n_{0}$, computed from the model QSDE1, corresponds exactly to the particle density of the standard drift-diffusion model, and the spin vector vanishes for all time. This situation is different in the model QSDE2 since the equations are fully coupled. For the model QSDE2, Figure 1 shows the particle density $n_{0}$ and the components $n_{j}$ of the spin vector versus position at various times. The solution at $t=1$ corresponds to the steady state. We observe a charge built-up of $n_{0}$ in the low-doped region of the diode. The spin vector components vary only slightly in this region but their gradients are significant in the high-doped regions close to the contacts. Clearly, the components $n_{j}$ do not need to be positive and, in fact, they even do not have a sign.

The models QSDE1 and QSDE2 are well defined only if $|\vec{n}|/n_{0}<1$. We plot this ratio in Figure 2 at various times for the model QSDE2. In all the presented cases, the quotient stays below one. This indicates that $b_{k}[\vec{n}/n_{0}]$ is well defined also in this model.

We have shown in Theorem 5 that the spin vector of the model QSDE1 converges to zero if the electric potential satisfies certain conditions. In Figure 3, the relative difference $\|n_{0}(t)-n_{0}(\infty)\|_{2}/\|n_{0}(\infty)\|_{2}$ versus time is depicted (semilogarithmic plot), where $n_{0}(\infty)$ denotes the steady-state particle density of model QSDE1 or QSDE2, respectively. The norm $\|\cdot\|_{2}$ is the Euclidean norm. The stationary solution is approximated by $n_{0}(t^{*})$ with $t^{*}=1$. Whereas the decay of the solution to the model QSDE1 is numerically of exponential type (in agreement with the theoretical results), the decay for the model QSDE2 seems to be exponential only for small times.

In the final Figure 4, we present the current-voltage characteristics for the models QSDE1 and QSDE2, i.e. the relation between $J_{0}$ at $x=1$ and the applied bias $V_{A}$. The characteristics of model QSDE1 correspond to the current-voltage curve of the standard drift-diffusion model. We observe that the additional terms in the definition of $J_{0}$ lead to an increase of the particle current density.