Energy-transport systems for optical lattices

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

Marcel Braukhoff Mathematisches Institut, Universität zu Köln, Weyertal 86-90, 50931 Köln, Germany  and  Ansgar Jüngel Institute for Analysis and Scientific Computing, Vienna University of Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria
July 6, 2019

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.

Key words and phrases:
Energy-transport models, optical lattice, degenerate equations, quadratic gradient, existence of weak solutions, finite differences.
2000 Mathematics Subject Classification:
35K59, 35K65, 35Q20, 82B40.
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 ,


where is the spatial variable, is the crystal momentum defined on the -dimensional torus with unit measure, and is the time. The velocity is defined by with the energy , is the potential, and is the collision operator. Compared to the semiconductor Boltzmann equation, there are two major differences.

First, the band energy is given by the periodic dispersion relation


The constant 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, [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 . This has important consequences regarding the integrability of the equilibrium distribution (see below).

Second, the potential is given by , where is the particle density and models the strength of the on-site interaction between spin-up and spin-down components [23]. In semiconductor physics, is the electric potential which is a given function or determined self-consistently from the (scaled) Poisson equation [15]. The definition 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

where is the relaxation time and 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

where are the Lagrange multipliers resulting from the mass and energy constraints. For , we obtain the Fermi-Dirac distribution, while for , equals the Maxwell-Boltzmann distribution. We may consider as a function of and write .

The variable can be interpreted as the negative inverse (absolute) temperature, while is related to the so-called chemical potential [15]. Since the band energy is bounded, the equilibrium is integrable even when , 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 and energy density solve the energy-transport equations


for , , where the particle and energy current densities are given by


and the diffusion coefficients depend nonlocally on and hence on ; see Proposition 1. The structure of system (3)-(4) is similar to the semiconductor case [17] but the diffusion coefficients are different. For , the Joule heating term contains the squared gradient , while in the semiconductor case it contains which is generally smoother than .

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

see Proposition 3. Here, the functions and are called the dual entropy variables, and the coefficients are defined in (20). In the dual entropy variables, the potential terms are eliminated, leading to the “symmetric” problem

where the matrix 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 and the entropy production . For instance, the determinant of the Jacobi matrix 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 ). 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 [23, Suppl.]. As can be interpreted as the temperature, the high-temperature limit corresponds to the limit . Expanding around up to zeroth order leads to the diffusion equation (see Section 3)


In the case , we obtain the logarithmic diffusion equation 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 . This motivates us to compute the next-order expansion. It turns out that at first order and with , is solving the (rescaled) energy-transport equations


where and is the “reverted” energy. The case corresponds to the maximal energy . Taking into account the periodic lattice structure, we solve (6)-(7) on the torus , together with the initial conditions , in .

The structure of the diffusion equation (6) is similar to (5), but the diffusion coefficient contains 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 . 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 (or ) is not conserved but the total energy . In fact, in terms of , the squared gradient term is eliminated,


Unfortunately, this formulation does not help for the analysis since the treatment of is delicate as lies in the dual space but is generally not an element of because of the degeneracy (we have only ).

The analysis of system (6)-(7) is very challenging since the first equation is degenerate in , 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

close to the constant equilibrium has been shown in [1]. The idea is that in such a situation, the temperature is strictly positive which removes the degeneracy in the term . 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

in a bounded domain, where . The global existence of weak solutions to the corresponding initial-boundary-value problem was proved. Again, the idea is a positivity bound for 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 . Since may vanish, we can expect a gradient estimate for only on . Although the quadratic gradient term also possesses 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


for , where and are given functions. We show the existence of a weak solution satisfying , and , ; see Theorem 9. In one space dimension and under a smallness assumption on the variance of and , the strict positivity of 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 since some estimates in the proof of Proposition 8 are not uniform in ; 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 is constant and sufficiently large, the solution converges to a constant steady state. However, if the constant 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 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 and energy .

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 in the diffusive scaling:


where is the Knudsen number [3], are the phase-space variables (space and crystal momentum), and is the time. We recall that the velocity equals , where the energy is given by (2). The potential is defined by . In the physical literature [23], the collision operator is given by the relaxation-time approximation

where the function is determined by maximizing the free energy (18) associated to (11) under the constraints


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

where and are the Lagrange multipliers and is a parameter which may take the values (Maxwell-Boltzmann statistics) or (Fermi-Dirac statistics). The relaxation time 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 be a (smooth) solution to the Boltzmann equation (11). We assume that the formal limits , , and exist. Then the particle and energy densities

are solutions to (3)-(4), and the diffusion coefficients are defined by

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.


To derive the balance equations, we multiply the Boltzmann equation (11) by and , respectively, and integrate over :


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 (which in fact defines ) in (13)-(14) and observe that the function is odd for any such that its integral over vanishes. This leads to

Passing to the formal limit gives the balance equations (3) with


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

and perform the formal limit ,


A straightforward computation shows that

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

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

In the following we write , where

Proposition 2 (Diffusion matrix).

The diffusion matrix is symmetric and positive definite.


The proof is similar to Proposition 3 in [17]. Let with , . Then

Since is symmetric in and , the symmetry of 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

The entropy density can be reformulated as


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

Proposition 3 (Entropy structure).

It holds that


where and are the so-called dual entropy variables and


Identity (18) implies that

and consequently,

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

which proves the identity in (19). Using the positive definiteness of , a computation shows that is positive definite too, and the inequality in (19) follows. ∎

2.3. Singularities and degeneracies in the energy-transport system

We denote by and the particle and energy densities depending on the dual entropy variable . We have the (implicit) formulation

Lemma 4.

Let , . Then


We differentiate

This gives after a rearrangement

In a similar way, we obtain

and with

the conclusion follows. ∎

Since can be positive and , the expression may vanish, so the determinant of may be not finite. Moreover, the numerator of the determinant may vanish, and the function may be not invertible. This is made more explicit in the following remark.

Remark 5 (Case ).

In the Maxwell-Boltzmann case, we can make the numerator of the determinant in Lemma 4 explicit. Indeed, it is clear that and . For the computation of , we observe first that


where, by symmetry,

Now, we have

Expanding the exponentials in and using (21)-(22), we find that

and this expression is independent of . For , we use the identity and integration by parts:

Summing over gives

We conclude that and consequently,

For certain values of or , this expression may vanish such that at these values. This shows that the relation between and 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

where , , are functions depending on , defined in Lemma 4, and on

The above formula shows that we lose the gradient estimate if . ∎

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 be given (to simplify). According to Example 6.8 in [15], we have


which is nonzero as long as . Furthermore, by Remark 8.12 in [15], it holds that , , and , and so

This expression is degenerate only at the boundary of the domain of definition (i.e. at ). 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 is interpreted as the negative inverse temperature, so high temperatures correspond to small values of . In this section, we perform a high-temperature expansion of (3)-(4), i.e., we expand around for small up to first order. Our ansatz is


3.1. Zeroth-order expansion

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

Therefore, up to order , we infer that

At high temperature, the relaxation time depends on the particle density in a nonlinear way, [20]. At low densities, i.e. , we obtain the logarithmic diffusion equation

where . 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),

Therefore, by (23),

and (16) and give, up to order ,