A reaction network approach to the convergence to equilibrium of quantum Boltzmann equations for Bose gases

When the temperature of a trapped Bose gas is below the Bose-Einstein transition temperature and above absolute zero, the gas is composed of two distinct components: the Bose-Einstein condensate and the cloud of thermal excitations. The dynamics of the excitations can be described by quantum Boltzmann models. We establish a connection between quantum Boltzmann models and chemical reaction networks. We prove that the discrete differential equations for these quantum Boltzmann models converge to an equilibrium point. Moreover, this point is unique for all initial conditions that satisfy the same conservation laws. In the proof, we then employ a toric dynamical system approach, similar to the one used to prove the global attractor conjecture, to study the convergence to equilibrium of quantum kinetic equations, derived in [49,50].


Introduction
Several years after the invention of the Boltzmann-Nordheim equation, which is the quantum version of the classical Boltzmann one, to describe the evolution of dilute quantum gases (cf. [40,51]), a renewal in the kinetic theory of bosons has started by the pioneering work of Kirkpatrick and Dorfman [34,35]. This work of Kirkpatrick and Dorfman was later extended by Zaremba, Nikuni and Griffin [54], in which the full coupling system of a quantum Boltzmann equation for the density function of the normal fluid/thermal cloud and a Gross-Pitaevskii equation for the wavefunction of the BEC has been introduced. In an independent work, the same model was derived by Pomeau, Brachet, Métens and Rica in [52]. We prefer to [26,41] for further discussions on the topic. In the models by Zaremba, Nikuni and Griffin and Pomeau, Brachet, Métens and Rica , there are two type of collisional processes.
• The 1 ↔ 2 interactions between the condensate and the excited atoms, described by the C 12 collision operator.
• The C 22 collision operator describes The 2 ↔ 2 interactions between the excited atoms themselves, described by the C 22 collision operator.
A third collisional process, previously missing, was proposed by Reichl and Gust [28,44]. This process takes into account 1↔3 type collisions between the excitations and is described by the collision operator C 31 . However, the derivation of the new collision operator C 31 was very complicated, since it involves the computations of around 40000 individual terms. As a result, a concise mathematical justification for the existence of the missing collision operator C 31 had been open for many years, and has been solved only until recently in [49]. The spatial homogeneous kinetic equation for the evolution of the density function f (t, p) of the thermal cloud, derived in [49][Section I], takes the form ∂ t f (p) = C 12 [f ](p) + C 22 [f ](p) + C 31 [f ](p), (1.1) in which the forms of C 12 , C 22 , C 31 are given explicitly below − f (p 1 )(f (p 2 ) + 1)(f (p 3 ) + 1)(f (p 4 ) + 1) , (1.4) in which n is the density of the condensate, t ∈ R + is the time variable, p ∈ (Z/L) d \{O} is the d-dimensional non-zero momentum variable, V is proportional to the volume of the periodic box − L 2 , L 2 d , m is the particle mass, ω is the Bogoliubov dispersion relation defined as (1. 5) and g is the interacting constant. We have normalized the Plank constant to be 1. In the above collision operators, the kernels are defined as follows and with u p and v p being defined as In the setting of [49], we could fix n as a constant, under the assumption that the thermal could fraction is quite small, in comparison to the condensate. Moreover, in the sum on the momenta p =0 , the origin is removed due to the fact that the condensate has been factored out in the Bogoliubov diagonalization (cf. [49,50]).

Remark 1.1
As it has been discussed in [49], the BEC is in a cubic box with periodic boundary conditions, the quantum Boltzmann equation is then in the discrete form. In order for the conservations of momentum and energy to be satisfied, the following system needs to have solutions on the lattice (1.10) At the first sign, the system does have solutions due to the complicated form of the Bogoliubov dispersion relation (1.5). However, it has been pointed out in [49,50] that when the temperature of the system is lower but closed to the Bose-Einstein condensation transition temperature, the Bogoliubov dispersion relation can be replaced by the Hatree-Fock energy (ω(p) ≈ c|p| 2 ). In this regime, the two collision operators C 12 and C 22 dominate the collisional processes. The contribution of third collision operator C 31 becomes non-trivial when both u p and v p are large, corresponding to significantly low temperatures. In this low temperature regime, the excitations are phonon-like and the Bogoliubov dispersion relation (1.5) can be replaced by the phonon dispersion relation (1.16). The replacement of (1.5) by (1.16) guarantees the existence of solutions to (1.10), and thus, the conservation laws are satisfied.
Simplified Quantum Boltzmann model of the thermal cloud. In our work, we try to provide a deeper understanding of the property of the system derived in [49] by studying a simplified version of it. If we denote then our simplified system for f 1 writes ]dp 2 dp 3 dp 4 , (1.13) ]dp 2 dp 3 , and ]dp 2 dp 3 dp 4 , The quantities K 22 p 1 ,p 2 ,p 3 ,p 4 , K 12 p 1 ,p 2 ,p 3 ≥ 0 are the collision kernels, which are radially symmetric, and symmetric with respect to the permutation of p 1 , p 2 , p 3 , and p 4 : Reaction networks and a toric dynamical system approach for the relaxation to equilibrium problem. The study of the relaxation of BECs to thermodynamic equilibrium has played very important role in the theory of Bose gases [43,29,28,25,54]. Our main tool is to convert these equations into chemical reaction systems and use an extension of the theory of toric dynamical systems (cf. [15]).
In general, there is great interest in understanding the qualitative behavior of deterministically modeled chemical reaction systems, including the existence of positive equilibria, stability properties of equilibria, and the non-extinction, or persistence, of species, which are the constituents of these systems [13,53,21,22,31,2,4,24,8,3,15]. Toric dynamical systems -originally called complex-balanced systems (cf. [15,32]) -are models used to describe an important class of chemical kinetics. The complex-balanced condition was first introduced by Boltzmann [11] for modeling collisions in kinetic gas theory. Based on this condition, it was shown by Horn and Jackson [32,30,20,27] that a complex-balanced system has a unique locally stable equilibrium within each linear invariant subspace. To underline the tight connection to the algebraic study of toric varieties, the name "toric dynamical system" was proposed in [15]. The most important problem in the theory of toric dynamical systems is the Global Attractor Conjecture, which says that the complex balanced equilibrium of a toric dynamical system is a globally attracting point within each linear invariant subspace. This global attractor question is strongly related to the convergence to equilibrium problem in the study of kinetic equations. A proof to the Global Attractor Conjecture for small dimensional systems has been supplied in [17], for strongly connected networks in [2], and a complete proof has been proposed in [14].
Our goal is to use the tools developed in [17,14] to prove the relaxation to equilibrium of Discrete Velocity Models of a model of (1.11), whose collision operator is C 12 . Similarly, we will prove the relaxation to equilibrium of another model of (1.11), whose collision operator is C 12 +C 22 , and modified quantum Boltzmann model of the thermal cloud (1.11), whose collision operator is C 12 +C 22 +C 13 . A related approach for the study of acoustic wave turbulence has been used in [48]. Let us also mention that some mathematical results of similar kinetic models have been obtained in [1,6,7,5,9,10,12,19,33,36,37,38,18,47,23,39,46,45].
The plan of our paper is the following: • In section 2, we show that the discrete version of a simplified version of (1.11), that contains only C 12 , could be rewritten as a chemical reaction network. By using an approach inspired by the theory of toric dynamical system, we prove in Theorem 2.1. that the solution of the discrete version of a simplified version of (1.11), that contains only C 12 , converges to the equilibrium exponentially in time.
• In section 3, we generalize Theorem 2.1 to collision operators of the forms C 13 and C 22 . We prove that the solutions of the discrete versions of these equations, associated with the collision operators C 13 and C 22 converge to equilibria exponentially in Theorems 3.1 and 3.2. In the case of C 22 , we consider a one-dimensional version of the model.
• In Theorem 4.1 of Section 4, we extend Theorem 3.2 to a simplified version of (1.11), that contains only C 12 + C 22 , and the modified quantum Boltzmann model of the thermal cloud (1.11), that contains only C 12 + C 22 + C 31 .
2 A reaction network approach for the case of C 12 2.1 The dynamical system associated to C 12 As mentioned in the introduction, the model derived from physics to describe the system that couples BEC-excitations at very low temperature is the discrete version of a simplified version of (1.11), that contains only C 12 , described below.
Let L R denote the lattice of integer points The discrete version of the simplified version of (1.11), that contains only C 12 , readṡ

Decoupling the quantum Boltzmann equation associated to C 12
Note that when p 1 = 0, K 12 p 1 ,p 2 ,p 3 is also 0, and therefore, we geṫ which says that f 0 (t) is a constant for all time t. Moreover, f p 1 does not depend on f 0 for all p 1 = 0. Therefore, without loss of generality, we can suppose that f 0 (0) = 0, which leads to f 0 (t) = 0 for all t.
Taking into account the fact E(p) = c|p|, note that if p 1 , p 2 , p 3 ∈ L R are different from 0 and p 3 = p 1 + p 2 and |p 3 | = |p 1 | + |p 2 | (like in the second sum of (2.1)), then p 1 , p 2 , p 3 must be collinear and on the same side of the origin. Therefore, we infer that there exists a vector P and k 1 , k 2 , k 3 > 0, k 1 , k 2 , k 3 ∈ Z such that Since L R is bounded, it follows that k 1 , k 2 , k 3 belong to a finite set of integer indices I = {1, . . . , I}. Arguing similarly for the first sum in (2.1), we deduce that (2.1) is equivalent with the following system for k 1 ∈ İ 3) Note that the system of equations (2.3) shows a decoupling of the system of equations (2.1) along a ray {kP 0 } with k > 0 (see Figure 1). As a consequence, it is sufficient to study the system of equations (2.3) for a fixed value of P 0 , instead of the system of equations (2.1).
If we denote f k 1 P 0 byf k 1 (with k 1 ∈ I) and K 12 k 1 P 0 ,k 2 P 0 ,k 3 P 0 by K 12 k 1 ,k 2 ,k 3 , we obtain the following new system for the ray {k 1 P 0 |k 1 > 0}: We denote this discrete version of C 12 by (2.7)

The chemical reaction network associated to C 12
For x ∈ R n >0 and α ∈ R n ≥0 , we denote by x α the monomial Π n i=1 x α i i .
where K is a positive parameter, called reaction rate constant. Then the mass-action dynamical system generated by this reaction iṡ

in which
x i is the concentration of the chemical species X i . For the case of a network that contains several reactions for 1 ≤ j ≤ m, its associated mass-action dynamical system is given bẏ In this section, we will show that the system (2.4) has the form (2.9) for a well-chosen set of reactions.
If y → y and y → y are reactions, we combine them together into a "reversible" reaction y ↔ y .
We will derive the system (2.4) from the network of chemical reactions of the form: for all k 1 , k 2 , k 3 in I such that k 2 + k 3 = k 1 . If we denote by F k the concentration of the species X k , we will show that, for appropriate choices of the reaction rate constants in (2.10) and (2.11), the differential equations satisfied by F k according the mass-action kinetics are exactly the same as (2.4).
In order to describe the connection between the mass-action system given by reactions of the form (2.10)-(2.11) and our system (2.4), we need to consider several cases.
and for the reversible reaction (2.12) the forward and backward rate constants are the same, i.e., we choose the reaction rate constants of the three reactions For example, consider the reversible reaction (2.12): in this reaction, X k 1 is created from Therefore, the rate of change of the species X k 1 due to this reaction is 2K 12 For the irreversible reaction (2.13), X k 1 is lost with the rate −2K 12 Therefore the rate of change of the species X k 1 due to this reaction is By exchanging the roles of X k 2 and X k 3 in (2.13), we obtain the rate −2K 12 . Therefore, the total rate of change of X k 1 due to the reactions in (2.12)-(2.13) is (2.14) Case 2: For 2k 2 = k 1 , k 1 , k 2 ∈ I, we consider We choose the reaction rate constant of 2X k 2 → X k 1 and the reaction rate constant of X k 1 → 2X k 2 to be K 12 k 1 ,k 2 ,k 3 . Also, we choose the reaction rate constant of X k 2 +X k 1 → 3X k 2 to be 2K 12 k 1 ,k 2 ,k 3 . Consider the first reaction (2.15): In this reaction, X k 1 is created from 2X k 2 with the rate K k 1 ,k 2 ,k 2 F 2 k 2 and X k 1 is decomposed into 2X k 2 with the rate −K 12 k 1 ,k 2 ,k 2 F k 1 . The rate of change of the species X k 1 is K 12 . For the second reaction (2.16): X k 1 is lost with the rate −2K 12 As a result, the rate of change of X k 1 due to the reactions (2.15)-(2.16) is Case 3: Next, for k 2 = k 3 + k 1 , k 1 = k 3 , k 1 , k 2 , k 3 ∈ I, let us look at the rate of change of X k 1 in For (2.18), the rate of change of X k 1 is 2K 12 For (2.19), the rate of change of X k 1 is 2K 12 By exchanging the roles of X 1 and X 3 , we obtain the Therefore, the rate of change of X k 1 due to reactions in (2.18)-(2.20) is Case 4: Now, for k 2 = 2k 1 , k 1 , k 2 ∈ I, let us look at the rate of change of X k 1 in For (2.22), the rate of change of X k 1 is 2K 12 Therefore, the rate of change of X k 1 due to the reactions (2.22)- From (2.14), (2.17), (2.21), (2.24), the total rate of change of X k 1 is which can be written aṡ which shows that the system of differential equations satisfied by the concentrations F k is exactly the same as the system of differential equations (2.4) satisfied by the densities f k .

A change of variables
In this section, we introduce a change of variables that will help us to investigate the dynamics of the system (2.26). .
The system (2.26) is converted intȯ , ∀k 1 ∈ I. (2.27) Suppose that G represents the column vector (G 1 , . . . , G I ) T . Let us also denote byX k , the vector  in which the only element that different from 0 is the k-th one.
Also, for k 1 = k 2 , we denote , Using these notations, the system (2.27) could be rewritten as: Equivalently, we can also writė where y ↔ y belongs to the set of reversible reactions with k 1 + k 2 = k 3 . Moreover, the solution f (t) of (2.1) converges to f * exponentially fast in the following sense: there exist positive constants C 1 , C 2 such that

Convergence to equilibrium
Proof By using the decoupling and the change of variables discussed in the previous sections, for each ray {kP 0 } k≥1 , we can reduce the study of f to F , which satisfies (2.26). From F , we can switch to study G, which is the solution of (2.29).
Step 1: The Lyapunov function. We recall that (2.27) could be rewritten under the forṁ (2.31) We define the function where G * k = 1 e kρ , for some ρ > 0, and we will show that L is a Lyapunov function for the system (2.27).
We have then H y,y = H y ,y for y and y as in (2.30). Moreover, we have since log is an increasing function. Also, note that the above inequality is strict unless for all reactions y ↔ y .
Since (G * ) y = (G * ) y for all reactions y ↔ y , this implies G * k 1 · G * k 2 = G * k 1 +k 2 for all k 1 and k 2 such that k 1 + k 2 ≤ I. As a consequence G * k = e −ρk , for some positive constant ρ. Moreover, (2.37) implies that at equilibrium (G) y = (G) y for all reactions y ↔ y , which leads to G k = e −ρ k , for some positive constant ρ . By the conservation relation By the monotonicity of the function ρ → e −ρk 1−e −ρk , we conclude that ρ = ρ , i.e., G * is the only equilibrium point that satisfies the same conservation relation as the initial condition. Now, we will prove that there exists exactly one critical point of the Lyapunov function L within each invariant set the projection of ∇L on the tangent space to the set S c is 0 if and only if there exists a constant such that A direct consequence of the above is the following system of identities Moreover, since G k and G * k satisfy the same conservation law then it follows that G = G * . This implies that G * is the only critical point of L on the invariant set S c .
Step 2: Differential inclusions and persistence. Now let us observe that (2.4) could be regarded as a K-variable mass-action system for the reversible network (2.30). For this we write and note that 1 + F k + F k is bounded below by 1 and above by 1 + 2C, where Therefore, the results of [14] about persistence of K-variable reversible mass-action systems can be applied and we conclude that the system is persistent. Alternatively, we can also use the Petri net argument of [4], to prove that the system is persistent, as follows. Note that F k is the density function of the species X k . It is straightforward that each siphon is {X 1 , X 2 , · · · , X I }, which contains the support of the P -semiflow (see [4] for the definition of siphons and P-semiflows) given by I k=1 kF k = constant.
As a result, the Petri net theory developed in [4] can be applied and it follows that the system is persistent.
Therefore, by using the existence of the globally defined strict Lyapunov function L, and the LaSalle invariance principle, it follows that all trajectories converge to the unique positive equilibrium G * that we discussed in Step 1.
Step 3: Exponential rate of convergence. Define and define Following [16], we compute the Jacobian of S at the equilibrium point G * , applied to an arbitrary vector δ = 0 that belongs to the span of the vectors y − y Jac(S(G * ))δ = y↔y K y↔y (G * ) y ((y − y ) * δ)H y,y (G * )(y − y ), (2.39) in which the inner product * is defined as Now, we compute the Jacobian of R at the equilibrium point G * , where the second equality is due to the fact that since G * is an equilibrium we have that S(G * ) = 0. Since is a diagonal matrix and A := Jac(S(G * )) is negative definite, then D 1/2 AD 1/2 is also negative definite with respect to this inner product. Since det(DA − λId) = det(D 1/2 AD 1/2 − λId), ∀λ ∈ R, it follows that D 1/2 AD 1/2 and DA have the same eigenvectors, so DA is negative definite. In other words, Jac(R(G * )) is negative definite. The exponential rate of convergence then follows from the fact that the Jacobian above is negative definite. This leads to the conclusion of the theorem.

Remark 2.1
The Lyapunov function (2.32) in the variable F reads and it is a strictly convex function.

Remark 2.2
If the intersection between the ray {kP 0 } k≥1 and L R contains a single point, then the solution f (t) of (2.1) has f P 0 ≡ 0, so f P 0 ≡ constant.
3 A reaction network approach for the case of C 13 and C 22 3.1 The dynamical system associated to C 13 As we discussed in the Introduction, we are also interested in the dynamics given by the discrete model of the collision operator C 13 , described in (1.15).
Let L R denote the lattice of integer points The discretized quantum Boltzmann equation for C 13 readṡ Similar to the C 12 case, when p = 0, K 13 p 1 ,p 2 ,p 3 ,p 4 = 0, and we obtaiṅ which means f 0 (t) is a constant for all time t, and we can assume f 0 (t) = 0 for all t.
Since in the first sum of (3.1), we consider (p 1 , p 2 , p 3 , p 4 ) satisfying we infer that there exists a vector P and k 1 , k 2 , k 3 , k 4 ≥ 0, k 1 , k 2 , k 3 , k 4 ∈ Z such that Using the same arguments as the case of C 12 , we can deduce that Equation (3.1) for C 13 is equivalent with the following family of decoupled systems for k 1 ∈ I = {1, 2, . . . , I} where P is the closest point to the origin among the lattice points on its ray: Denoting f kP by F k (with k ∈ I) and K 12 k 1 P,k 2 P,k 3 P,k 4 P by K 12 k 1 ,k 2 ,k 3 ,k 4 , we obtaiṅ In order to ensure that all the variables F k are coupled with each other, let us assume that I ≥ 4. We have the following conservation of energy for C 13 I k=1 kḞ k = 0, (3.5) or equivalently I k=1 kF k = const. (3.6) Similar to the case of C 12 , we define and then we have Note that, similar to the previous section, 0 < F k < ∞ and 0 < G k < 1.
The system (3.4) can be now writteṅ . (3.7) This system can also be rewritten aṡ (3.8) whereX k is, as mentioned earlier, the vector in which the only element that is 1 is the k-th one, and .
We can also writė where y ↔ y rang over the reversible reactions shown above.
Moreover, the solution f (t) of (3.1) converges to f * exponentially fast in the following sense: there exists positive constants C 1 , C 2 such that Proof The proof of Theorem 3.2 then follows exactly from the same Lyapunov function (2.32) and arguments as in Theorem 2.1.

The dynamical system associated to C 22
Let us consider a discretized version of the quantum Boltzmann model associated to the collision operator given by C 22 : Let L R denote the lattice of integer points The discretized quantum Boltzmann equation associated to C 22 reads ∀p 1 ∈ L Ṙ Similar to the C 12 case, when p = 0, K 22 p 1 ,p 2 ,p 3 ,p 4 = 0, and we obtaiṅ f 0 = 0, which means f 0 (t) is a constant for all time t. As a consequence, we can suppose that f 0 (0) = 0, which implies f 0 (t) = 0 for all t. In (3.9), the sums for C 22 are taken over (p 1 , p 2 , p 3 , p 4 ) satisfying In this case, unlike in the case of C 12 and C 13 , we cannot infer from (3.10) that there exists a vector P and k 1 , k 2 , k 3 , k 4 ≥ 0, k 1 , k 2 , k 3 , k 4 ∈ Z such that However, let us consider the following simplified version of (3.9) for C 22 (3.11) Recall that I = {1, · · · , I}. We also suppose that I ≥ 3. We have the following conservation of energy I k=1 kḞ k = 0, (3.12) or equivalently I k=1 kF k = const. (3.13) For C 22 , the following "conservation of mass" also holds I k=1Ḟ k = 0, (3.14) or equivalently I k=1 F k = const. (3.15) Similar to the case of C 12 , define and the system (3.9) can be now writteṅ . (3.16) This system can be rewritten aṡ in which the only element that is 1 is the k-th one, and .
We can proceed like in [32] to obtainρ =ρ andρ =ρ , yielding ρ 1 = ρ 1 and ρ 2 = ρ 2 . We can still use the Petri net argument of [4] or the result in [14], to prove that the system is persistent. For example, to use the method from [4], we note that we have two siphons {X 1 , X 2 , · · · , X I }, {X 2 , · · · , X I }. However, we also have the conservations of mass and energy k=1 F k = constant, k=1 kF k = constant, that leads to the P -semiflow k=2 (k − 1)F k = constant. Therefore, similar to the case of C 12 , it follows that the system is persistent, and we can use the same Lyapunov function as in the proof of Theorem 2.1 to obtain the desired convergence result.

Conclusion
In this work, we point out a connection between quantum Boltzmann models derived in [49] and chemical reaction network models. We prove that the discrete, simplified versions of some differential equations for these quantum Boltzmann models relax to an equilibrium point, by a toric dynamical system approach, similar to the one used in a recently proposed proof of the global attractor conjecture [14].