Energy-transport systems for optical lattices: derivation, analysis, simulation

TL;DR

This paper derives and analyzes energy-transport equations for fermions in optical lattices, exploring their mathematical properties and numerical behavior, including existence of solutions and steady-state convergence.

Contribution

It provides a formal derivation of energy-transport equations from Boltzmann equations and proves existence of weak solutions with numerical validation.

Findings

Particle density converges to a constant steady state with high initial energy density.

Nonconstant steady states can occur depending on initial conditions.

The first-order expansion improves the physical realism of the model.

Abstract

Energy-transport equations for the transport of fermions in optical lattices are formally derived from a Boltzmann transport equation with a periodic lattice potential in the diffusive limit. The limit model possesses a formal gradient-flow structure like in the case of the energy-transport equations for semiconductors. At the zeroth-order high temperature limit, the energy-transport equations reduce to the whole-space logarithmic diffusion equation which has some unphysical properties. Therefore, the first-order expansion is derived and analyzed. The existence of weak solutions to the time-discretized system for the particle and energy densities with periodic boundary conditions is proved. The difficulties are the nonstandard degeneracy and the quadratic gradient term. The main tool of the proof is a result on the strong convergence of the gradients of the approximate solutions. Numerical simulations in one space dimension show that the particle density converges to a constant steady state if the initial energy density is sufficiently large, otherwise the particle density converges to a nonconstant steady state.

Equations

$${\partial }_{t}f+u\cdot {\nabla }_{x}f+{\nabla }_{x}V\cdot {\nabla }_{p}f=Q(f),$$

$$\epsilon (p)=-2{\epsilon }_0\sum _{i=1}^d\cos (2\pi p_i),p\in {\mathbb{T}}^d.$$

$$Q(f)=\frac{1}{\tau }({\mathcal{F}}_f-f),$$

$${\mathcal{F}}_{f}(x,p,t)=\big{(}\eta+\exp(-\lambda_{0}(x,t)-\lambda_{1}(x,t)\varepsilon(p))\big{)}^{-1},\quad x\in{\mathbb{R}}^{d},\ p\in{\mathbb{T}}^{d},\ t>0,$$

$${\partial }_{t}n+\mathrm{div}{J}_n=0,{\partial }_{t}E+\mathrm{div}{J}_E-J_n\cdot \nabla V=0$$

$$J_n=-\sum _{i=0}^1{D}_{0i}\nabla {\lambda }_i-{\lambda }_1{D}_{00}\nabla V,J_E=-\sum _{i=0}^1{D}_{1i}\nabla {\lambda }_i-{\lambda }_1{D}_{10}\nabla V,$$

$$\frac{dH}{dt}=-{\int }_{{\mathbb{R}}^d}\sum _{i,j=0}^1\nabla {\mu }_i^{\top }{L}_{ij}\nabla {\mu }_{j}dx\le 0;$$

$$\partial_{t}\begin{pmatrix}n\\ E\end{pmatrix}=\operatorname{div}\bigg{(}\begin{pmatrix}L_{00}&L_{01}\\ L_{10}&L_{11}\end{pmatrix}\nabla\begin{pmatrix}\mu_{0}\\ \mu_{1}\end{pmatrix}\bigg{)},$$

$$\partial_{t}n=\operatorname{div}\bigg{(}\frac{\tau_{0}\nabla n}{n(1-\eta n)}\bigg{)}\quad\mbox{in }{\mathbb{R}}^{d}.$$

$${\partial }_{t}n$$

$${\partial }_{t}W$$

$$\partial_{t}W_{\rm tot}=\operatorname{div}\bigg{(}\frac{\nabla W}{2n(1-\eta n)}+\frac{UW}{1-\eta n}\nabla n\bigg{)}.$$

$$J_n=-\nabla (n\theta )+n\nabla V,J_E=-{\kappa }_0\nabla \theta +\frac{5}{2}\theta J_n,E=\frac{3}{2}n\theta ,$$

$${\partial }_{t}n=\Delta ({\theta }^{\alpha }n),{\partial }_t(\theta n)=\Delta ({\theta }^{\alpha +1}n)+\frac{n}{\tau }(1-\theta )$$

$$\frac{1}{△t}(n^k-n^{k-1})$$

$$\frac{1}{△t}(W^k-W^{k-1})$$

$$\alpha {\partial }_t{f}_{\alpha }+u\cdot {\nabla }_x{f}_{\alpha }+\nabla V_{\alpha }\cdot {\nabla }_p{f}_{\alpha }=\frac{1}{\alpha }{Q}_{\alpha }(f_{\alpha }),(x,p)\in {\mathbb{R}}^d\times {\mathbb{T}}^d,t>0,$$

$$Q_{\alpha }(f)=\frac{1}{{\tau }_{\alpha }}({\mathcal{F}}_f-f),$$

$${\int }_{{\mathbb{T}}^d}({\mathcal{F}}_f-f)dp=0,{\int }_{{\mathbb{T}}^d}({\mathcal{F}}_f-f)\epsilon (p)dp=0,$$

$${\mathcal{F}}_f(x,p,t)=\frac{1}{\eta +\exp (-{\lambda }_0(x,t)-{\lambda }_1(x,t)\epsilon (p))},$$

$$n=n[{\mathcal{F}}_f]={\int }_{{\mathbb{T}}^d}{\mathcal{F}}_{f}dp,E=E[{\mathcal{F}}_f]={\int }_{{\mathbb{T}}^d}{\mathcal{F}}_f\epsilon (p)dp$$

$$D_{ij}^{k\ell }=\tau {\int }_{{\mathbb{T}}^d}{u}_k{u}_{\ell }{\mathcal{F}}_f(1-\eta {\mathcal{F}}_f)\epsilon (p{)}^{i+j}dp,i,j=0,1,k,\ell =1,\dots ,d.$$

$${\partial }_{t}n[f_{\alpha }]+\frac{1}{\alpha }{\mathrm{div}}_x{\int }_{{\mathbb{T}}^d}u{f}_{\alpha }dp=0,$$

$${\partial }_{t}E[f_{\alpha }]+\frac{1}{\alpha }{\mathrm{div}}_x{\int }_{{\mathbb{T}}^d}\epsilon u{f}_{\alpha }dp-\frac{1}{\alpha }{\nabla }_x{V}_{\alpha }\cdot {\int }_{{\mathbb{T}}^d}u{f}_{\alpha }dp=0.$$

$${\partial }_{t}n[{\mathcal{F}}_{f_{\alpha }}]+\alpha {\partial }_{t}n[g_{\alpha }]+{\mathrm{div}}_x{\int }_{{\mathbb{T}}^d}u{g}_{\alpha }dp=0,$$

$${\partial }_{t}E[{\mathcal{F}}_{f_{\alpha }}]+\alpha {\partial }_{t}E[g_{\alpha }]+{\mathrm{div}}_x{\int }_{{\mathbb{T}}^d}u\epsilon g_{\alpha }dp-{\nabla }_{x}V\cdot {\int }_{{\mathbb{T}}^d}u{g}_{\alpha }dp=0.$$

$$J_n={\int }_{{\mathbb{T}}^d}ugdp,J_E={\int }_{{\mathbb{T}}^d}u\epsilon gdp.$$

$$\alpha {\partial }_t({\mathcal{F}}_{f_{\alpha }}+\alpha g_{\alpha })+u\cdot {\nabla }_x({\mathcal{F}}_{f_{\alpha }}+\alpha g_{\alpha })+{\nabla }_{x}V\cdot {\nabla }_p({\mathcal{F}}_{f_{\alpha }}+\alpha g_{\alpha })=-\frac{g_{\alpha }}{{\tau }_{\alpha }},$$

$$u\cdot {\nabla }_x{\mathcal{F}}_f+{\nabla }_{x}V\cdot {\nabla }_p{\mathcal{F}}_f=-\frac{g}{\tau }.$$

$${\nabla }_x{\mathcal{F}}_f={\mathcal{F}}_f(1-\eta {\mathcal{F}}_f)({\nabla }_x{\lambda }_0+\epsilon {\nabla }_x{\lambda }_1),{\nabla }_p{\mathcal{F}}_f={\mathcal{F}}_f(1-\eta {\mathcal{F}}_f)u{\lambda }_1,$$

$$g=-\tau{\mathcal{F}}_{f}(1-\eta{\mathcal{F}}_{f})\big{(}u\cdot\nabla_{x}\lambda_{0}+\varepsilon u\cdot\nabla_{x}\lambda_{1}+\lambda_{1}\nabla_{x}V\cdot u\big{)}.$$

$$\mathcal{F}(\lambda )=\frac{1}{\eta +\exp (-{\lambda }_0-{\lambda }_1\epsilon (p))},\lambda =({\lambda }_0,{\lambda }_1)\in {\mathbb{R}}^2,p\in {\mathbb{T}}^d.$$

Full text

The authors acknowledge partial support from the Austrian Science Fund (FWF), grants P27352, P30000, and W1245.

1. Introduction

An optical lattice is a spatially periodic structure that is formed by interfering optical laser beams. The interference produces an optical standing wave that may trap neutral atoms [4]. The lattice potential mimics the crystal lattice in a solid, while the trapped atoms mimic the valance electrons in a solid state crystal. In contrast to solid materials, it is easily possible to adjust the geometry and depth of the potential of an optical lattice. Another advantage is that the dynamics of the atoms, if cooled down to a few nanokelvin, can be followed on the time scale of milliseconds. Therefore, optical lattices are ideal systems to study physical phenomena that are difficult to observe in solid crystals. Moreover, they are promising candidates to realize quantum information processors [14] and extremely precise atomic clocks [2].

The dynamics of ultracold fermionic clouds in an optical lattice can be modeled by the Fermi-Hubbard model with a Hamiltonian that is a result of the lattice potential created by interfering laser beams and short-ranged collisions [11]. In the semi-classical picture, the interplay between diffusive and ballistic regimes can be described by a Boltzmann transport equation [13], which is able to model qualitatively the observed cloud shapes [23].

In this paper, we investigate moment equations which are formally derived from a Boltzmann equation in the diffusive regime. The motivation of our work is the observation that in the (relative) high-temperature limit, the lowest-order diffusion approximation of the Boltzmann equation leads to a logarithmic diffusion equation [20] which has nonphysical properties in the whole space (for instance, it loses mass). Our aim is to derive the next-order approximation, leading to energy-transport equations for the particle and energy densities, and to analyze and simulate the resulting system of degenerate parabolic equations under periodic boundary conditions [20].

The starting point is the (scaled) Boltzmann equation for the distribution function $f(x,p,t)$,

[TABLE]

where $x\in{\mathbb{R}}^{d}$ is the spatial variable, $p$ is the crystal momentum defined on the $d$-dimensional torus ${\mathbb{T}}^{d}$ with unit measure, and $t>0$ is the time. The velocity $u$ is defined by $u(p)=\nabla_{p}\varepsilon(p)$ with the energy $\varepsilon(p)$, $V(x,t)$ is the potential, and $Q(f)$ is the collision operator. Compared to the semiconductor Boltzmann equation, there are two major differences.

First, the band energy $\varepsilon(p)$ is given by the periodic dispersion relation

[TABLE]

The constant $\varepsilon_{0}$ is a measure for the tunneling rate of a particle from one lattice site to a neighboring one. In semiconductor physics, usually a parabolic band structure is assumed, $\varepsilon(p)=\frac{1}{2}|p|^{2}$ [15]. This formula also appears in kinetic gas theory as the (microscopic) kinetic energy. The band energy in optical lattice is bounded, while it is unbounded when $\varepsilon(p)=\frac{1}{2}|p|^{2}$. This has important consequences regarding the integrability of the equilibrium distribution (see below).

Second, the potential $V$ is given by $V=-U_{0}n$, where $n=\int_{{\mathbb{T}}^{d}}fdp$ is the particle density and $U_{0}>0$ models the strength of the on-site interaction between spin-up and spin-down components [23]. In semiconductor physics, $V$ is the electric potential which is a given function or determined self-consistently from the (scaled) Poisson equation [15]. The definition $V=-U_{0}n$ leads to unexpected “degeneracies” in the moment equations, see the discussion below.

The collision operator is given as in [23] by the relaxation-time approximation

[TABLE]

where $\tau>0$ is the relaxation time and ${\mathcal{F}}_{f}$ is determined by minimizing the free energy for fermions associated to (1) under the constraints of mass and energy conservation (see Section 2 for details), leading to

[TABLE]

where $(\lambda_{0},\lambda_{1})$ are the Lagrange multipliers resulting from the mass and energy constraints. For $\eta=1$, we obtain the Fermi-Dirac distribution, while for $\eta=0$, ${\mathcal{F}}_{f}$ equals the Maxwell-Boltzmann distribution. We may consider ${\mathcal{F}}_{f}$ as a function of $(\lambda_{0},\lambda_{1})$ and write ${\mathcal{F}}(\lambda_{0},\lambda_{1};p)=[\eta+\exp((-\lambda_{0}-\lambda_{1}\varepsilon(p))]^{-1}$.

The variable $\lambda_{1}$ can be interpreted as the negative inverse (absolute) temperature, while $\lambda_{0}$ is related to the so-called chemical potential [15]. Since the band energy is bounded, the equilibrium ${\mathcal{F}}_{f}$ is integrable even when $\lambda_{1}>0$, which means that the absolute temperature may be negative. In fact, negative absolute temperatures can be realized in experiments with cold atoms [22]. Negative temperatures occur in equilibrated (quantum) systems that are characterized by an inverted population of energy states. The thermodynamical implications of negative temperatures are discussed in [21].

In the following, we detail the main results of the paper.

Formal derivation and entropy structure

Starting from the Boltzmann equation (1), we derive formally moment equations in the limit of large times and dominant collisions. More precisely, the particle density $n=\int_{{\mathbb{T}}^{d}}{\mathcal{F}}dp$ and energy density $E=\int_{{\mathbb{T}}^{d}}{\mathcal{F}}\varepsilon dp$ solve the energy-transport equations

[TABLE]

for $x\in{\mathbb{R}}^{d}$, $t>0$, where the particle and energy current densities are given by

[TABLE]

and the diffusion coefficients depend nonlocally on ${\mathcal{F}}$ and hence on $(\lambda_{0},\lambda_{1})$; see Proposition 1. The structure of system (3)-(4) is similar to the semiconductor case [17] but the diffusion coefficients $D_{j}$ are different. For $V=-U_{0}n$, the Joule heating term $J_{n}\cdot\nabla V$ contains the squared gradient $|\nabla n|^{2}$, while in the semiconductor case it contains $|\nabla V|^{2}$ which is generally smoother than $|\nabla n|^{2}$.

System (3)-(4) possesses a formal gradient-flow or entropy structure. Indeed, the entropy $H$, defined in Section 2.2, is nonincreasing in time,

[TABLE]

see Proposition 3. Here, the functions $\mu_{0}=\lambda_{0}+\lambda_{1}V$ and $\mu_{1}=\lambda_{1}$ are called the dual entropy variables, and the coefficients $L_{ij}$ are defined in (20). In the dual entropy variables, the potential terms are eliminated, leading to the “symmetric” problem

[TABLE]

where the matrix $(L_{ij})$ is symmetric and positive definite. This formal gradient-flow structure allows for the development of an existence theory but only for uniformly positive definite diffusion matrices [10]. A general existence result (including electric potentials) is still missing.

A further major difficulty comes from the fact that the system possesses certain “degeneracies” in the mapping $(n,E)\mapsto\mu=(\mu_{0},\mu_{1})$ and the entropy production $-dH/dt$. For instance, the determinant of the Jacobi matrix $\partial(n,E)/\partial\mu$ may vanish at certain points. Such a situation also occurs for the semiconductor energy-transport equations but only at the boundary of the domain of definition (namely at $E=0$). In the present situation, the degeneracy may occur at points in the interior of the domain of definition. In view of these difficulties, an analysis of the general energy-transport model (3)-(4) is currently out of reach. This motivates our approach to introduce a simplified model.

Analysis of high-temperature energy-transport models

We show the existence of weak solutions to a simplified energy-transport model. It is argued in [20] that the temperature is large (relative to the nanokelvin scale) in the center of the atomic cloud for long times. Therefore, we simplify (3)-(4) by performing the high-temperature limit. For high temperatures, the relaxation time may be approximated by $\tau(n)=\tau_{0}/(n(1-n))$ [23, Suppl.]. As $\theta=-1/\lambda_{1}$ can be interpreted as the temperature, the high-temperature limit corresponds to the limit $\lambda_{1}\to 0$. Expanding ${\mathcal{F}}(\lambda_{0},\lambda_{1})$ around $(\lambda_{0},0)$ up to zeroth order leads to the diffusion equation (see Section 3)

[TABLE]

In the case $\eta=0$, we obtain the logarithmic diffusion equation $\partial_{t}n=\tau_{0}\Delta\log n$ which predicts a nonphysical behavior. Indeed, in two space dimensions, it can be shown that the particle number is not conserved and the unique smooth solution exists for finite time only; see, e.g., [9, 24]. We expect a similar behavior when $\eta>0$. This motivates us to compute the next-order expansion. It turns out that at first order and with $V=-U_{0}n$, $(n,E)$ is solving the (rescaled) energy-transport equations

[TABLE]

where $U=U_{0}/(2d\varepsilon_{0}^{2})$ and $W=1-UE$ is the “reverted” energy. The case $W=0$ corresponds to the maximal energy $E=1/U_{0}$. Taking into account the periodic lattice structure, we solve (6)-(7) on the torus ${\mathbb{T}}^{d}$, together with the initial conditions $n(0)=n^{0}$, $W(0)=W^{0}$ in ${\mathbb{T}}^{d}$.

The structure of the diffusion equation (6) is similar to (5), but the diffusion coefficient contains $W$ as a factor, adding a degeneracy to the singular logarithmic diffusion equation. It is an open problem whether this factor removes the unphysical behavior of the solution to (5) in ${\mathbb{R}}^{d}$. We avoid this problem by solving (6)-(7) in a bounded domain and by looking for strictly positive particle densities. Is is another open problem to prove the existence of solutions to (6)-(7) in the whole space.

Because of the squared gradient term in (7), the energy $E$ (or $W$) is not conserved but the total energy $W_{\rm tot}=W-(U/2)n^{2}$. In fact, in terms of $W_{\rm tot}$, the squared gradient term is eliminated,

[TABLE]

Unfortunately, this formulation does not help for the analysis since the treatment of $\partial_{t}(n^{2})=2n\partial_{t}n$ is delicate as $\partial_{t}n$ lies in the dual space $H^{1}({\mathbb{T}}^{d})^{\prime}$ but $n$ is generally not an element of $H^{1}({\mathbb{T}}^{d})$ because of the degeneracy (we have only $W^{1/2}\nabla n\in L^{2}({\mathbb{T}}^{d})$).

The analysis of system (6)-(7) is very challenging since the first equation is degenerate in $W$, and the second equation contains a quadratic gradient term. In the literature, there exist existence results for degenerate equations with quadratic gradient terms [8, 12], but the degeneracy is of porous-medium type. A more complex degeneracy was investigated in [7]. In our case, the degeneracy comes from another variable, which is much more delicate to analyze.

Related problems appear in semiconductor energy-transport theory, but only partial results have been obtained so far. Let us review these results. The existence of stationary solutions to (3) with the current densities

[TABLE]

close to the constant equilibrium has been shown in [1]. The idea is that in such a situation, the temperature $\theta$ is strictly positive which removes the degeneracy in the term $\nabla(n\theta)$. The parabolic system was investigated in [18, 19], and the global existence of weak solutions was shown without any smallness condition but for a simplifed energy equation. Again, the idea was to prove a uniform positivity bound for the temperature, which removes the degeneracy. A more general result (but without electric potential) was achieved in [25] for the system

[TABLE]

in a bounded domain, where $0<\alpha<1$. The global existence of weak solutions to the corresponding initial-boundary-value problem was proved. Again, the idea is a positivity bound for $\theta$ but this bound required a nontrivial cut-off procedure and several entropy estimates.

In this paper, we make a step forward in the analysis of nonlinear parabolic systems with nonstandard degeneracies by solving (6)-(7) without any positive lower bound for $W$. Since $W$ may vanish, we can expect a gradient estimate for $n$ only on $\{W>0\}$. Although the quadratic gradient term also possesses $W$ as a factor, the treatment of this term is highly delicate, because of low time regularity. Therefore, we present a result only for a time-discrete version of (6)-(7), namely for its implicit Euler approximation

[TABLE]

for $x\in{\mathbb{T}}^{d}$, where $\triangle t>0$ and $(n^{k-1},W^{k-1})$ are given functions. We show the existence of a weak solution $(n^{k},W^{k})$ satisfying $n^{k}\geq 0$, $W^{k}\geq 0$ and $W^{k}n^{k}$, $W^{k}\in H^{1}({\mathbb{T}}^{d})$; see Theorem 9. In one space dimension and under a smallness assumption on the variance of $W^{k-1}$ and $n^{k-1}$, the strict positivity of $W^{k}$ can be proved; see Theorem 11.

The existence proof is based on the solution of a regularized and truncated problem by means of the Leray-Schauder fixed-point theorem. Standard elliptic estimates provide bounds uniform in the approximation parameters. The key step is the proof of the strong convergence of the gradient of the particle density. For this, we show a general result for degenerate elliptic problems; see Proposition 8. This result seems to be new. Standard results in the literature need the ellipticity of the differential operator [5]. Unfortunately, we are not able to perform the limit $\triangle t\to 0$ since some estimates in the proof of Proposition 8 are not uniform in $\triangle t$; also see Remark 10 for a discussion.

Numerical simulations

The time-discrete system (9)-(10) is discretized by finite differences in one space dimension and solved in an semi-implicit way. The large-time behavior exhibits an interesting phenomenon. If the initial energy $W^{0}$ is constant and sufficiently large, the solution $(n(t),W(t))$ converges to a constant steady state. However, if the constant $W^{0}$ is too small, the stationary particle density is nonconstant. In both cases, the time decay is exponential fast, but the decay rate becomes smaller for smaller constants $W^{0}$ since the diffusion coefficient in (6) is smaller too.

The paper is organized as follows. Section 2 is devoted to the formal derivation of the general energy-transport model and its entropy structure, similar to the semiconductor case [3]. The high-temperature expansion is performed in Section 3, leading to the energy-transport system (6)-(7). The strong convergence of the gradients is shown in Section 4. In Section 5 the existence result is stated and proved. The numerical simulations are presented in Section 6, and the Appendix is concerned with the calculation of some integrals involving the velocity $u(p)$ and energy $\varepsilon(p)$.

2. Formal derivation and entropy structure

2.1. Derivation from a Boltzmann equation

We consider the following semiclassical Boltzmann transport equation for the distribution function $f(x,p,t)$ in the diffusive scaling:

[TABLE]

where $\alpha>0$ is the Knudsen number [3], $(x,p)$ are the phase-space variables (space and crystal momentum), and $t>0$ is the time. We recall that the velocity equals $u(p)=\nabla_{p}\varepsilon(p)$, where the energy $\varepsilon(p)$ is given by (2). The potential $V_{\alpha}$ is defined by $V_{\alpha}=-U_{0}n_{\alpha}$. In the physical literature [23], the collision operator $Q_{\alpha}$ is given by the relaxation-time approximation

[TABLE]

where the function ${\mathcal{F}}_{f}$ is determined by maximizing the free energy (18) associated to (11) under the constraints

[TABLE]

which express mass and energy conservation during scattering events. The solution of this problem is given by

[TABLE]

where $\lambda_{0}$ and $\lambda_{1}$ are the Lagrange multipliers and $\eta\geq 0$ is a parameter which may take the values $\eta=0$ (Maxwell-Boltzmann statistics) or $\eta=1$ (Fermi-Dirac statistics). The relaxation time $\tau_{\alpha}\geq 0$ generally depends on the particle density but at this point we do not need to specifiy the dependence.

We show the following result.

Proposition 1** (Derivation).**

Let $f_{\alpha}$ be a (smooth) solution to the Boltzmann equation (11). We assume that the formal limits $f=\lim_{\alpha\to 0}f_{\alpha}$, $g=\lim_{\alpha\to 0}(f_{\alpha}-{\mathcal{F}}_{f_{\alpha}})/\alpha$, and $\tau=\lim_{\alpha\to 0}\tau_{\alpha}$ exist. Then the particle and energy densities

[TABLE]

are solutions to (3)-(4), and the diffusion coefficients $D_{ij}=(D_{ij}^{k\ell})\in{\mathbb{R}}^{d\times d}$ are defined by

[TABLE]

The proof of the proposition is similar to those of Propositions 1 and 2 in [17]. For the convenience of the reader, we present the (short) proof.

Proof.

To derive the balance equations, we multiply the Boltzmann equation (11) by $1$ and $\varepsilon$, respectively, and integrate over ${\mathbb{T}}^{d}$:

[TABLE]

The integrals involving the collision operator vanish in view of mass and energy conservation; see (12). We have integrated by parts in the last integral on the left-hand side of (14). Next, we insert the Chapman-Enskog expansion $f_{\alpha}={\mathcal{F}}_{f_{\alpha}}+\alpha g_{\alpha}$ (which in fact defines $g_{\alpha}$) in (13)-(14) and observe that the function $p\mapsto u(p)\varepsilon(p)^{j}{\mathcal{F}}_{f_{\alpha}}(p)$ is odd for any $j\in{\mathbb{N}}_{0}$ such that its integral over ${\mathbb{T}}^{d}$ vanishes. This leads to

[TABLE]

Passing to the formal limit $\alpha\to 0$ gives the balance equations (3) with

[TABLE]

To specify the current densities, we insert the Chapman-Enskog expansion in (11),

[TABLE]

and perform the formal limit $\alpha\to 0$,

[TABLE]

A straightforward computation shows that

[TABLE]

and inserting this into (16) gives an explicit expression for $g$:

[TABLE]

Therefore, the current densities (15) lead to (4). This finishes the proof. ∎

In the following we write ${\mathcal{F}}_{f}={\mathcal{F}}(\lambda)$, where

[TABLE]

Proposition 2** (Diffusion matrix).**

The diffusion matrix ${\mathcal{D}}=(D_{ij})\in{\mathbb{R}}^{2d\times 2d}$ is symmetric and positive definite.

Proof.

The proof is similar to Proposition 3 in [17]. Let $z=(w,y)\in{\mathbb{R}}^{2d}$ with $w$, $y\in{\mathbb{R}}^{d}$. Then

[TABLE]

Since $D_{ij}^{k\ell}$ is symmetric in $(i,j)$ and $(k,\ell)$, the symmetry of ${\mathcal{D}}$ is clear. ∎

2.2. Entropy structure

The entropy structure of (3)-(4) follows from the abstract framework presented in [17]. In the following, we make this framework explicit. First, we introduce the entropy

[TABLE]

The entropy density $h$ can be reformulated as

[TABLE]

The following result shows that the entropy is nonincreasing in time.

Proposition 3** (Entropy structure).**

It holds that

[TABLE]

where $\mu_{0}=\lambda_{0}+\lambda_{1}V$ and $\mu_{1}=\lambda_{1}$ are the so-called dual entropy variables and

[TABLE]

Proof.

Identity (18) implies that

[TABLE]

and consequently,

[TABLE]

Therefore, using (3)-(4) and integration by parts,

[TABLE]

which proves the identity in (19). Using the positive definiteness of ${\mathcal{D}}$, a computation shows that $(L_{ij})$ is positive definite too, and the inequality in (19) follows. ∎

2.3. Singularities and degeneracies in the energy-transport system

We denote by $n$ and $E$ the particle and energy densities depending on the dual entropy variable $\mu=(\mu_{0},\mu_{1})=(\lambda_{0}+V\lambda_{1},\lambda_{1})$. We have the (implicit) formulation

[TABLE]

Lemma 4**.**

Let $\omega_{i}(\mu):=\int_{{\mathbb{T}}^{d}}{\mathcal{F}}(1-\eta{\mathcal{F}})\varepsilon(p)^{i}dp$, $i\in{\mathbb{N}}_{0}$. Then

[TABLE]

Proof.

We differentiate

[TABLE]

This gives after a rearrangement

[TABLE]

In a similar way, we obtain

[TABLE]

and with

[TABLE]

the conclusion follows. ∎

Since $\mu_{1}$ can be positive and $\omega_{0}>0$, the expression $1-U_{0}\mu_{1}\omega_{0}$ may vanish, so the determinant of $\partial(n,E)/\partial\mu$ may be not finite. Moreover, the numerator of the determinant may vanish, and the function $\mu\mapsto(n,E)$ may be not invertible. This is made more explicit in the following remark.

Remark 5** (Case $\mu=0$).**

In the Maxwell-Boltzmann case, we can make the numerator of the determinant in Lemma 4 explicit. Indeed, it is clear that $\omega_{0}=n$ and $\omega_{1}=E$. For the computation of $\omega_{2}$, we observe first that

[TABLE]

where, by symmetry,

[TABLE]

Now, we have

[TABLE]

Expanding the exponentials in ${\mathcal{F}}$ and using (21)-(22), we find that

[TABLE]

and this expression is independent of $i\neq j$. For $i=j$, we use the identity $\sin(2\pi p_{i}){\mathcal{F}}=(\partial{\mathcal{F}}/\partial p_{i})/(4\pi\varepsilon_{0}\mu_{1})$ and integration by parts:

[TABLE]

Summing over $i=1,\ldots,n$ gives

[TABLE]

We conclude that $\omega_{2}=4\varepsilon_{0}^{2}dn-E/\mu_{1}+(d-1)E^{2}/(dn)$ and consequently,

[TABLE]

For certain values of $\mu_{1}$ or $(n,E)$, this expression may vanish such that $\det\partial(n,E)/\partial\mu=0$ at these values. This shows that the relation between $(n,E)$ and $\mu$ needs to treated with care. ∎

Remark 6** (Degeneracy in the entropy production).**

A tedious computation, detailed in [6, Chapter 8.4], shows that the entropy production can be written as

[TABLE]

where $g_{i}(\lambda)$, $i=1,2,3$, are functions depending on $\omega_{i}$, defined in Lemma 4, and on

[TABLE]

The above formula shows that we lose the gradient estimate if $1-U_{0}\mu_{1}\omega_{0}=0$. ∎

Remark 7** (Comparison with the semiconductor case).**

For the semiconductor energy-transport equations in the parabolic band approximation, we do not face the singularities and degeneracies occuring in the model for optical lattices. Indeed, let the potential $V$ be given (to simplify). According to Example 6.8 in [15], we have

[TABLE]

Then

[TABLE]

which is nonzero as long as $E>0$. Furthermore, by Remark 8.12 in [15], it holds that $\omega_{0}=n$, $\omega_{1}=E$, and $\omega_{2}=5E^{2}/(3n)$, and so

[TABLE]

This expression is degenerate only at the boundary of the domain of definition (i.e. at $E=0$). Often, such kind of degeneracies may be handled; an important example is the porous-medium equation. In the case of optical lattices, the degeneracy may occur in the interior of the domain of definition, which is much more delicate. ∎

3. High-temperature expansion

The Lagrange multiplier $\lambda_{1}$ is interpreted as the negative inverse temperature, so high temperatures correspond to small values of $|\lambda_{1}|$. In this section, we perform a high-temperature expansion of (3)-(4), i.e., we expand ${\mathcal{F}}(\lambda)$ around $(\lambda_{0},0)$ for small $|\lambda_{1}|$ up to first order. Our ansatz is

[TABLE]

3.1. Zeroth-order expansion

At zeroth-order, we have by (15), (16), and using formula (67) from the appendix,

[TABLE]

Therefore, up to order $O(\lambda_{1})$, we infer that

[TABLE]

At high temperature, the relaxation time depends on the particle density in a nonlinear way, $\tau=\tau_{0}/(n(1-\eta n))$ [20]. At low densities, i.e. $\eta=0$, we obtain the logarithmic diffusion equation

[TABLE]

where $\varepsilon_{1}=\frac{1}{2}\tau_{0}(4\pi\varepsilon_{0})^{2}$. We already mentioned in the introduction that the (smooth) solution to this equation in two space dimensions loses mass, which is unphysical. Therefore, we compute the next-order expansion.

3.2. First-order expansion

We calculate, using (66),

[TABLE]

Therefore, by (23),

[TABLE]

and (16) and $\nabla_{p}\varepsilon=u$ give, up to order $O(\lambda_{1}^{2})$,

[TABLE]

Then, by (15)-(16) and taking into account (67)–(69), we infer that, again up to first order,

[TABLE]

Therefore, the first-order expansion leads to

[TABLE]

We rescale the time by $t_{s}=(8\pi^{2}\varepsilon_{0}^{2}\tau_{0})t$ and introduce $U=U_{0}/(2d\varepsilon_{0}^{2})$ and $W=1-UE$. Then, writing again $t$ instead of $t_{s}$, system (24)-(25) becomes

[TABLE]

The existence of weak solutions to a time-discrete version of (26), together with periodic boundary conditions, is shown in Section 5.

4. A strong convergence result for the gradient

The key tool of the existence analysis of Section 5 is the following result on the strong convergence of the gradients of certain approximate solutions for the following equation. Let $\triangle t>0$, $\overline{n}\in L^{\infty}(\Omega)$, $y\in L^{\infty}(\Omega)\cap H^{1}(\Omega)$, and $\Psi\in C^{1}({\mathbb{R}})$. We consider the equation

[TABLE]

for $n\in L^{\infty}(\Omega)$ such that $y\Psi(n)\in H^{1}(\Omega)$.

Proposition 8**.**

Let $\Omega\subset{\mathbb{R}}^{d}$ be a bounded domain, $\varepsilon>0$, $\triangle t>0$, let $\overline{n}\in L^{\infty}(\Omega)$ be such that $\overline{n}\geq 0$ in $\Omega$, and let $\Psi\in C^{1}({\mathbb{R}})$ satisfy $\Psi^{\prime}>0$. Let $(y_{\varepsilon})$ be a bounded sequence in $H^{1}(\Omega)$ satisfying $y_{\varepsilon}\geq C(\varepsilon)>0$ for some $C(\varepsilon)>0$, $y_{\varepsilon}\to y$ strongly in $L^{2}(\Omega)$ and weakly in $H^{1}(\Omega)$ as $\varepsilon\to 0$. Furthermore, let $n_{\varepsilon}\in L^{2}(\Omega)$ with $\Psi(n_{\varepsilon})\in H^{1}(\Omega)$ be a weak solution to

[TABLE]

for all $\phi\in H^{1}(\Omega)$. Then there exist a function $n\in L^{\infty}(\Omega)$ such that $y\Psi(n)\in H^{1}(\Omega)$, being a weak solution of (27), and a subsequence of $(n_{\varepsilon})$, which is not relabeled, such that, as $\varepsilon\to 0$,

[TABLE]

Proof.

Step 1. First, we derive some uniform bounds. Set $M=\|\overline{n}\|_{L^{\infty}(\Omega)}$. Taking $(\Psi(n_{\varepsilon})-\Psi(M))_{+}=\max\{0,\Psi(n_{\varepsilon})-\Psi(M)\}$ as a test function in (28), we find that

[TABLE]

Since $-(\overline{n}-M)(\Psi(n_{\varepsilon})-\Psi(M))_{+}\geq 0$, it follows that

[TABLE]

and hence, $n_{\varepsilon}\leq M$ in $\Omega$. In a similar way, using $\Psi(n_{\varepsilon})_{-}=\min\{0,\Psi(n_{\varepsilon})\}$ as a test function and using $\overline{n}\geq 0$, we infer that $n_{\varepsilon}\geq 0$. This shows that $(n_{\varepsilon})$ is bounded in $L^{\infty}(\Omega)$. Hence, there exists a subsequence which is not relabeled such that, as $\varepsilon\to 0$,

[TABLE]

for some function $n\in L^{\infty}(\Omega)$. Since $\Phi^{\prime}$ is positive and continuous on ${\mathbb{R}}$, the boundedness of $(n_{\varepsilon})$ in $L^{\infty}(\Omega)$ implies that there exists a constant $C>0$ such that $1/C\leq\Psi^{\prime}(n_{\varepsilon})\leq C$ for all $\varepsilon>0$. Next, we choose the test function $\phi=\Psi(n_{\varepsilon})$ in (28) and use $\Psi(n_{\varepsilon})\leq\Psi(M)$ to find that

[TABLE]

We deduce that $(y_{\varepsilon}^{1/2}\nabla\Psi(n_{\varepsilon}))$ is bounded in $L^{2}(\Omega)$ and, for a subsequence,

[TABLE]

for some function $\xi\in L^{2}(\Omega)$. Now, let $\phi$ be a smooth test function. Then we can take the limit, for a subsequence, in (28) and obtain

[TABLE]

Note that this equation holds also for all $\phi\in H^{1}(\Omega)$.

Step 2. As $(y_{\varepsilon})$ is strongly converging in $L^{2}(\Omega)$, we deduce from (29) that $y_{\varepsilon}^{1/2}(n_{\varepsilon}-n)\rightharpoonup 0$ weakly in $L^{2}(\Omega)$. We claim that this convergence is even strong. Indeed, the sequence

[TABLE]

is uniformly bounded in $L^{2}(\Omega)$. Thus, $(y_{\varepsilon}n_{\varepsilon}^{2})$ is bounded in $H^{1}(\Omega)$ and by compactness, for a subsequence, $y_{\varepsilon}n_{\varepsilon}^{2}\to\zeta\geq 0$ strongly in $L^{2}(\Omega)$ or $y_{\varepsilon}^{1/2}n_{\varepsilon}\to\zeta^{1/2}$ strongly in $L^{4}(\Omega)$. The strong convergence of $(y_{\varepsilon}^{1/2})$ in $L^{2}(\Omega)$ and the weak* convergence of $(n_{\varepsilon})$ in $L^{\infty}(\Omega)$ imply that $y_{\varepsilon}^{1/2}n_{\varepsilon}\rightharpoonup y^{1/2}n$ weakly in $L^{2}(\Omega)$. Therefore, $\zeta^{1/2}=y^{1/2}n$ and

[TABLE]

This proves the claim.

Step 3. The next goal is to show that

[TABLE]

Taking into account (32), the strong convergence of $(y_{\varepsilon}^{1/2})$ in $L^{2}(\Omega)$, and the $L^{\infty}$ bound for $(n_{\varepsilon})$, it follows that

[TABLE]

converges to zero. Thus, for a subsequence, $y^{1/2}(n_{\varepsilon}-n)\to 0$ a.e. in $\Omega$ and consequently, $n_{\varepsilon}-n\to 0$ a.e. in $\{y>0\}$. Then the a.e. pointwise convergence $\mathrm{1}_{\{y>0\}}(n_{\varepsilon}-n)\to 0$ and the dominated convergence theorem show (33).

Step 4. We wish to identify $\xi$ in (30) and (31). The bounds for $(\nabla y_{\varepsilon})$ and $(y_{\varepsilon}\nabla\Psi(n_{\varepsilon}))$ in $L^{2}(\Omega)$ show that $(\nabla(y_{\varepsilon}\Psi(n_{\varepsilon})))$ is bounded in $L^{2}(\Omega)$ and so, $(y_{\varepsilon}\Psi(n_{\varepsilon}))$ is bounded in $H^{1}(\Omega)$. By compactness, for a subsequence, $\nabla(y_{\varepsilon}\Psi(n_{\varepsilon}))\rightharpoonup\nabla(y\Psi(n))$ weakly in $L^{2}(\Omega)$ and $y_{\varepsilon}\Psi(n_{\varepsilon})\to y\Psi(n)$ strongly in $L^{2}(\Omega)$. We can identify the limit since $y_{\varepsilon}\to y$ strongly in $L^{2}(\Omega)$ and $\Psi(n_{\varepsilon})\rightharpoonup^{*}\theta$ weakly* in $L^{\infty}(\Omega)$ with $\theta=\Psi(n)$ in $\{y>0\}$ lead to $y_{\varepsilon}\Psi(n_{\varepsilon})\rightharpoonup y\Psi(n)$ weakly in $L^{2}(\Omega)$. Moreover, since $\Psi(n_{\varepsilon})\to\Psi(n)$ strongly in $L^{2}(\Omega)$ and $\nabla y_{\varepsilon}\to\nabla y$ weakly in $L^{2}(\Omega)$, we infer that

[TABLE]

Here, we have used additionally that $(y_{\varepsilon}\nabla\Psi(n_{\varepsilon}))$ is bounded in $L^{2}(\Omega)$ and that $L^{1}(\Omega)$ is dense in $L^{2}(\Omega)$. Similarly as above, we deduce that

[TABLE]

converges to zero. Therefore, $y^{1/2}\xi=\nabla(y\Psi(n))-\Psi(n)\nabla y$ in $\Omega$.

Step 5. We obtain from (34) that

[TABLE]

and consequently,

[TABLE]

applying the dominated convergence theorem to the last integral. Furthermore, using the test function $\phi=y(\Psi(n_{\varepsilon})-\Psi(n))$ in (28),

[TABLE]

By (33), we have $\|y(\Psi(n_{\varepsilon})-\Psi(n))\|_{L^{2}(\Omega)}\to 0$ and by (32), $y_{\varepsilon}(\Psi(n_{\varepsilon})-\Psi(n))^{2}\to 0$ in $\Omega$ for a subsequence. Then, by dominated convergence, $\int_{\Omega}y_{\varepsilon}(\Psi(n_{\varepsilon})-\Psi(n))^{2}|\nabla y|^{2}dx\to 0$. We have proved that

[TABLE]

Subtracting (36) from (37), we conclude that

[TABLE]

Taking into account this convergence and (34), it follows again by the dominated convergence theorem that

[TABLE]

This shows the first part of the proposition.

Step 6. It remains to show that the limit $n$ solves (27). Let $\phi\in W^{1,\infty}(\Omega)$. Since (a subsequence of) $(n_{\varepsilon})$ converges weakly* to $n$ in $L^{\infty}(\Omega)$, we have

[TABLE]

Furthermore,

[TABLE]

By Step 5, the first integral converges to

[TABLE]

while the second integral converges to zero since

[TABLE]

We conclude that (27) holds in the weak sense for test functions in $W^{1,\infty}(\Omega)$ but a density argument shows that it is sufficient to take test functions in $H^{1}(\Omega)$. This finishes the proof. ∎

5. Existence of solutions to the high-temperature model

We prove the existence of weak solutions to (26) in ${\mathbb{T}}^{d}$. We recall the definition of the total (“reverted”) energy

[TABLE]

and introduce the total variance

[TABLE]

The main result is as follows.

Theorem 9** (Existence of weak solutions).**

Let $\triangle t>0$, $U>0$, $\eta\in(0,1]$, $0<\delta<1/(1+\eta)$ and let

[TABLE]

Then there exists a weak solution $(n^{k},W^{k})$ to (9)-(10) in the following sense: It holds $\delta\leq n^{k}\leq\|n^{k-1}\|_{L^{\infty}({\mathbb{T}}^{d})}\leq(1-\delta)/\eta$, $0\leq W^{k}\leq\|W^{k-1}\|_{L^{\infty}({\mathbb{T}}^{d})}$ in ${\mathbb{T}}^{d}$, $W^{k}n^{k}$, $W^{k}\in H^{1}({\mathbb{T}}^{d})$, as well as

[TABLE]

for all $\phi_{0}\in H^{1}({\mathbb{T}}^{d})$ and $\phi_{1}\in H^{1}({\mathbb{T}}^{d})\cap L^{\infty}({\mathbb{T}}^{d})$, where $g(n^{k})=n^{k}(1-\eta n^{k})$. For this solution, the following monotonicity properties hold:

[TABLE]

where $W_{\rm tot}^{k}$ and $V^{k}$ are defined in (38) and (39), respectively. Moreover, if

[TABLE]

holds then $W^{k}\not\equiv 0$.

Remark 10** (Comments).**

1. The existence result holds true for more general functions $g(n)$ under the assumption that $g(n)$ is strictly positive for $\delta\leq n\leq(1-\delta)/\eta$.

2. One may interpret $W^{k}$ as a “renormalized” solution since we need test functions of the form $W^{k}\phi_{1}$ in order to avoid vacuum sets $W^{k}=0$. Such an idea has been used, for instance, for the compressible quantum Navier-Stokes equations to avoid vacuum sets in the particle density [16]. Test functions of the type $W^{k}\phi$ allow for the trivial solution $n^{k}=n^{k-1}$ and $W^{k}=0$ but assumption (43) excludes this situation. It means that no constant steady state with $W^{k}=0$ exists if the variance of $n^{k-1}$ is small compared to the energy $\int_{{\mathbb{T}}^{d}}W^{k-1}dx$.

3. The second inequality in (42) involves $W^{k-2}$ which makes sense when the equations are solved iteratively, starting from $k=1$. Also (43) can be iterated. Indeed, if (43) holds for $(n^{k-1},W^{k-1})$, the monotonicity property (42) and mass conservation $\int_{{\mathbb{T}}^{d}}n^{k}dx=\int_{{\mathbb{T}}^{d}}n^{k-1}dx$ imply that

[TABLE]

4. We are not able to perform the limit $\triangle t\to 0$. The reason is that we cannot perform the limit in the quadratic gradient term $|\nabla(W^{k}n^{k})-n^{k}\nabla W^{k}|^{2}$, since we cannot prove the strong convergence of $\nabla(W^{k}n^{k})-n^{k}\nabla W^{k}$. Proposition 8 provides such a result for the time-discrete elliptic case. The key step is to show that

[TABLE]

where $n$, $W$ are the (weak) limits of $(n^{k})$, $(W^{k})$, respectively, and $\langle\cdot,\cdot\rangle$ is the dual product between $H^{1}({\mathbb{T}}^{d})^{\prime}$ and $H^{1}({\mathbb{T}}^{d})$. It is possible to show that $\triangle t^{-1}(n^{k}-n^{k-1})$ is bounded in $H^{1}({\mathbb{T}}^{d})^{\prime}$, but the limit $W^{k}n^{k}-Wn\to 0$ strongly in $H^{1}({\mathbb{T}}^{d})$ (more precisely: the limit of the piecewise constant in time construction of $W^{k}n^{k}$ in $L^{2}(0,T;H^{1}({\mathbb{T}}^{d}))$) cannot be expected. ∎

In the one-dimensional case and under the smallness condition (44) below, we can show that $W^{k}$ is positive, which allows us to define the weak solution to (9)-(10) in the standard sense (with test functions $\phi_{1}$ instead of $W^{k}\phi_{1}$). We set

[TABLE]

Theorem 11** (One-dimensional case).**

Let the assumptions of Theorem 9 hold, let $d=1$, $G=\max_{\delta\leq s\leq\|n^{k-1}\|_{L^{\infty}({\mathbb{T}})}}g(s)$, and let $(n^{k},W^{k})$ for $k\geq 0$ be the solution given by Theorem 9. We assume that

[TABLE]

Then $W^{k}$ is strictly positive, $n^{k}\in H^{1}({\mathbb{T}})$, and (9)-(10) hold in the sense of $H^{1}({\mathbb{T}})^{\prime}$.

We proceed to the proof of Theorems 9 and 11. In this section, $\varepsilon$ denotes a positive parameter and not the band energy. Since we are not concerned with the kinetic equations, no notational confusion will occur. Let $\triangle t>0$, $\alpha$, $\gamma$, $\delta$, $\varepsilon>0$ satisfying $\gamma<1$ and $\delta<1/(1+\eta)$. Define the truncations

[TABLE]

and $g_{\delta}(s)=[s]_{\delta}(1-\eta[s]_{\delta})$ for $s\in{\mathbb{R}}$. Then $g_{\delta}$ is continuous and strictly positive. Given $W^{k-1}$, $n^{k-1}\in L^{\infty}({\mathbb{T}}^{d})$ satisfying $\delta\leq n^{k-1}\leq(1-\delta)/\eta$, we solve the regularized and truncated nonlinear problem in ${\mathbb{T}}^{d}$

[TABLE]

Remark 12**.**

Let us explain the approximation (45)-(46). The truncation of $W^{k}$ with parameter $\gamma$ ensures that the coefficients are bounded, while the truncation $g_{\delta}(n^{k})$ with parameter $\delta>0$ guarantees that the denominator is always positive. The regularization parameter $\varepsilon$ gives strict ellipticity for (45), since generally the first term on the right-hand side of (45) without $\varepsilon$ is degenerate. Finally, the approximation of the quadratic gradient term with parameter $\alpha$ avoids regularity issues since it holds $|\nabla n^{k}|^{2}\in L^{1}({\mathbb{T}}^{d})$ only. ∎

5.1. Solution of an approximated problem

First, we prove the existence of solutions to (45)-(46).

Lemma 13** (Existence for the approximated problem).**

There exists a weak solution $(n^{k},$ $W^{k})\in H^{1}({\mathbb{T}}^{d})^{2}$ to (45)-(46).

Proof.

We define the fixed-point operator $S:L^{2}({\mathbb{T}}^{d})^{2}\times[0,1]\to L^{2}({\mathbb{T}}^{d})^{2}$ by $S(n^{*},W^{*};\theta)=(n,W)$, where $(n,W)\in H^{1}({\mathbb{T}}^{d})^{2}$ is the unique solution to the linear problem

[TABLE]

where

[TABLE]

The approximation and truncation ensure that these forms are bounded on $H^{1}({\mathbb{T}}^{d})$. The bilinear forms $a_{0}$ and $a_{1}$ are coercive. By the Lax-Milgram lemma, there exists a unique solution $(n,W)\in H^{1}({\mathbb{T}}^{d})^{2}$ to (47). Thus, the fixed-point operator is well defined (and has compact range). Furthermore, $S(n^{*},W^{*};0)=0$. Standard arguments show that $S$ is continuous. Let $(n,W)$ be a fixed point of $S(\cdot,\cdot;\theta)$, i.e., $(n,W)$ solves (45)-(46) with $(n^{k},W^{k})$ replaced by $(n,W)$. With the test functions $\phi_{0}=n$ and $\phi_{1}=W$ and the inequality

[TABLE]

we find that

[TABLE]

The last integral is nonnegative since $[W]_{\gamma}W\geq 0$. Therefore,

[TABLE]

where $C(\varepsilon)$ and $C(\delta)$ are positive constants independent of $(n,W)$. This provides the necessary uniform bound for all fixed points of $S(\cdot,\cdot;\theta)$. We can apply the Leray-Schauder fixed-point theorem to infer the existence of a fixed point for $S(\cdot,\cdot;1)$, i.e. of a weak solution to (45)-(46). ∎

5.2. Removing the truncation

The following maximum principle holds.

Lemma 14** (Maximum principle).**

Let $(n^{k},W^{k})$ be a weak solution to (45)-(46). Then

[TABLE]

where $\gamma\leq 1/\|W^{k-1}\|_{L^{\infty}({\mathbb{T}}^{d})}$.

Proof.

We choose $(n^{k}-\delta)_{-}=\min\{0,n^{k}-\delta\}$ as a test function in (45):

[TABLE]

Since $-(n^{k-1}-\delta)(n^{k}-\delta)_{-}\geq 0$, this gives

[TABLE]

and hence, $n^{k}\geq\delta$ in ${\mathbb{T}}^{d}$. In a similar way, the test function $(n^{k}-N)_{+}=\max\{0,n^{k}-N\}$ with $N:=\|n^{k-1}\|_{L^{\infty}({\mathbb{T}}^{d})}$ leads to $n^{k}-N\leq 0$ in ${\mathbb{T}}^{d}$. Next, we use $W^{k}_{-}\leq 0$ as a test function in (46):

[TABLE]

We deduce that $W^{k}\geq 0$. The proof of $W^{k}\leq\|W^{k-1}\|_{L^{\infty}({\mathbb{T}}^{d})}\leq 1/\gamma$ is similar, using the test function $(W^{k}-\|W^{k-1}\|_{L^{\infty}({\mathbb{T}}^{d})})_{+}$. ∎

We have shown that $(n^{k},W^{k})$ solves

[TABLE]

where $g(n)=n(1-\eta n)$.

5.3. The limit $\alpha\to 0$

Let $(n_{\alpha}^{k},W_{\alpha}^{k})$ be a weak solution to (48)-(49). We use the test function $n^{k}_{\alpha}$ in (48),

[TABLE]

and the test function $W^{k}_{\alpha}$ in (49),

[TABLE]

which provides immediately uniform $H^{1}({\mathbb{T}}^{d})$ estimates since $g(n^{k}_{\alpha})\geq C(\delta)>0$:

[TABLE]

where the constants are independent of $\alpha$. By compactness, this implies the existence of a subsequence which is not relabeled such that, as $\alpha\to 0$,

[TABLE]

This shows that, maybe for a subsequence, $W^{k}_{\alpha}/g(n^{k}_{\alpha})\to W^{k}/g(n^{k})$ and $1/g(n^{k}_{\alpha})\to 1/g(n^{k})$ a.e. in ${\mathbb{T}}^{d}$, and by dominated convergence, strongly in $L^{2}({\mathbb{T}}^{d})$.

We claim that $n^{k}_{\alpha}\to n^{k}$ strongly in $H^{1}({\mathbb{T}}^{d})$. Let $y_{\alpha}:=(W^{k}_{\alpha}+\varepsilon)/g(n^{k}_{\alpha})$. Then $y_{\varepsilon}\geq\varepsilon/\sup_{s\in[\delta,N]}g(s)>0$, where $N=\|n^{k-1}\|_{L^{\infty}({\mathbb{T}}^{d})}$, and $y_{\alpha}\to y:=W^{k}/g(n^{k})\geq 0$ strongly in $L^{2}(\Omega)$. Thus, $y_{\alpha}^{1/2}\nabla n^{k}_{\alpha}\rightharpoonup y^{1/2}\nabla n^{k}$ weakly in $L^{2}(\Omega)$, and it follows that

[TABLE]

Taking $n^{k}_{\alpha}-n^{k}$ as a test function in (48), we obtain

[TABLE]

Subtraction of these integrals leads to

[TABLE]

Since $y_{\alpha}\geq\varepsilon/\sup_{s\in[\delta,N]}g(s)>0$, this proves the claim. In particular, $1/(1+\alpha|\nabla n^{k}_{\alpha}|^{2})\to 1$ in $L^{2}({\mathbb{T}}^{d})$. From this, we can directly deduce that

[TABLE]

The above convergence results are sufficient to pass to the limit $\alpha\to 0$ in (48)-(49), showing that $(n^{k},W^{k})$ solves

[TABLE]

5.4. The limit $\varepsilon\to 0$

This limit is the delicate part of the proof. We first state a lemma concerning weak and strong convergence.

Lemma 15**.**

Let $(f_{n})$ be a weakly and $(g_{n})$ be a strongly converging sequence in $L^{2}(\Omega)$ which have the same limit. If $|f_{n}(x)|\leq|g_{n}(x)|$ for all $n\in\mathbb{N}$ and a.e. $x\in{\mathbb{T}}^{d}$, then $(f_{n})$ converges strongly in $L^{2}(\Omega)$.

Proof.

Let $f$ denote the weak limit of $(f_{n})$ and $(g_{n})$. Due to the weak lower semi-continuity of the norm,

[TABLE]

Thus, the limes inferior and superior coincide and $\|f_{n}\|_{L^{2}(\Omega)}\to\|f\|_{L^{2}(\Omega)}$. Together with the weak convergence of $(f_{n})$, we deduce the strong convergence. ∎

Let $(n^{k}_{\varepsilon},W^{k}_{\varepsilon})$ be a weak solution to (52)-(53). Inequalities (50) and (51) show the following bounds uniform in $\varepsilon$:

[TABLE]

By compactness, there exists a subsequence (not relabeled) such that, as $\varepsilon\to 0$,

[TABLE]

Again, we need strong convergence for $\nabla n^{k}_{\varepsilon}$. Since equation (52) is degenerate, we obtain a weaker result. For this, let $\Psi\in C^{2}({\mathbb{R}})$ be strictly monotonically increasing and satisfy $\Psi^{\prime}(t)=1/g(t)$ for $\delta\leq t\leq(1-\delta)/\eta$. Thus, we can apply Proposition 8 for $y_{\varepsilon}=W^{k}_{\varepsilon}+\varepsilon$ and $y=W^{k}$ to conclude that, up to a subsequence,

[TABLE]

The latter convergence implies that $n^{k}_{\varepsilon}-n^{k}\to 0$ a.e. in $\{W^{k}>0\}$ and, by dominated convergence,

[TABLE]

Now, let $y_{\varepsilon}=(W^{k}_{\varepsilon}+\varepsilon)/g(n^{k}_{\varepsilon})$ and $y=W^{k}/g(n^{k})$. We know that $y_{\varepsilon}$, $y\in L^{\infty}({\mathbb{T}}^{d})\cap H^{1}({\mathbb{T}}^{d})$ and $y_{\varepsilon}\to y$ strongly in $L^{2}({\mathbb{T}}^{d})$. Thus, we can again apply Proposition 8 with $\Psi=\mathrm{Id}$ and infer that $W^{k}n^{k}\in H^{1}({\mathbb{T}}^{d})$ as well as

[TABLE]

for all $\phi_{0}\in H^{1}({\mathbb{T}}^{d})$.

Let $\phi_{1}$ be a smooth test function. We use the test function $W^{k}\phi_{1}$ in the weak formulation of (53):

[TABLE]

We pass to the limit $\varepsilon\to 0$ term by term. By (57),

[TABLE]

For the integral $I^{2}_{\varepsilon}$, we use the strong convergence (59) and the weak convergence of $(\nabla W^{k}_{\varepsilon})$ in $L^{2}({\mathbb{T}}^{d})$ to infer that

[TABLE]

The remaining integral $I^{3}_{\varepsilon}$ requires some work. As a preparation, using Proposition 8, we infer similarly to (58) that

[TABLE]

strongly in $L^{2}({\mathbb{T}}^{d})$. Let

[TABLE]

Then $(\xi_{\varepsilon})$ is bounded in $L^{2}(\Omega)$ and admits a weakly convergent subsequence, i.e. $\xi_{\varepsilon}\rightharpoonup\xi$ for some $\xi\in L^{2}({\mathbb{T}}^{d})$. Similarly as in the proof of Proposition 8, i.e. with an argument as in (35), we can find that, up to a subsequence, $(W^{k}_{\varepsilon}+\varepsilon)^{1/2}\xi_{\varepsilon}\rightharpoonup h$ weakly in $L^{2}({\mathbb{T}}^{d})$ implying $(W^{k})^{1/2}\xi=h$. In particular,

[TABLE]

Since $|f_{\varepsilon}(x)|\leq|h_{\varepsilon}(x)|$ for a.e. $x\in{\mathbb{T}}^{d}$ and all $\varepsilon>0$, we can apply Lemma 15 and obtain that, up to a subsequence, $f_{\varepsilon}$ converges strongly in $L^{2}(\Omega)$. Thus,

[TABLE]

Hence, passing to the limit $\varepsilon\to 0$ in (60), we infer that $(n^{k},W^{k})$ solves (9)-(10).

5.5. Energy estimate

We claim that the total energy $\int_{{\mathbb{T}}^{d}}W_{\rm tot}^{k}dx$ is nondecreasing in $k$. Let $(n^{k}_{\varepsilon},W^{k}_{\varepsilon})$ be a weak solution to (52)-(53). Then

[TABLE]

Taking the test functions $\phi_{1}=U$ in (53) and $\phi_{0}=n_{\varepsilon}^{k}$ in (52) and subtracting both equations, the above integral becomes

[TABLE]

Thus, with the lower semi-continuity of the norm, we have

[TABLE]

In view of mass conservation $\int_{{\mathbb{T}}^{d}}n^{k}dx=\int_{{\mathbb{T}}^{d}}n^{k-1}dx$, it follows by Jensen’s inequality that

[TABLE]

This shows the energy inequality in (42). Finally, assumption (43) gives $\int_{{\mathbb{T}}^{d}}W^{k}dx>0$ and consequently $W^{k}\not\equiv 0$.

5.6. An estimate for the variance

We claim that the total variance

[TABLE]

is nonincreasing in $k$. For the proof, we observe that, taking the test function $\phi_{1}=1$ in the weak formulation of (53) and performing the limit $\varepsilon\to 0$,

[TABLE]

Thus, by the energy estimate (61),

[TABLE]

We employ (61) again to find that

[TABLE]

In view of (63), the second bracket on the right-hand side is nonnegative, such that (62) leads to

[TABLE]

We take the test function $\phi_{1}=2\triangle t$ in (41):

[TABLE]

Combining the previous two inequalities, we arrive at

[TABLE]

Since the measure of ${\mathbb{T}}^{d}$ is one, we have

[TABLE]

Thus, taking into account mass conservation $\int_{{\mathbb{T}}^{d}}n^{k}dx=\int_{{\mathbb{T}}^{d}}n^{k-1}dx$ and $\int_{{\mathbb{T}}^{d}}W^{k-1}dx\leq\int_{{\mathbb{T}}^{d}}W^{k-2}dx$ (see (62)),

[TABLE]

and the claim follows after using the lower-semicontinuity of the $L^{2}$-norm.

5.7. Proof of Theorem 11

Let $d=1$. The second inequality in (42) implies that

[TABLE]

where $G=\max_{\delta\leq s\leq\{n^{k-1}\|_{L^{\infty}({\mathbb{T}})}\}}g(s)\geq\|g(n^{k})\|_{L^{\infty}({\mathbb{T}})}$. By the mean-value theorem, there exists $x_{0}\in{\mathbb{T}}$ such that

[TABLE]

Then, using Jensen’s inequality and the energy estimate in (42),

[TABLE]

By definition (39) of $V^{k}$, the right-hand side is positive if

[TABLE]

which is our assumption. Since $W^{k}\in H^{1}({\mathbb{T}})\hookrightarrow C^{0}([0,1])$, we conclude that $W^{k}>0$ in $[0,1]$. Then we can use $\phi_{1}=\phi/W^{k}$ as a test function in (41) and obtain the standard weak formulation of (41) for test functions $\phi\in H^{1}({\mathbb{T}})$. Furthermore, for $\phi\in H^{1}({\mathbb{T}})$,

[TABLE]

showing that $n^{k}\in H^{1}({\mathbb{T}})$ and finishing the proof.

6. Numerical simulations

We solve the one-dimensional equations (6) and (8) on the torus in conservative form, i.e. for the variables $n$ and $W_{\rm tot}$. The equations are discretized by the implicit Euler method and solved in a semi-implicit way:

[TABLE]

where $W^{k}_{\rm tot}=W^{k}-(U/2)(n^{k})^{2}$, $x\in{\mathbb{T}}=(0,1)$. The spatial derivatives are discretized by centered finite differences with constant space step $\triangle x>0$. For given $(W^{k-1},n^{k-1})$, the first equation (64) is solved for $n^{k}$. This solution is employed in the second equation (65) which is solved for $W^{k}_{\rm tot}$. Finally, we define $W^{k}=W^{k}_{\rm tot}+(U/2)(n^{k})^{2}$. We choose the parameters $U=10$, $\eta=1$, $\triangle t=10^{-5}$, and $\triangle x=10^{-2}$. The initial energy $W^{0}$ is constant and the initial density equals

[TABLE]

The time evolution of the particle density and energy is shown in Figure 1 for initial energy $W^{0}=1$. The variables converge to the constant steady state $(n^{\infty},W^{\infty})$ as $t\to\infty$, which is almost reached after time $t=0.1$. Since equations (64)-(65) are conservative, the total particle number $\int_{0}^{1}n(x,t)dx$ and the total energy $\int_{0}^{1}W_{\rm tot}(x,t)dx$ are constant in time. Consequently, the values for the steady state can be computed explicitly. We obtain for $W^{0}=1$,

[TABLE]

The energy stays positive for all times, so the high-temperature equations are strictly parabolic, and the convergence to the (constant) steady state is quite natural.

The situation is different in Figure 2, where the particle density converges to a nonconstant steady state $n^{\infty}$ (we have chosen $W^{0}=1/4$). This can be understood as follows. By contradiction, let both the particle density and energy be converging to a constant steady state. Then $n^{\infty}=1/2$ (see the above calculation) and

[TABLE]

However, this contradicts the fact that the energy $W$ is nonnegative which follows from the maximum principle. Therefore, it is plausible that either $n$ or $W$ cannot converge to a constant. If $n^{\infty}$ is not constant, $W^{\infty}\nabla n^{\infty}$ is constant only if $W^{\infty}=0$. Thus, it is reasonable that the energy converges to zero, while $n^{\infty}$ is not constant. One may say that there is not sufficient initial “reverted” energy to level the particle density.

Another difference between Figure 1 and Figure 2 is the time scale. For larger initial energies, the convergence to equilibrium is faster. In fact, Figure 3 shows that the decay of the $\ell^{2}$ norm of $n(t)-n^{\infty}$ and $W(t)-W^{\infty}$ is exponential. Here, we have chosen the initial particle density $n^{0}(x)=\frac{1}{4}$ for $0\leq x<\frac{1}{2}$, $n^{0}(x)=\frac{3}{4}$ for $\frac{1}{2}\leq x<1$, and the initial energy $W^{0}\in\{\frac{1}{4},\frac{1}{2},1\}$. For $W^{0}\in\{1,\frac{1}{2}\}$, we have $n^{\infty}=\frac{1}{2}$ and $W^{\infty}=\max\{0,W^{0}-U/32\}$. For $W^{0}=\frac{1}{4}$, it holds that $W^{\infty}=0$ and we have set $n^{\infty}(x)=n(x,2)$.

Finally, we compute the numerical convergence rates for different space and time step sizes $\triangle x$ and $\triangle t$, respectively. Since there is no explicit solution available, we choose as reference solution the solution to (64)-(65) with $\triangle x=1/1680$ (for the computation of the spatial $\ell^{2}_{x}$ error) and $\triangle t=1/5040$ (for the computation of the $\ell^{2}_{t}\ell^{2}_{x}$ error). Figure 4 shows that the temporal error is linear in $\triangle t$, and Figure 5 indicates that the spatial error is quadratic in $\triangle x$. These values are expected in view of our finite-difference discretization and they confirm the validity of the numerical scheme.

Appendix A Calculation of some integrals

We recall that $\varepsilon(p)=-2\varepsilon_{0}\sum_{k=1}^{d}\cos(2\pi p_{i})$. Then $u_{i}(p)=(\partial\varepsilon/\partial p_{i})(p)=4\pi\varepsilon_{0}\sin(2\pi p_{i})$, and we calculate

[TABLE]

since $\partial^{2}\varepsilon/\partial p_{i}\partial p_{j}=0$ for $i\neq j$. We compute the integral $\int_{{\mathbb{T}}^{d}}\varepsilon^{2}u_{i}^{2}dp$. First, let $d=1$. Then

[TABLE]

Furthermore, for $d>1$,

[TABLE]

We conclude that

[TABLE]