Alexandre Mouton

Thomas Rey

ON DETERMINISTIC NUMERICAL METHODS FOR THE QUANTUM BOLTZMANN-NORDHEIM EQUATION. I. SPECTRALLY ACCURATE APPROXIMATIONS, BOSE-EINSTEIN CONDENSATION, FERMI-DIRAC SATURATION

Keywords: Boltzmann-Nordheim equation, quantum, Bose-Einstein condensation, Fermi-Dirac saturation, spectral method, large time behavior 2010 Mathematics Subject Classification: 76P05, 65N35, 82C40

Spectral methods, thanks to their high accuracy and the possibility to use fast algorithms, represent an effective way to approximate the collisional kinetic equations of Boltzmann type, such as the Boltzmann-Nordheim equation. This equation, modeled on the seminal Boltzmann equation, describes using a statistical physics formalism the time evolution of a gas composed of bosons or fermions. Using the spectral-Galerkin algorithm introduced in [F. Filbet, J. Hu, and S. Jin, ESAIM: Math. Model. Numer. Anal., 2011], together with some novel parallelization techniques, we investigate some of the conjectured properties of the large time behavior of the solutions to this equation. In particular, we are able to observe numerically both Bose-Einstein condensation and Fermi-Dirac relaxation.

Introduction

The Boltzmann-Nordheim equation (BNE), also known as the quantum Boltzmann equation, models the time evolution of a phase-space density f = f (t, x, k), describing the probability to find at time t ≥ 0 a quantum particle localized at the infinitesimal position dx dk, where x ∈ Ω ⊂ R 3 and k ∈ R 3 . The quantity k designates the energy level of the particle and corresponds to the usual "velocity" variable in a kinetic Boltzmann-like equation. Its total momentum is then given by v = k, where is the reduced Planck constant, and we shall keep the v notation throughout the paper. We are interested in this work in efficient, deterministic numerical simulation of this equation.

This equation was first formulated by Uehling and Uhlenbeck in the seminal paper [START_REF] Uehling | Transport phenomena in Einstein-Bose and Fermi-Dirac gases[END_REF], starting from a classical Boltzmann equation with heuristic arguments. It can be used to model both bosons and fermions gases, possibly trapped by a confining potential. Due to its high dimensionality, the study of the full quantum Boltzmann equation is still widely open. The equation is indeed posed in the usual six-dimensional kinetic phase space. Moreover, the collision operator is defined by a sevenfold integral itself. This causes a lot of difficulty, both theoretically and numerically. For example, in the space homogeneous bosonic case, it is well known that the particle distribution function f can develop finite time blow-up (weak convergence towards Dirac deltas even if the kinetic energy is conserved), the so-called Bose-Einstein condensation.

The study of the main mathematical properties of the collision operator has been done in [START_REF] Dolbeault | Kinetic models and quantum effects: a modified boltzmann equation for fermi-dirac particles[END_REF][START_REF] Escobedo | On a quantum Boltzmann equation for a gas of photons[END_REF][START_REF] Lu | The Boltzmann equation for Bose-Einstein particles: regularity and condensation[END_REF][START_REF] Briant | On the Cauchy problem for the homogeneous Boltzmann-Nordheim equation for bosons: local existence, uniqueness and creation of moments[END_REF], allowing to understand almost completely the Cauchy problem (existence, stability, creation and propagation of moments, convergence towards the equilibrium) in the space homogeneous setting, without confining potential. The recent breakthrough [START_REF] Escobedo | Finite time blow-up for the boson Nordheim equation[END_REF] proved rigorously the Bose-Einstein condensation in the bosonic case, under very mild hypotheses on the moments of the initial condition. Nevertheless, the precise blowup behavior (localization of the explosion, rates, etc.) is still open. Finally, the theoretical study of the full space dependent problem is also mostly open, apart from some recent results concerning the anisotropic setting and Bose-Einstein condensates that can be found in [START_REF] Li | Global existence of solutions of the Boltzmann equation for Bose-Einstein particles with anisotropic initial data[END_REF].

The numerical study of this equation is also difficult for all these reasons. The first attempt to compute numerically the collision operator, in a simplified setting (taking advantage of some very specific symmetry properties) was done in the work [START_REF] Pareschi | Fast, conservative and entropic numerical methods for the Bosonic Boltzmann equation[END_REF]. Extensions of this idea were introduced in [START_REF] Hu | A fast spectral algorithm for the quantum Boltzmann collision operator[END_REF] by exploiting the convolution-like structure of the collision operator. This work is following along the lines of the so-called classical and fast spectral methods for the Boltzmann equation, that was investigated by many authors in the series of works [START_REF] Pareschi | A Fourier Spectral Method for Homogeneous Boltzmann Equations[END_REF][START_REF] Mouhot | Fast algorithms for computing the Boltzmann collision operator[END_REF][START_REF] Filbet | Solving the Boltzmann Equation in N log2 N[END_REF][START_REF] Filbet | Analysis of Spectral Methods for the Homogeneous Boltzmann Equation[END_REF][START_REF] Filbet | On steady-state preserving spectral methods for homogeneous Boltzmann equations[END_REF][START_REF] Pareschi | On the stability of equilibrium preserving spectral methods for the homogeneous Boltzmann equation[END_REF][START_REF] Hu | A new stability and convergence proof of the Fourier-Galerkin spectral method for the spatially homogeneous Boltzmann equation[END_REF][START_REF] Pareschi | Moment preserving Fourier-Galerkin spectral methods and application to the Boltzmann equation[END_REF]. The full extent of this convolution idea was used in [START_REF] Filbet | A numerical scheme for the quantum boltzmann equation efficient in the fluid regim[END_REF], allowing to write a fast (spectral) method able to compute the full collision operator with the "reasonably low" numerical cost of O(N 5 log(N )) operations, where N is the number of unknowns in each velocity dimension, and with spectral accuracy. Extension of this method to the full inhomogeneous case, and to the asymptotic preserving setting was then done in the series of works [START_REF] Hu | A fast algorithm for the energy space boson Boltzmann collision operator[END_REF][START_REF] Hu | An asymptotic-preserving scheme for the semiconductor Boltzmann equation with two-scale collisions: a splitting approach[END_REF][START_REF] Hu | An asymptotic-preserving scheme for the semiconductor Boltzmann equation toward the energy-transport limit[END_REF] for the fermionic (electrons) case only. Nevertheless, these works only deal with the simplified 1-dimensional (1d) in x, 2-dimensional in v setting and use an approximation of the collision operator. Even the slightly more realistic 1d in x, 3d in v case was never tackled up to the authors knowledge, not to mention the full 3d in x, 3d in v case. One should also note the earlier attempts from [START_REF] Jin | Discretization of the Multiscale Semiconductor Boltzmann Equation by Diffusive Relaxation Schemes[END_REF][START_REF] Klar | A Numerical Method for Kinetic Semiconductor Equations in the Drift-Diffusion Limit[END_REF], where a diffusive relaxation system was used, as an approximation of the full equation, and the Fokker-Planck like approach recently seen in [START_REF] Carrillo | Numerical study of Bose-Einstein condensation in the Kaniadakis-Quarati model for bosons[END_REF].

The paper is organized as follows: After a small introduction on the BNE equation and its main mathematical features, we detail in Section 1.2 its large time behavior, which depends on the type of considered particles: classical ones, Bosons, or Fermions, and on some critical parameters which insure condensation in the latter cases. We then present the fast algorithm from [START_REF] Filbet | A numerical scheme for the quantum boltzmann equation efficient in the fluid regim[END_REF] along with how we decided to implement it in Section 2. Some more information on this topic is also available in the Appendix. Section 3 is then devoted to the presentation of the so-called rescaling-velocity method, which allows to improve the accuracy of the fast spectral algorithm. Section 4 present some numerical results on the method, emphasizing first on the spectral accuracy of the algorithm in Section 4.1, and then on its long time behavior in Section 4.2, for the both 2D and 3D fermions and bosons cases.

1.1. The Boltzmann-Nordheim equation. We are then interested in the following kinetic equation in the d x × d v + 1-dimensional phase space, for d x ≤ d v ≤ 3:

(1.1)

∂ t f + τ v • ∇ x f = c Q α (f ) , f (0, x, v) = f 0 (x, v),
where

• x ∈ Ω x ⊂ R dx is the dimensionless position variable,
• v ∈ R dv is the dimensionless velocity variable,

• t ∈ R + is the dimensionless time variable,

• f : R + × Ω x × R dv → R is the distribution function of the particles,

• τ ∈ {0, 1} is a dimensionless parameter for indicating if the transport part is taken into account or not, • c ≥ 0 is a dimensionless scaling parameter for the collisions (usually taken as the inverse of the characteristic Knudsen number), • Q α (f ) : R + ×Ω x ×R dv → R is the collision operator, and α characterizes the type of particles, to be specified later.

In the present case, the collision operator Q α (f ) is defined as follows:

(1.2)

Q α (f )(t, x, v) := R dv ×S dv-1 B (|v -v * |, θ) f (t, x, v ) f (t, x, v * ) (1 -α f (t, x, v) -α f (t, x, v * )) -f (t, x, v) f (t, x, v * ) 1 -α f (t, x, v ) -α f (t, x, v * ) dv * dσ , with (v , v * , θ) ∈ R dv × R dv × [0, π] linked to (v, v * , σ) ∈ R dv × R dv × S dv-1 by the relations        v = v + v * 2 + |v -v * | 2 σ , v * = v + v * 2 - |v -v * | 2 σ ,
where

cos θ = v -v * |v -v * | • σ,
and the collision kernel B : R dv × [0, π] → R of variable hard sphere (VHS) type: [START_REF] Allemand | The incompressible Euler limit of the Boltzmann equation for a gas of fermions[END_REF] being given constants. Finally, α is defined as follows: α = + dv , for Fermi-Dirac particles, -dv , for Bose-Einstein particles.

B(r, θ) = Φ(r) b(cos θ) , Φ(r) = C Φ r γ , C Φ > 0, γ ∈ [0,
Remark 1 . Note that the case γ = 0 and d v = 2 represents the toy model of maxwellian pseudomolecules, whereas γ = 1 with d v = 3 models the physically relevant case of a hard sphere collision dynamics.

We also assume that the microscopic dynamics is subject to the so-called Grad's cut-off assumption, namely that the function b :

[-1, 1] → R + is such that b L ∞ (-1,1) < +∞.
Remark 2 . The Boltzmann-Nordheim collision operator, although being quadrilinear in its full form (1.2), is actually only trilinear in the two cases of interests α ∈ {0, 1} using cancellations in (1.3) .

Weak form and macroscopic properties. Note that, if we define

ω = v -v |v -v | , we have v = v -((v -v * ) • ω) ω , v * = v * + ((v -v * ) • ω) ω , v = v -((v -v * ) • ω) ω , v * = v * + ((v -v * ) • ω) ω ,
meaning that the transformation (v, v * ) → (v , v * ) is involutive. Thanks to this remark, one has the following weak form for the collision operator: for any test function Ψ such that the left hand side is defined,

(1.3) R dv Q α (f )(t, x, v) ψ(v) dv = R dv ×R dv ×S dv-1 B (|v -v * |, θ) f f * (1 -αf -αf * ) -f f * 1 -αf -αf * Ψ(v * ) + Ψ(v ) -Ψ(v) -Ψ(v * ) dv * dv dσ,
where we used the shorthand notations

f = f (t, x, v), f * = f (t, x, v * ), f = f (t, x, v ), f * = f (t, x, v * ).
In particular, one has using (1.1) that the collisior operator preserves the global mass, momentum and kinetic energy:

(1.4) R dv Q α (f )(t, x, v) ψ(v) dv = 0 , for ψ(v) = 1, ψ(v) = v, and ψ(v) = |v| 2 .
Finally, let us assume that the distribution function f satifies the following properties:

• f (t, x, v) > 0 for any t, x, v;

• For the Fermi-Dirac case (α > 0), f (t, x, v) < 1 α for any t, x, v.

Under these assumptions, we can define the local entropy associated to f as

H[f ](t, x) = - R dv f (t, x, v) ln (f (t, x, v)) + 1 -αf (t, x, v) α ln (1 -αf (t, x, v)) dv . Taking ψ(t, x, v) = ln f (t, x, v) 1 -α f (t, x, v) in (1.
3), we get the following expansion of the entropy dissipation functional:

D[f ](t, x) := R dv Q α (f )(t, x, v) ψ(t, x, v) dv = - 1 4 R dv R dv S dv-1 B(|v -v * |, θ) f f * (1 -αf -αf * ) -f f * 1 -αf -αf * log f f * (1 -αf -αf * ) f f * (1 -αf -αf * ) dσ dv * dv .
In addition, we have

R dv ∂ t f (t, x, v) ψ(t, x, v) dv = -∂ t H[f ](t, x, v) ,
and

R dv v • ∇ x f (t, x, v) ψ(t, x, v) dv = ∇ x • 1 α R dv v (ln (1 -α f (t, x, v))) dv .
We finally deduce that the global entropy

R dx H[f ](t, x) dx increases as t → +∞ since Ωx ∂ t H[f ](t, x) dx = - Ωx D[f ](t, x) dx ≥ 0 , ∀ t ≥ 0 . 1.2.
Convergence towards a steady state. In this section, we focus on the homogeneous equation associated to (1.1) that is characterized by τ = 0:

(1.5)

∂ t f = c Q α (f ) , f (0, v) = f 0 (v).
In such case, using (1.4) the momenta of order 0, 1, and 2 of f are constant:

(1.6)

ρ(t) = ρ 0 := R dv f 0 dv , u(t) = u 0 := 1 ρ 0 R dv v f 0 dv , e(t) = e 0 := 1 2ρ 0 R dv |v -u 0 | 2 f 0 dv .
It is well known that the distribution function f converges to a limit state f ∞ : R dv → R + that depends only on these d v + 2 quantities as t → +∞, see e.g. [START_REF] Briant | On the Cauchy problem for the homogeneous Boltzmann-Nordheim equation for bosons: local existence, uniqueness and creation of moments[END_REF][START_REF] Escobedo | Finite time blow-up for the boson Nordheim equation[END_REF][START_REF] Li | Global existence of solutions of the Boltzmann equation for Bose-Einstein particles with anisotropic initial data[END_REF]. To characterize this limit state, we proceed as follows:

The classical case α = 0. The limit state is the following maxwellian function [START_REF] Cercignani | The Mathematical Theory of Dilute Gases[END_REF] f

∞ (v) = M c (v) = ρ 0 (2πT ) dv/2 exp - |v -u 0 | 2 2T ,
where the temperature T is obtained from the internal energy e 0 thanks to the relation

T = 2 d v e 0 .
The quantum case α = 0. We introduce the so-called indicator function χ

(1.7) χ(v) = ρ 0 |D(u 0 , A)| 1 D(u 0 ,A) (v) , A = 2 e 0 (d v + 2) d v ,
and the following quantum maxwellian function [START_REF] Briant | On the Cauchy problem for the homogeneous Boltzmann-Nordheim equation for bosons: local existence, uniqueness and creation of moments[END_REF] (1.8)

M q (v) = 1 |α| 1 z -1 exp |v-u 0 | 2 2T + sgn(α)
, where z and T are linked to ρ 0 and e 0 as follows: since the mass and the internal energy are preserved we integrate M q (v) and |v -u 0 | 2 M q (v) to get

             ρ 0 = (2πT ) dv 2 |α| K dv 2 (z) , (1.9a) e 0 = d v T 2 K dv 2 +1 (z) K dv 2 (z) , (1.9b)
where K ν (z) is defined as

K ν (z) = 1 Γ(ν) +∞ 0 x ν-1 z -1 e x + sgn(α) dx = F ν (z) , if α > 0 (Fermi-Dirac), B ν (z) , if α < 0 (Bose-Einstein),
and Γ is the Gamma function

Γ(ν) = +∞ 0 x ν-1 e -x dx .
The quantities F ν and B ν are known as complete Fermi-Dirac and Bose-Einstein integrals respectively.

Remark 3 . Setting Fν (µ) = F ν (e µ ) for any µ ∈ R and ν > -1, one can show (see e.g [START_REF] Boersma | Asymptotic Expansion of a Class of Fermi-Dirac Integrals[END_REF]) that

       Fν (µ) = e µ + o(1) , as µ → -∞, Fν (µ) = µ ν Γ(ν + 1) 1 + π ν (ν -1) 6µ 2 + O 1 µ 4 , as µ → +∞.
One also has for ν > 1 that

B ν (z) = z + o(1) , as z → 0 + , B ν (z) = ζ(ν) + o(1) , as z → 1 -, where ζ(ν) = +∞ k = 1 1 k ν is Riemann's zeta function. Finally B 1 (z) → +∞ , as z → 1 -.
If we inject (1.9b) in (1.9a) in order to get rid of T , we obtain the following nonlinear equation in z:

(1.10) |α| ρ 0 ν 2πe 0 ν = K ν (z) ν+1 K ν+1 (z) ν , with ν = d v 2 .
Now, to distinguish between the fermionic and bosonic cases, let us define F ν and B ν as

F ν = F ν (z) ν+1 F ν+1 (z) ν , B ν (z) = B ν (z) ν+1 B ν+1 (z) ν .
We observe numerically that F ν and B ν are nondecreasing, continuous functions 1 on their respective definition domain (see Figure 1.1). In addition, thanks to the properties of

F ν and B ν in Remark 3 , one can show if d v = 2, 3 that          F dv 2 (z) = z + o(1) , as z → 0 + , F dv 2 (e µ ) = ( dv 2 + 1) dv 2 Γ( dv 2 + 1) 1 + O(µ -2 ) 1 + O(µ -2 ) , as µ → +∞,
1 Unfortunately, we were not able to prove this assertion. and that

B dv 2 (z) = z + o(1) , as z → 0 + .
More precisely, the following expansions can be obtained:

             B 1 (z) = log(1 -z) 2 Li 2 (1) = 6 log(1 -z) 2 π 2 → +∞ , as z → 1 -, B 3 2 (z) = Li 3 2 (z) 5 2
Li 5

2

(z)

3 2 → ζ( 3 2 ) 5 2 ζ( 5 2 ) 3 2 , as z → 1 -.
Finally, let us define the set I eq as (1.11)

I eq =                                    ]0, 2[ , for 2D Fermi-Dirac particles,   0, 5 3 10 π   , for 3D Fermi-Dirac particles, ]0, +∞[ , for 2D Bose-Einstein particles,   0, ζ( 3 2 ) 5 2 ζ( 5 2 ) 3 2
  , for 3D Bose-Einstein particles.

Then, we have shown that there exists a unique solution of (1.10) if and only if there holds

|α| ρ 0 d v 4πe 0 dv 2 ∈ I eq .
In such case, (1.10) can be solved numerically to get z, then T thanks to (1.9) and identify the limit state f ∞ as M q . On the other hand, if

|α| ρ 0 d v 4πe 0 dv 2 / ∈ I eq ,
we handle a quantum degenerative case and the limit state f ∞ cannot be identified as M q . Recent works (see e.g. [START_REF] Briant | On the Cauchy problem for the homogeneous Boltzmann-Nordheim equation for bosons: local existence, uniqueness and creation of moments[END_REF][START_REF] Lu | The Boltzmann equation for Bose-Einstein particles: regularity and condensation[END_REF][START_REF] Escobedo | Finite time blow-up for the boson Nordheim equation[END_REF]) indicates that, for the 3D degenerate Bose-Einstein case which describes a quantum condensation, the limit state f ∞ is identified as

(1.12) M q (v) = m 0 δ 0 (v -u 0 ) + 1 |α| 1 exp |v-u 0 | 2 2T -1
, where δ 0 is the usual Dirac function, the critical mass m 0 and the temperature T are linked to the density ρ 0 and the internal energy e 0 as (1.13)

T = 2e 0 ζ(3/2) 3ζ(5/2) ,
and

(1.14) m 0 = ρ 0 - (2πT ) 3/2 |α| ζ(3/2) = ρ 0 - 1 |α| 4πe 0 3 3/2 ζ(3/2) 5/2 ζ(5/2) 3/2 .
We can remark here that

|α| ρ 0 d v 4πe 0 dv 2 / ∈ I eq ⇐⇒ m 0 > 0 .
Concerning Fermi-Dirac degenerative cases, we have to deal with the so-called "maximum assumption" max f ≤ |α| -1 in order to define the entropy. It was shown in [START_REF] Escobedo | Entropy maximisation problem for quantum relativistic particles[END_REF][START_REF] Allemand | The incompressible Euler limit of the Boltzmann equation for a gas of fermions[END_REF] that, in the critical case where ρ 0 and e 0 satisfy (1.15)

ρ 0 |D(0, A)| = 1 |α| , A = 2 e 0 (d v + 2) d v ,
the distribution function f converges to the limit state χ defined in (1.7). If the maximum assumption is broken, the entropy cannot be defined anymore as H[f ] so nothing clear can be said about the limit state in this context.

Practical computations. Let us classify the limit states according to the initial parameters: Mass ρ 0 and internal energy e 0 given. We want to compute the fugacity z and the temperature T . For this purpose, we test the assertion

|α| ρ 0 d v 4πe 0 dv 2 ∈ I eq .
• If it holds, we solve numerically (1.10) with Brent's method [START_REF] Brent | Algorithms for Minimization without Derivatives[END_REF] to get an approximation of z. When Fermi-Dirac particles are considered, we use the scale z = e µ and run the iterative method to find µ. This allows to get large values of z when |α| ρ 0 dv 4πe 0 dv2 is close to the upper bound of I eq . The limit distribution is then set 2 as f ∞ = M q . • If it does not hold, -when 3D Bose-Einstein particles are considered, we use (1.13)-(1.14) to set T, m 0 and we set z = 1 and the limit state is defined again as f ∞ = M q . -For 2D and 3D Fermi-Dirac cases, we test the equality (1.15).

* If it holds, the limit state is identified as f ∞ = χ and the relaxation is entropic. * If it does not hold, we cannot say anything about the behaviour of f (t, •) as t → +∞ since it is non-entropic.

Mass ρ and temperature T given. We want to compute the fugacity z and the internal energy e 0 .

• If the 3D Bose-Einstein case is considered, we test the inequality

|α| ρ 0 ≤ (2πT ) 3 2 ζ( 3 2
).

-If it is true, then we can solve numerically (1.10) with Brent's method to get an approximation of z. We use it to get e 0 from (1.9b) and set the limit distribution as f ∞ = M q . -If it is false, then we use (1.13)-(1.14) to set T, m 0 and we set z = 1 when Bose-Einstein particles are considered and the limit state is defined as f ∞ = M q . • If 2D Bose-Einstein or Fermi-Dirac cases are considered, we can solve numerically (1.10) with Brent's method to get an approximation of z, then set the limit distribution as f ∞ = M q and compute e 0 thanks to (1.9b).

Fast deterministic approximations of the Boltzmann-Nordheim collision operator

Let us now introduce the numerical method that we shall study and use in this paper. This section is strongly inspired from [START_REF] Filbet | A numerical scheme for the quantum boltzmann equation efficient in the fluid regim[END_REF]. For the sake of simplicity, all the functions described in this section will only depend on the velocity variable, denoted by ξ.

2.1. The Carleman representation. Let us decompose the BNE operator (1.2) using the following bilinear and trilinear operators (see also Remark 2 ):

Q + c (F, G)(ξ) = R dv S dv-1 B (|ξ -ξ * |, θ) F G * dσ dξ * , Q - c (F, G)(ξ) = F (ξ) R dv S dv-1 B (|ξ -ξ * |, θ) G * dσ dξ * , Q 1,q (F, G, H)(ξ) = R dv S dv-1 B (|ξ -ξ * |, θ) F G * H * dσ dξ * , Q 2,q (F, G, H)(ξ) = F (ξ) R dv S dv-1 B (|ξ -ξ * |, θ) G H * dσ dξ * , Q 3,q (F, G, H)(ξ) = F (ξ) R dv S dv-1 B (|ξ -ξ * |, θ) G * H dσ dξ * , Q 4,q (F, G, H)(ξ) = F (ξ) R dv S dv-1 B (|ξ -ξ * |, θ) G * H * dσ dξ * .
(2.1)

Hence the complete collision operator Q(F, F ) writes as

Q α (F ) = Q + c (F, F ) -Q - c (F, F ) -α [Q 1,q (F, F, F ) + Q 2,q (F, F, F ) -Q 3,q (F, F, F ) -Q 4,q (F, F, F )] .
Remark 4 . One can notice in (2.1) that the only term that consists on a full integral on the three distributions F , G and H is Q 1,q . This term will then logically be the most computationally expensive.

By completing the square, one has the following useful lemma from [START_REF] Mouhot | Fast algorithms for computing the Boltzmann collision operator[END_REF]:

Lemma 2.1. For any function F : R dv → R n , we have

1 2 S dv-1 F (|w| σ -w) dσ = 1 |w| dv-2 R dv δ 2x • w + |x| 2 F (x) dx .
Thanks to this lemma, one has the following Carleman-like representation for the classical collision operator:

Q + c (F, G)(ξ) -Q - c (F, G)(ξ) = R 2dv B(z, y) δ(z • y) [F (ξ + z) G(ξ + y) -F (ξ) G(ξ + z + y)] dz dy ,
with B defined as

B(z, y) = 2 dv-1 C Φ (|z| 2 + |y| 2 ) γ-dv+2 2 b   |y| |y| 2 + |z| 2 , |z| |y| 2 + |z| 2   .
Similar computations lead to the following Carleman-like representations for the quantum, trilinear collision operators (with the same effect on computational complexity as in Remark 4 , namely that Q 1,q will be the more computationally intensive operator)

Q 1,q (F, G, H)(ξ) = R 2dv B(z, y) δ(z • y) F (ξ + z) G(ξ + y) H(ξ + z + y) dz dy , Q 2,q (F, G, H)(ξ) = F (ξ) R 2dv B(z, y) δ(z • y) G(ξ + z) H(ξ + y) dz dy , Q 3,q (F, G, H)(ξ) = F (ξ) R 2dv B(z, y) δ(z • y) G(ξ + z + y) H(ξ + z) dz dy , Q 4,q (F, G, G)(ξ) = F (ξ) R 2dv B(z, y) δ(z • y) G(y) H(ξ + y) dz dy .
2.2. Truncation and periodization of the BNE operator. In order to develop a fast spectral approximation for the Boltzmann-Nordheim operator, one has first to truncate it, in order to make it deal with function with finite support. We shall follow along the lines of [START_REF] Filbet | A numerical scheme for the quantum boltzmann equation efficient in the fluid regim[END_REF].

Theorem 2.1. Let the function G be compactly supported in the d v -dimensional ball B(0, S). Then, the support of

ξ → Q α (G)(ξ) is included in B(0, √ 2S).
Proof. Considering z, y ∈ R dv such that z • y = 0, we have

|ξ| 2 + |ξ + z + y| 2 = |ξ + z| 2 + |ξ + y| 2 .
Hence, if ξ, ξ + z, ξ + y, ξ + z + y ∈ B(0, S), it is possible to have

G(ξ+z) G(ξ+y) [1 -α (G(ξ) + G(ξ + z + y))]-G(ξ) G(ξ+z+y) [1 -α (G(ξ + z) + G(ξ + y))] = 0 .
In this situation, we have

|ξ| 2 ≤ |ξ| 2 + |ξ + z + y| 2 ≤ 2S 2 , so we can conclude that the support of Q α (G) is included in the ball B 0, √ 2S if the support of G is in B(0, S) in ξ.
We now define a truncated collision operator Q R α with R > 0 as follows: given a function G :

R dv → R with compact support in B(0, S), we define G : [-π, π] dv → R such that G(ξ) = G(ξ) , if ξ ∈ B(0, S), 0 , else.
At this point, we have several candidates for Q R α , each one of them being written under the following generic way

Q R α (G)(ξ) = R 2dv B(x, y) δ(x • y) G(ξ + x) G(ξ + y) 1 -α G(ξ) + G(ξ + x + y) -G(ξ) G(ξ + x + y) 1 -α G(ξ + x) + G(ξ + y) I(x, y, ξ) dx dy ,
where I is defined as follows:

• the classical spectral method from [29]

I(x, y, ξ) = 1 B(0,S) (ξ + x + y) 1 B(0,2S) (x -y) 1 B(0, √ 2S) (ξ) , • the fast spectral method from [26] I(x, y, ξ) = 1 B(0,2R) (x -y) 1 B(0, √ 2S) (ξ) , R ≥ S ,
• the fast quantum Boltzmann spectral method from [START_REF] Filbet | A numerical scheme for the quantum boltzmann equation efficient in the fluid regim[END_REF]:

(2.2) I(x, y, ξ) = 1 B(0,R) (x) 1 B(0,R) (y) 1 B(0, √ 2S) (ξ) , R ≥ S ,
In order to simplify the computations of the following section, we consider the definition (2.2) for I.

It gives the following expression of

Q R α (G): Q R α (G)(ξ) = B(0,R) 2 B(x, y) δ(x • y) G(ξ + x) G(ξ + y) 1 -α G(ξ) + G(ξ + x + y) -G(ξ) G(ξ + x + y) 1 -α G(ξ + x) + G(ξ + y) dx dy .
We now define the periodized function G as

G(ξ + 2kπ) = G(ξ) , ∀ ξ ∈ [-π, π] dv , ∀ k ∈ Z dv .
Choosing R > 0 and S > 0 such that R ≥ S and S is small enough allows to write

Q R α ( G)(ξ) = Q α ( G)(ξ) ∀ξ ∈ [-π, π] dv .
More precisely, since G is supported in B(0, S), G is supported on the balls B(2kπ, S) for any k ∈ Z dv . Consequently, if a particle is equipped with a velocity ξ ∈ B(2kπ, S), it may interact with any particle equipped with a velocity ξ * ∈ B 2kπ, 2R + √ 2S . Such candidate particle should not be in any B(2lπ, S) for all l ∈ Z dv for avoiding artificial interactions. In order to guarantee this constraint, we must impose 2π -S ≥ 2R + √ 2S, namely

(2.3) S ≤ 2π 2R + (1 + √ 2) .
This condition is generally referred as the anti-aliasing condition on S and R.

Fast spectral approximation of the BNE operator.

For now on, us assume that Supp ξ G ⊂ B(0, S) ⊂ [-π, π] dv and consider a truncation radius R > 0 satisfying R ≥ 2S and (2.3). Hence

Q α (G)(v) = Q R α (G)(v) . Given a function f ∈ L 2 per [-π, π] dv )
, we define its Fourier series representation as

(2.4) f (v) = k∈Z dv f k e ik•v , f k = 1 (2π) dv [-π,π] dv f (v)e -ik•v dv,
where we used the multi-index notation k = (k 1 , . . . , k dv ) to denote the d v -dimensional sums over the indexes k j , j = 1, . . . , d v .

We introduce the truncated terms

Q R 1,c , Q R 2,c , Q R 1,q , Q R 2,q , Q R 3,q , Q R 4,q as Q R 1,c (F, G)(ξ) = B(0,R) 2 B(z, y) δ(z • y) F (ξ + z) G(ξ + y) dz dy , Q R 2,c (F, G)(ξ) = F (ξ) B(0,R) 2 B(z, y) δ(z • y) G(ξ + z + y) dz dy , Q R 1,q (F, G, H)(ξ) = B(0,R) 2 B(z, y) δ(z • y) F (ξ + z) G(ξ + y) H(ξ + z + y) dz dy , Q R 2,q (F, G, H)(ξ) = F (ξ) B(0,R) 2 B(z, y) δ(z • y) G(ξ + z) H(ξ + y) dz dy , Q R 3,q (F, G, H)(ξ) = F (ξ) B(0,R) 2 B(z, y) δ(z • y) G(ξ + z + y) H(ξ + z) dz dy , Q R 4,q (F, G, H)(ξ) = F (ξ) B(0,R) 2 B(z, y) δ(z • y) G(y) H(ξ + y) dz dy .
Using orthogonality properties of orthogonal polynomials, one has the following expansion of the quantum kernel modes. Lemma 2.2. For any N ∈ N * and any multi-index n, one has the following representation of the Fourier modes of the terms of the BNE collision operator:

Q R 1,c (F, G) n = k,l ∈ I N k+l = n β R (k, l) F k G l , Q R 2,c (F, G) n = k,l ∈ I N k+l = n β R (l, l) F k G l , Q R 1,q (F, G, H) n = k,l,m ∈ I N k+l+m = n β R (k + m, l + m) F k G l H m , Q R 2,q (F, G, H) n = k,l,m ∈ I N k+l+m = n β R (k, l) F k G l H m , Q R 3,q (F, G, H) n = k,l,m ∈ I N k+l+m = n β R (k + l, l) F k G l H m , Q R 4,q (F, G, H) n = k,l,m ∈ I N k+l+m = n β R (l, k + l) F k G l H m ,
with the so-called kernel modes β R (k, l) being defined as

β R (k, l) = B(0,R) 2 B(z, y) δ(z • y) e i (z•k+y•l) dz dy ,
and the set I N given by

I N = dv d = 1 - N 2 , . . . , N 2 -1 ⊂ Z dv .
Such an expansion of the kernel mode yields an asymptotic complexity for computing the BNE operator of O(N 2dv ) for the classical terms and O(N 3dv ) for the term Q R 1,q . We shall now expand as in [START_REF] Mouhot | Fast algorithms for computing the Boltzmann collision operator[END_REF] the collision kernel in order to introduce a convolution sum, allowing to improve the overall computational complexity by almost d v orders of magnitude.

β R (k, l) = S dv-1 S dv-1 δ(e • e ) R 0 a (|ρ|) ρ dv-2 e i ρ k•e dρ R 0 b (|ρ|) (ρ ) dv-2 e i ρ l•e dρ de de . 2.4.1. 2D case. In such case, defining e θ = (cos θ, sin θ), the expression of β R (k, l) is simplified into β R (k, l) = 2π 0 2π 0 δ(e θ • e θ ) R 0 a (|ρ|) e i ρ k•e θ dρ R 0 b (|ρ|) e i ρ l•e θ dρ dθ dθ .
Taking advantage of symmetries of the kernel, one then has

β R (k, l) = π 0 π 0 δ(e θ • e θ ) R -R a (|ρ|) e i ρ k•e θ dρ R -R b (|ρ|) e i ρ l•e θ dρ dθ dθ Since e θ • e θ = 0 with θ, θ ∈ [0, π] is equivalent to θ ≡ θ + π 2 (modπ), we finally get β R (k, l) = π 0 φ 2 R,a (k • e θ ) φ 2 R,b (k • e θ+ π 2 (modπ) ) dθ , with φ 2 R,a and φ 2 R,b defined as φ 2 R,a (s) = R -R a (|ρ|) e iρs dρ , φ 2 R,b (s) = R -R b (|ρ|) e iρs dρ . Since φ 2 R,a and φ 2 R,b are even, we deduce that θ → φ 2 R,a (k • e θ ) φ 2 R,b (k • e θ+ π 2 (modπ)
) is π-periodic. Hence, the rectangular integration formula is of infinite order and is optimal [START_REF] Shen | Spectral methods: algorithms, analysis and applications[END_REF], so one can write

(2.5) β R (k, l) ≈ π M M -1 p = 0 α p (k) α p (l) , with (2.6) α p (k) = φ 2 R,a (k • e θp ) , α p (l) = φ 2 R,b (l • e θp+ π 2 (modπ) ) , θ p = p π M ,
and a spectrally small error.

Remark 5 . In the specific case where a = b, the function

θ → φ 2 R,a (k • e θ ) φ 2 R,b (k • e θ+ π 2 (modπ) ) is π 2 -periodic. Consequently, we have (2.7) β R (k, l) = 2 π 2 0 φ 2 R,a (k • e θ ) φ 2 R,a (k • e θ+ π 2 (modπ) ) dθ ,
and the rectangular integration formula gives

(2.8) β R (k, l) ≈ π M M -1 p = 0 α p (k) α p (l) , with (2.9) α p (k) = φ 2 R,a (k • e θp ) , α p (l) = φ 2 R,a (l • e θp+ π 2 (modπ) ) , θ p = p π 2M ,
and a spectrally error.

2.4.2. 3D case. In the 3d case, tedious (but similar to the previous section) computations yield

β R (k, l) = 1 4 [0,π] 4 sin θ sin θ δ(e θ,ϕ • e θ ,ϕ ) φ 3 R,a (k • e θ,ϕ ) φ 3 R,b (l • e θ ,ϕ ) + δ(e θ,ϕ+π • e θ ,ϕ ) φ 3 R,a (k • e θ,ϕ+π ) φ 3 R,b (l • e θ ,ϕ ) + δ(e θ,ϕ • e θ ,ϕ +π ) φ 3 R,a (k • e θ,ϕ ) φ 3 R,b (l • e θ ,ϕ +π ) + δ(e θ,ϕ+π • e θ ,ϕ +π ) φ 3 R,a (k • e θ,ϕ+π ) φ 3 R,b (l • e θ ,ϕ +π ) dϕ dϕ dθ dθ .
with φ 3 R,a and φ 3 R,b defined as

φ 3 R,a (s) = R -R a (|ρ|) |ρ|e iρ s dρ , φ 3 R,b (s) = R -R b (|ρ|) |ρ|e iρ s dρ .
We can remark that both φ 3 R,a and φ 3 R,b functions are even. Hence, one has

β R (k, l) = 1 4 [0,π] 4 sin θ sin θ δ(e θ,ϕ • e θ ,ϕ ) φ 3 R,a (k • e θ,ϕ ) φ 3 R,b (l • e θ ,ϕ ) + δ(-e π-θ,ϕ • e θ ,ϕ ) φ 3 R,a (-k • e π-θ,ϕ ) φ 3 R,b (l • e θ ,ϕ ) + δ(-e θ,ϕ • e π-θ ,ϕ ) φ 3 R,a (k • e θ,ϕ ) φ 3 R,b (-l • e π-θ ,ϕ ) + δ(e π-θ,ϕ • e π-θ ,ϕ ) φ 3 R,a (-k • e π-θ,ϕ ) φ 3 R,b (-l • e π-θ ,ϕ
) dϕ dϕ dθ dθ , so we finally get

β R (k, l) = [0,π] 4 sin θ sin θ δ(e θ,ϕ • e θ ,ϕ ) φ 3 R,a (k • e θ,ϕ ) × φ 3 R,b (l • e θ ,ϕ ) dϕ dϕ dθ dθ .
For any e ∈ S 2 , we define l → Π e ⊥ (l) as

Π e ⊥ (l) = l -(l • e) e .
Hence, β R (k, l) is reformulated as follows:

β R (k, l) = [0,π] 2 sin θ φ 3 R,a (k • e θ,ϕ )   (θ ,ϕ ) ∈ [0,π] 2 e θ,ϕ •e θ ,ϕ =0 sin θ φ 3 R,b (l • e θ ,ϕ ) dϕ dθ   dϕ dθ .
Denoting with S 2 + the half-sphere, we can write

(θ ,ϕ ) ∈ [0,π] 2 e θ,ϕ •e θ ,ϕ =0 sin θ φ 3 R,b (l • e θ ,ϕ ) dϕ dθ = S 2 + ∩e ⊥ θ,ϕ φ 3 R,b (l • e ) de .
We define e θ l ,

ϕ l = Π e ⊥ θ,ϕ (l) Π e ⊥ θ,ϕ (l)
. Hence the 2D plan e ⊥ θ,ϕ is provided with the direct orthonormal basis {e θ l ,ϕ l , e θ,ϕ × e θ l ,ϕ l }. In addition, the intersection S 2 + ∩ e ⊥ θ,ϕ is the unit half-circle S 1 + and

l • e = Π e ⊥ θ,ϕ (l) e θ l ,ϕ l • e ,
where e θ l ,ϕ l •e is no more than the cosinus of the angle between e θ l ,ϕ l and e ∈ S 2 + ∩e ⊥ θ,ϕ . Consequently, we obtain

(θ ,ϕ ) ∈ [0,π] 2 e θ,ϕ •e θ ,ϕ =0 sin θ φ 3 R,b (l • e θ ,ϕ ) dϕ dθ = π 0 sin θ φ 3 R,b Π e ⊥ θ,ϕ (l) cos θ dθ =: ψ 3 R,b (Π e ⊥ θ,ϕ (l)) ,
and

β R (k, l) = [0,π] 2 sin θ φ 3 R,a (k • e θ,ϕ )ψ 3 R,b (Π e ⊥ θ,ϕ (l)) dϕ dθ .
We can notice that the function κ(θ, ϕ)

= sin θ φ 3 R,a (k • e θ,ϕ )ψ 3 R,b (Π e ⊥ θ,ϕ (l)
) is π-periodic in θ and in ϕ. Indeed, reminding that e θ,ϕ+π = e π-θ,ϕ firstly gives

κ(θ, ϕ + π) = κ(π -θ, ϕ) ,
so the π-periodicity in ϕ will be proved along with the π-periodicity in θ. This last property can be easily proved since sin is π-antiperiodic, e θ+π,ϕ = -e θ,ϕ , φ 3 R,a is even, and

Π e ⊥ θ+π,ϕ (l) = Π e ⊥ θ,ϕ (l) ,
which makes the basis {e θ l ,ϕ l , e θ+π,ϕ × e θ l ,ϕ l } indirect. Consequently, we have

ψ 3 R,b (Π e ⊥ θ+π,ϕ (l)) = -ψ 3 R,b (Π e ⊥ θ,ϕ (l)) .
Since κ is π-periodic in θ and ϕ, the rectangular integration rule on [0, π] 2 is of infinite order and optimal, so we can write

(2.10) β R (k, l) ≈ π 2 M 1 M 2 M 1 -1 p = 0 M 2 -1 q = 0 sin θ p φ 3 R,a (k • e θp,ϕq )ψ 3 R,b (Π e ⊥ θp,ϕq (l)) ,
with e θp,ϕq = (sin θ p cos ϕ q , sin θ p sin ϕ q , cos θ p ),

θ p = p π M 1 , ϕ q = π q M 2
, and

(2.11)

φ 3 R,a (s) = R -R a (|ρ|) |ρ|e iρ s dρ , φ 3 R,b (s) = R -R b (|ρ|) |ρ|e iρ s dρ , ψ 3 R,b (Π e ⊥ θ,ϕ (l)) = π 0 sin θ φ 3 R,b Π e ⊥ θ,ϕ
(l) cos θ dθ .

Fast computation of the collision modes.

It is possible to gather the formula (2.5) and (2.10) under the common expression

(2.12) β R (k, l) = C P -1 p = 0 α R,p (k) α R,p (l) ,
where C, P , α p , α p are defined as follows:

2D case. Let C = π M
, P = M and, for any p = 0, . . . , M -1, (2.13)

               α R,p (k) = φ 2 R,a (k • e θp ) , α R,p (l) = φ 2 R,b (l • e θp+ π 2 (modπ) ) , θ p = p π M , e θp = (cos θ p , sin θ p ) , φ 2 R,a (s) = R -R a (|ρ|) e iρ s dρ , φ 2 R,b (s) = R -R b (|ρ|) e iρ s dρ ,
In the specific case of 2D Maxwellian molecules (γ = 0, b ≡ 1) so we get

B(z, y) = 2C Φ a(|z|) = 2C Φ , b(|y|) = 1 .
This gives

β R (k, l) = π M P -1 p = 0 α R,p (k) α R,p (l) ,
with α R,p (k) and α R,p (l) computed as follows:

α R,p (k) = 2C Φ φ 2 R (k • e θp ) , α R,p (l) = φ 2 R (l • e θp+ π 2 (modπ) ) , φ 2 R (s) = 2R Sinc(Rs) . 3D case. Let C = π 2 M 1 M 2 , P = M 1 M 2 ,
and, for any p = 0, . . . , M 1 -1, q = 0, . . . , M 2 -1,

(2.14)                          α R,q+M 2 p (k) = sin θ p φ 3 R,a (k • e θp,ϕq ) , α R,q+M 2 p (l) = ψ 3 R,b (Π e ⊥ θp,ϕq (l)) , θ p = p π M 1 , ϕ q = π q M 2
, e θp,ϕq = (sin θ p cos ϕ q , sin θ p sin ϕ q , cos θ p ) ,

φ 3 R,a (s) = R -R a (|ρ|) |ρ|e iρ s dρ , φ 3 R,b (s) = R -R b (|ρ|) |ρ|e iρ s dρ , ψ 3 R,b (Π e ⊥ θ,ϕ (l)) = π 0 φ 3 R,b Π e ⊥ θ,ϕ (l) cos θ dθ .
In the specific case of 3D hard spheres (γ = 1, b ≡ 1), we get

B(z, y) = 4C Φ , a(|z|) = 4C Φ , b(|y|) = 1 .
This gives

β R (k, l) = π 2 M 1 M 2 P -1 p = 0 α R,p (k) α R,p (l) , with M = M 1 M 2 , α R,p (k) and α R,p (l) computed as follows: α R,q+M 2 p (k) = 4C Φ sin θ p φ 3 R (k • e θp,ϕq ) , α R,q+M 2 p (l) = ψ 3 R Π e ⊥ θp,ϕq (l) , φ 3 R (s) = R 2 2 Sinc(sR) -Sinc 2 sR 2 , ψ 3 R (s) = 2R 2 Sinc 2 Rs 2 .
Both these formula allow for an overall computational cost of O(N 2dv log 2 N ). The practical details of the implementation are summarized in the appendix.

A rescaling velocity model of the BNE

We have seen in Section 1.2 that the BNE can exhibits blowup in the quantum degenerative fermionic case, but with no a priori knowledge on the rate of explosion. In order to simulate numerically such result we propose to adopt the rescaling velocity framework that was developed in the series of papers [START_REF] Filbet | A rescaling velocity method for kinetic equations: the homogeneous case[END_REF][START_REF] Filbet | A Rescaling Velocity Method for Dissipative Kinetic Equations -Applications to Granular Media[END_REF][START_REF] Rey | A Rescaling Velocity Method for Collisional Kinetic Equations -Application to Flocking[END_REF]. It consists in rescaling dynamically the equation, hence allowing to follow very accurately the spread and concentration of the solution to the equation.

3.1.

A general scaling framework. We now introduce a given function ω : R + × Ω x → R * + . Such function has to be positive, and represents accurately the support in velocity of the distribution f (t, x, •) solution to (1.1). Candidates for ω will be presented later. At present time, we focus on the function G : R + × Ω x × R dv that is introduced as follows: for any (t,

x, v) ∈ R + × Ω x × R dv , (3.1) f (t, x, v) = µ π ω(t, x) L dv G t, x, π ω(t, x) L (v -λ u(t, x)) ,
with µ > 0, L > 0 and λ ∈ {0, 1} to be precised later. We can remark that

ρ(t, x) = R dv f (t, x, v) dv = µ R dv G(t, x, ξ) dξ , R dv ξ G(t, x, ξ) dξ = 0 ,
for all (t, x) ∈ R + × Ω x . We shall exhibit the evolution equation that is satisfied by G by injecting (3.1) in (1.1). We first define the change of variables

ξ(t, x, v) = π ω(t, x) L (v -λ u(t, x)) , v(t, x, ξ) = L ξ π ω(t, x)
+ λ u(t, x) , so we have

f (t, x, v) = µ π ω(t, x) L dv G (t, x, ξ(t, x, v)) ⇔ G(t, x, ξ) = µ -1 π ω(t, x) L -dv f (t, x, v(t, x, ξ)) .
The time derivative of f is rewritten as

∂ t f (t, x, v) = µ π ω(t, x) L dv d v ∂ t ω(t, x) ω(t, x) G (t, x, ξ(t, x, v)) + (∂ t G) (t, x, ξ(t, x, v)) + ∂ t ω(t, x) ω(t, x) ξ(t, x, v) • (∇ ξ G) (t, x, ξ(t, x, v)) -λ π ω(t, x) L ∂ t u(t, x) • (∇ ξ G) (t, x, ξ(t, x, v)) . (3.2)
The free transport term is rewritten as

v • ∇ x f (t, x, v) = L π ω(t, x) ξ(t, x, v) + λ u(t, x) µ π ω(t, x) L dv • d v ω(t, x) ∇ x ω(t, x) G (t, x, ξ(t, x, v)) + (∇ x G) (t, x, ξ(t, x, v)) + ∂ x j ω(t, x) ω(t, x) ξ(t, x, v) • (∇ ξ G) (t, x, ξ(t, x, v)) j = 1,...,dv - π ω(t, x) L λ ∂ x j u(t, x) • (∇ ξ G) (t, x, ξ(t, x, v)) j = 1,...,dv .
(

Concerning the rewritting of the collision operator Q α (f ), we first remark that

B (|v -v * |, cos θ) = π ω(t, x) L -γ B (|ξ(t, x, v) -ξ(t, x, v * )| , cos θ) , for any (v, v * , θ) ∈ R dv × R dv × R.
In the same spirit, we can remark that, with the convention (1.1), we have

f (t, x, v ) = π ω(t, x) L dv µ G t, x, ξ(t, x, v) + ξ(t, x, v * ) 2 + |ξ(t, x, v) -ξ(t, x, v * )| 2 σ , f (t, x, v * ) = π ω(t, x) L dv µ G t, x, ξ(t, x, v) + ξ(t, x, v * ) 2 - |ξ(t, x, v) -ξ(t, x, v * )| 2 σ , .
Consequently, the collision operator writes as

(3.4) Q α (f )(t, x, v) = π ω(t, x) L dv-γ µ 2 Q αµ( πω(t,x) L ) dv (G) (t, x, ξ(t, x, v)) .
Gathering (3.2), (3.3) and (3.4), we get the evolution equation for G:

∂ t G + ∇ ξ • ∂ t ω ω ξ -λ π ω L ∂ t u G G + τ L π ω ξ + λ u • d v ω ∇ x ω G + ∇ x G + ∂ x j ω ω ξ • ∇ ξ G j = 1,...,dv - π ω L λ ∂ x j u • ∇ ξ G j = 1,...,dv = c µ π ω L -γ Q αµ( πω L ) dv (G) (3.5)
3.2. Application to the space homogeneous BNE. We now focus on the following homogeneous quantum Boltzmann equation written under divergence form:

           ∂ t G + ∇ ξ • ∂ t ω ω ξ -λ π ω L ∂ t u G = c µ π ω L -γ Q R αµ( πω L ) dv (G) , G(t = 0, ξ) = 1 µ L π ω 0 dv f 0 L π ω 0 ξ + u 0 ,
At this point, the classical method can be distinguished from the velocity rescaling method as follows:

The classical method. It is characterized by the following parameter choice:

λ = 0 , ω ≡ 1 , µ = L π dv , so the equation to be discretized on R + × [-π, π] dv is      ∂ t G = c π L -dv-γ Q R α (G) , G(t = 0, ξ) = f 0 L π ξ ,
and

f is obtained from G on R + × Ω v with Ω v = [-L, L] dv .
The velocity rescaling method. It is characterized by

λ = 1 , µ = 1 ,
and ω to be identified later, so the equation in G to be solved on

R + × [-π, π] dv is            ∂ t G + ∇ ξ • ∂ t ω ω ξ - π ω L ∂ t u G = c π ω L -γ Q R α( πω L ) dv (G) , G(t = 0, ξ) = L π ω 0 dv f 0 L π ω 0 ξ + u 0 , and f (t, •) is obtained from G(t, •) on Ω v (t) = u(t) + - L ω(t) , L ω(t) dv for any t ∈ R + .
3.3. Numerical scheme for the velocity rescaled BNE. We assume from now that the support of

G(t = 0, •) is included in B(0, S) ⊂ [-π, π] dv with S > 0 and R > 0 such that R ≥ 2S , π ≥ R + (1 + √ 2)S 2 .
We consider a time grid (t n ) n ∈ N ⊂ R + with ∆t n = t n+1 -t n > 0, t 0 = 0 and G n : [-π, π] dv → R being an approximation of G(t n , •). Finally, we define ∆v n and ∆ξ as

∆v n = 2L ω n N dv , ∆ξ = 2π N dv ,
and we consider a space discretization of G n with

G n (ξ j ) ≈ G n j , ξ j d ,d = 2j d π N , j d = - N 2 , . . . , N 2 -1 , d = 1, . . . , d v ,
and

G 0 j := G(0, ξ j ) = 1 µ L π ω 0 dv f 0 L π ω 0 ξ j + λ u 0 . Defining f n h = (f n j ) j ∈ I N with f n j := µ π ω n L dv G n j ,
and

v n j d ,d = λ u n d + 2L j d N ω n for any j d = - N 2 , . . . , N 2 -1, d = 1, . . . , d v (u n
d being an approximation of u d (t n ), to be defined later), we have

G(t n , ξ j ) ≈ G n j ⇐⇒ f (t n , v n j ) ≈ f n j
for any j ∈ I N and n ∈ N. This provides an approximation of f (t, •) on Ω v (t).

Remark 6 . In most cases, it is possible to identify the limit state

f ∞ such that f → f ∞ as t → +∞.
For this purpose, we refer to Section 1.2 in order to identify the expression of f ∞ . This function is characterized thanks to the moments of f 0 (mass ρ 0 , averaged velocity u 0 , internal energy e 0 ...).

Assuming that f ∞ is analytically known, we build

f ∞,n h = (f ∞,n j ) j ∈ I N for any n such that f ∞,n j = f ∞ (v n j ) , which gives an approximation of f ∞ on Ω n v . We can also introduce G ∞,n h = (G ∞,n j ) j ∈ I N as G ∞,n j = 1 µ L π ω n dv f ∞,n j , ∀ j ∈ I N .
We describe now the procedure associated to a forward Euler time-semi-discretization:

(1) We assume that G n h (or equivalently f n h ) is known along with the mass ρ n , the averaged velocity u n , the internal energy e n , and the velocity scale ω n at time step n.

(2) We compute Q R αµ( πω n L ) dv (G n h ) j for any j ∈ I N : (a) We compute the Fourier modes of G n h and get

G n k = 1 N j ∈ I N G n j e -2iπ k•j N , ∀k ∈ I N . (b) We compute Q R i,c (G n h , G n h ) j (i = 1, 2) and Q R i,q (G n h , G n h , G n h ) j (i = 1, 2, 3, 4) for any j ∈ I N . (c) We get Q R αµ( πω n L ) dv (G n h ) j as Q R αµ( πω n L ) dv (G n h ) j =Q R 1,c (G n h , G n h ) j -Q R 2,c (G n h , G n h ) j -αµ( πω n L ) dv Q R 1,q (G n h , G n h , G n h ) j + Q R 2,q (G n h , G n h , G n h ) j -Q R 3,q (G n h , G n h , G n h ) j -Q R 4,q (G n h , G n h , G n h ) j .
(3) If we do not assume ω ≡ 1 (velocity rescaling case), we do the following computations:

(a) We estimate ρ n+1 , u n+1 and e n+1 with

ρ n+1 -ρ n ∆t n = c ∆v n j ∈ I N Q R α (f n h ) j = c∆ξ πω n L -γ µ 2 j ∈ I N Q R αµ( πω n L ) dv (G n h ) j ; ρ n+1 u n+1 -ρ n u n ∆t n = c ∆v n j ∈ I N v n j Q R α (f n h ) j = c ∆ξ πω n L -γ µ 2 j ∈ I N L π ω n ξ j + λ u n Q R αµ( πω n L ) dv (G n h ) j ; ρ n+1 e n+1 -ρ n e n ∆t n = c 2 ∆v n j ∈ I N |v n j -u n | 2 Q R α (f n h ) j = c 2 ∆ξ πω n L -γ µ 2 j ∈ I N L π ω n ξ j + (λ -1) u n 2 Q R αµ( πω n L ) dv (G n h ) j .
(b) We then compute the temperature T n+1 . Since it is trivial for the classical case α = 0 (T n+1 = 2/d v e n+1 ), let us focus on the quantum case α = 0. The temperature T n+1 is obtained as follows: test the inequality

|α| ρ n+1 d v 4πe n+1 dv 2 ∈ I eq ,
where I eq is defined by (1.11). If it is true, T n+1 is solved along with the fugacity z n+1 by solving the following nonlinear system:

             ρ n+1 = (2πT n+1 ) dv 2 |α| K dv 2 (z n+1 ) , e n+1 = d v T n+1 2 K dv 2 +1 (z n+1 ) K dv 2 (z n+1
) .

If it is false, we distinguish 2 cases:

• For 3D Bose-Einstein case, condensation occurs, so the fugacity z n+1 is equal to 1 and T n+1 is explicitly computed as follows:

T n+1 = 2 |α| ρ n+1 e n+1 3 ζ(5/2) (2π) 3/2 2/5
, and the critical mass m n+1 0 is obtained as

m n+1 0 = ρ n+1 - (2πT n+1 ) 3/2 |α| ζ(3/2) .
• For Fermi-Dirac case (2D or 3D), we test the equality

|α| ρ n+1 d v 4πe n+1 dv 2 =      2 ,
for 2D Fermi-Dirac case, If it is true, then z n+1 = +∞ and T n+1 = 0. If it is false, we cannot say anything about the values of z n+1 and T n+1 . (c) We compute an estimate of ω n+1 :

• By using the temperature (quantum case α = 0):

ω n+1 = √ T n+1 ,
• By using the kinetic energy (classical case α = 0):

ω n+1 = T n+1 + |u n+1 | 2 d v . (4) We compute ω n+1 -ω n ∆t n ω n ξ j -λ π ω n L u n+1 -u n ∆t n
G n j for any j ∈ I N , then its Fourier modes, and finally get

∇ ξ • ω n+1 -ω n ∆t n ω n ξ -λ π ω n L u n+1 -u n ∆t n G n h j
for any j ∈ I N .

(5) We finally compute G n+1 h by solving

G n+1 j -G n j ∆t n + ∇ ξ • ω n+1 -ω n ∆t n ω n ξ -λ π ω n L u n+1 -u n ∆t n G n h j = cµ πω n L -γ Q R αµ( πω n L ) dv (G n h ) j ,
and can reconstruct the approximation

f n+1 h of f (t n+1 , •) on Ω n+1 v = u n+1 + - L ω n+1 , L ω n+1 thanks to the relation f n+1 j = µ πω n+1 L dv G n+1 j , ∀ j ∈ I N .
(6) If we assume that ω ≡ 1 and λ = 0 (velocity rescaling is not considered), we compute ρ n+1 , u n+1 and e n+1 as follows:

ρ n+1 = ∆v j ∈ I N f n+1 j = ∆ξ µ j ∈ I N G n+1 j , ρ n+1 u n+1 = ∆v j ∈ I N v n+1 j f n+1 j = ∆ξ µ j ∈ I N L π ξ j G n+1 j , ρ n+1 e n+1 = ∆v 2 j ∈ I N |v n+1 j -u n+1 | 2 f n+1 j = ∆ξ 2 µ j ∈ I N L π ξ j 2 G n+1 j .
It is also possible to consider a RK2 time semi-discretization by using cleverly some parts of the explicit Euler time loop above.

Numerical results

We shall now present some numerical simulations of our fast spectral method for the BNE. For a given discrete solution f n h of the usual variables or G n h in rescaled variables, we aim to compute the following quantities:

• The mass ρ n , the averaged velocity u n , the kinetic energy E n c , and the internal energy e n :

ρ n := ∆v n j ∈ I N f n j = ∆ξ µ j ∈ I N G n j , u n := 1 ρ n ∆v n j ∈ I N v n j f n j = 1 ρ n ∆ξ µ j ∈ I N L π ω n ξ j + λ u n G n j , E n c := ∆v n j ∈ I N |v n j | 2 f n j = ∆ξ µ j ∈ I N L π ω n ξ j + λ u n 2 G n j , e n := 1 2ρ n ∆v n j ∈ I N |v n j -u n | 2 f n j = 1 2ρ n ∆ξ µ j ∈ I N L π ω n ξ j + (λ -1) u n 2 G n j ;
• The stress-energy tensor

T n = (T n d,d ) d,d ∈ dv (R): T n d,d := 1 ρ n ∆v n j ∈ I N (v n j d ,d -u n d ) (v n j d ,d -u n d ) f n j = 1 ρ n ∆ξ µ j ∈ I N L π ω n ξ j d ,d + (λ -1) u n d L π ω n ξ j d ,d + (λ -1) u n d G n j ;
• The entropy for the classical case H n c :

H n c := ∆v n j ∈ I N f n j ln(f n j ) = µ ∆ξ j ∈ I N G n j ln   µ π ω n L dv G n j   ;
• The entropy for the quantum case H n q :

H n q := ∆v n

j ∈ I N f n j ln(f n j ) + 1 -α f n j α ln 1 -α f n j = ∆ξ µ j ∈ I N G n j ln   µ πω n L dv G n j   + µ -1 πω n L -dv -α G n j α ln   1 -α µ πω n L dv G n j   ;
• p norms in space, for p ∈ (1, ∞):

f n h 2 =   ∆v n j ∈ I N |f n j | p   1/p =   ∆ξ µ 2 π ω n L dv j ∈ I N G n j p   1/p . 4.1.
Validation of spectral accuracy. The first test sequence is dedicated to the validation of the spectral accuracy of the numerical methods presented in Section 3.2 in the non-degenerative case.

To do this, we proceed as follows:

(1) Set the initial state f 0 = f (t = 0, •) as a candidate for the limit state;

(2) Discretize it on a velocity grid to get f 0 h and compute the discrete moments ρ 0 h , u 0 h , e 0 h ; (3) If the velocity rescaling method is used, we set u 0 h = u 0 where u 0 is the exact mean velocity associated to f 0 , and deduce the discrete temperature T 0 h (and the discrete fugacity z 0 h for the quantum case) and the discrete velocity scale ω 0 h (default value is

ω 0 h = 1); (4) Set ρ ∞ = ρ 0 h , u ∞ = u 0 , T ∞ = T 0 h , z ∞ = z 0 h , ω ∞ = ω 0 h and define the associated limit state f ∞ ; (5) Discretize f ∞ to get f ∞ h .
To highlight the spectral accuracy of our scheme, we compute the value of the numerical residual Q(f ∞ h ) ∞ as the velocity grid is refined, over the simulation domain

Ω h = λ u 0 h + - L ω 0 h , L ω 0 h ,
where λ = 1 if the velocity rescaling is activated and λ = 0 in the contrary case. We also choose the truncation parameters R and S such that

R = S , π = (3 + √ 2) S 2 ,
so the artificial truncation of f 0 depends on the value of L provided by the user.

2D quantum non-degenerative case.

For testing the accuracy of the method in the 2D quantum case, we select the maxwellian function M q as the initial state. Since such tests were already performed in [START_REF] Filbet | A numerical scheme for the quantum boltzmann equation efficient in the fluid regim[END_REF], we chose to focus on test cases involving = 3 for highlighting the quantum effects in the collision operator. In order to avoid Fermi-Dirac quantum degeneracy, we define f 0 as

f 0 (v) = 1 |α| 1 z -1 exp( |v-u 0 | 2 2σ ) + sgn(α) , α = + dv , for Fermi-Dirac particles, -dv ,
for Bose-Einstein particles, where σ is set to 1 or 0.5. In addition, since we want to present the advantages of the velocity rescaling, we enrich the test list by considering 2 values of the initial average velocity u 0 :

ρ 0 = 1 , u 0 = 0 , or u 0 = u b = 8 3 √ 2 + 2 (1, 1) T ,
and several values of the limit L of the simulation bounding box. According to Tables 1234where the numerical residual Q R α (f ∞ h ) ∞ is computed for various values of L, σ, u 0 and grid sizes, it appears that our method is spectrally accurate if we consider a simulation bounding box with L well suited to the moments u 0 and T 0 = σ. This is highlighted for the 2D Bose-Einstein case with small values of σ: indeed, such choice for σ leads to a distribution f 0 that is very sharp and concentrated close to u 0 (see the profiles v x → f 0 (v x , 0) in Figures 4.2). Consequently, the spectral accuracy becomes very difficult to observe in such configuration even with the velocity rescaling method (see Table 3).

In addition, we remark that the velocity rescaling method provides much smaller residuals than the classical method for most cases that have been studied. This is not surprising for the case where u 0 = u b because of the support B(0, S) of the initial distribution that leads to the omission of a part of f 0 h . For the case where u 0 = 0, this can be explained by the fact that the rescaled distribution G 0 h is exactly centered on 0 without any discretization error which is not the case for the classical method.

4.1.2. 3D quantum non-degenerative case. We then perform a similar test sequence for the 3D nondegenerative quantum case. More precisely, for the Fermi-Dirac case, we chose the following parameters:

ρ 0 = 1 , u 0 = 0 or u 0 = u b = 4 3 + √ 2 1, 1, √ 2 T , σ ∈ {0.5, 1} ,
and for the Bose-Einstein case, we chose

u 0 = 0 or u 0 = u b = 4 3 + √ 2 1, 1, √ 2 T , σ ∈ {0.5, 1} .
In order to avoid Bose-Einstein condensation (cases with m 0 > 0), we chose ρ 0 = 0.5 in Tables 78and ρ 0 = 0.2 in Tables 9-10.

Since the computational cost is multiplied by 64 each time we multiply by 2 the number of grid points in each direction, we skipped the tests for 128 3 and 256 3 grids (see Tables 5-6-7-8-9-10).

Even without these highly refined grids, we can formulate some conclusions that complete those we have made for the 2D non-degenerative case: if we consider a case that is reasonably far from quantum degeneracy (Tables 5-6-9-10), and up to a good choice for L according to the prescribed values of ρ 0 , u 0 and σ, both classical and velocity rescaling methods are spectrally accurate when the velocity grid is refined. In addition, the velocity rescaling method provides smaller residuals Q R α (f ∞ h ) ∞ than the classical method for almost all tests. This last remark will be very useful for future 3D tests since it allows to reduce the size of velocity grid and, consequently, the computational cost of the 3D time-space numerical scheme.

L Grid u 0 = 0 u 0 =
Q R α (f ∞ h ) ∞ for 2D Fermi-
We can formulate an additional conclusion that is specific to the test case close to Bose-Einstein

L Grid u 0 = 0 u 0 =
Q R α (f ∞ h ) ∞ for 2D Fermi-Dirac non-degenerative case with = 3, ρ 0 = 1 and σ = 1 (z = 3.1887).

L

Grid u 0 = 0 u 0 = u b with rescaling without rescaling with rescaling without rescaling L = 4.0 σ = 0.5 16 2 1.26588e+02 1.10944e+03 1.59361e+02 6.13676e+01 32 2 2.25994e+00 4.33543e+01 2.24858e+00 1.25675e+02 64 2 6.00369e-02 1.86373e+00 6.00369e-02 7.76987e+00 128 2 8.25844e-04 7.93633e-03 8.25844e-04 1.82181e-01 256 2 6.20715e-08 4.30685e-08 5.94336e-08 1.62130e-05

L = 6.0 σ = 0.5 16 2
2.46280e+05 2.34093e+04 2.54047e+05 3.67843e+02 32 2 1.38535e+01 2.68651e+02 4.19147e+00 9.30524e+00 64 2 6.82833e-01 7.49992e+00 6.81247e-01 1.41910e+01 128 2 1.99389e-02 1.31040e-01 1.07615e-02 1.72672e-01 256 2 1.50878e-05 3.03279e-05 2.88911e-05 6.28249e-05 Table 3. Spectral accuracy. Values of the numerical residual Q R α (f ∞ h ) ∞ for 2D Bose-Einstein case with = 3, ρ 0 = 1 and σ = 0.5 (z = 0.943). condensation (see Tables 78): in such case, the discrete distribution f h is almost singular in a neighborhood of u 0 , so the computation of Q R α (f ∞ h ) ∞ may produce NaN values. This phenomenon is particularly significant when we set u 0 = 0 or when we activate the velocity rescaling method. Indeed, in such case, ξ = 0 is a grid point and G h (0) becomes very large. The only studied way to bypass this problem is to deactivate the velocity rescaling method and to choose u 0 outside of the velocity grid points. 4.2. Relaxation towards equilibrium. We have seen in Section 1.2 that the distribution function f (t, •) is expected to converges to an equilibrium state f ∞ as t → +∞ . We present here some tests

L Grid u 0 = 0 u 0 =
Q R α (f ∞ h ) ∞
for 2D Bose-Einstein case with = 3, ρ 0 = 1 and σ = 0.5 (z = 0.943) and σ = 1 (z = 0.7613).

L

Grid u 0 = 0 u 0 = u b with rescaling without rescaling with rescaling without rescaling L = 4.0 σ = 0.5 16 3 6.65068e-04 2.64044e-01 6.65068e-04 4.38601e-02 32 3 1.38282e-05 3.51620e-03 1.38282e-05 3.47688e-04 64 3 4.02472e-09 3.11989e-06 4.02472e-09 8.48441e-08 L = 6.0 σ = 0.5 16 3 1.46161e-02 4.93491e-01 1.46161e-02 6.27339e-01 32 3 4.64216e-04 1.27742e-02 4.64216e-04 1.44628e-02 64 3 5.17232e-07 3.23812e-06 5.17232e-07 1.34607e-05 L = 8.0 σ = 0.5 16 3 6.86816e-02 3.05672e+00 6.86816e-02 3.86836e+00 32 3 3.63615e-03 1.49763e-01 3.63615e-03 1.29305e-01 64 3 2.29417e-05 3.12185e-04 1.55549e-05 2.64884e-04 L = 10.0 σ = 0.5 16 3 6.02018e-02 5.16432e+01 6.02018e-02 3.06977e+01 32 3 1.57951e-02 9.41582e-01 1.57951e-02 6.93264e-01 64 3 2.09198e-04 6.20599e-03 2.09198e-04 5.89099e-03 Table 5. Spectral accuracy. Values of the numerical residual Q R α (f ∞ h ) ∞ for 3D Fermi-Dirac non-degenerative case with = 3, ρ 0 = 1 and σ = 0.5 (z = 2.43228e+01). that highlight this phenomenon for all treated cases and provide a numerical validation of the socalled Fermi-Dirac saturation and Bose-Einstein condensation. For this, let us associate to a given discrete initial distribution f 0 h (and its related macroscopic quantities) a threshold * defined as

* =   4πe 0 h d v dv 2 1 ρ 0 h   1 dv
According to Section 1.2, we expect the following behavior of the discrete distribution f h as t → +∞:

• 2D Bose-Einstein particles: f h converges to the quantum maxwellian defined as

f ∞ (v) = M q (v) = 1 2 1 1 z 0 h exp |v-u 0 h | 2 2T 0 h -1 , L
Grid u 0 = 0 u 0 = u b with rescaling without rescaling with rescaling without rescaling L = 4.0 σ = 1.0 16 3 9.28576e-04 4.02042e-02 9.28576e-04 4.32287e-03 32 3 5.11585e-05 8.89902e-05 5.11585e-05 9.91060e-06 64 3 2.93880e-08 1.19765e-07 2.27659e-08 1.72883e-07 L = 6.0 σ = 1.0 16 3 3.60402e-04 8.95561e-02 3.60402e-04 1.18694e-01 32 3 9.32472e-07 1.70466e-04 9.32472e-07 2.05708e-04 64 3 4.40082e-11 2.41131e-08 3.46700e-11 7.59905e-09 L = 8.0 σ = 1.0 16 3 4.72058e-04 4.90556e-01 4.72058e-04 2.76318e-01 32 3 1.79828e-06 1.73790e-03 1.79828e-06 2.09096e-03 64 3 5.38083e-11 5.71195e-08 5.82632e-11 1.47296e-07 L = 10.0 σ = 1.0 16 3 6.82783e-03 6.30273e+00 6.82783e-03 4.16489e+00 32 3 2.16863e-05 2.08866e-02 2.16863e-05 9.24323e-03 64 3 2.20701e-09 4.38071e-06 2.20701e-09 4.40874e-06 Table 6. Spectral accuracy. Values of the numerical residual Q R α (f ∞ h ) ∞ for 3D Fermi-Dirac non-degenerative case with = 3, ρ 0 = 1 and σ = 1 (z = 3.09922).

L

Grid u 0 = 0 u 0 = u b with rescaling without rescaling with rescaling without rescaling

L = 4.0 σ = 0.5 16 3 - - - 2.01861e+01 32 3 - - - 3.07997e+02 64 3 - - - 3.30064e+01 L = 6.0 σ = 0.5 16 3 - - - 1.33886e+03 32 3 - - - 5.39380e+02 64 3 - - - 1.83015e+02 L = 8.0 σ = 0.5 16 3 - - - 2.09609e+04 32 3 - - - 3.26936e+02 64 3 - - - 1.87903e+02 L = 10.0 σ = 0.5 16 3 - - - 5.58828e+03 32 3 - - - 2.73300e+01 64 3 - - - 5.21282e+04
Table 7. Spectral accuracy. Values of the numerical residual Q R α (f ∞ h ) ∞ for the 3D Bose-Einstein non-degenerative case with = 3, ρ 0 = 0.5 and σ = 0.5 (z = 0.99706).

where (z 0 h , T 0 h ) solves

       ρ 0 h = 2πT 0 h 2 B 1 (z 0 h ) , e 0 h = T 0 h B 2 (z 0 h ) B 1 (z 0 h ) ,
• 3D Bose-Einstein particles: f h converges to the limit state f ∞ defined as

f ∞ (v) =        M q (v) , if m 0,h = ρ 0 h - 1 3 4π e 0 h 3 3/2 ζ(3/2) 5/2 ζ(5/2) 3/2 ≤ 0, M q (v) , else,

L

Grid u 0 = 0 u 0 = u b with rescaling without rescaling with rescaling without rescaling L = 4.0 σ = 1.0 16 3 4.85493e-02 1.23198e+01 4.85493e-02 6.20351e-01 32 3 6.60721e-04 6.57365e-02 6.60721e-04 1.40535e-03 64 3 4.36371e-07

Q R α (f ∞ h ) ∞
for the 3D Bose-Einstein non-degenerative case with = 3, ρ 0 = 0.5 and σ = 1 (z = 0.63071).

L

Grid u 0 = 0 u 0 = u b with rescaling without rescaling with rescaling without rescaling L = 4.0 σ = 0.5 16 3 1.89242e-02 5.59497e+00 1.89242e-02 6.51227e-01 32 3 1.54426e-04 3.17405e-02 1.54426e-04 2.91398e-03 64 3 6.62168e-08 7.87246e-06 6.62168e-08 2.68118e-06 L = 6.0 σ = 0.5 16 3 6.16653e-01 2.43214e+01 6.16653e-01 9.89261e+00 32 3 3.25459e-03 2.18171e-01 3.25459e-03 4.97853e-01 64 3 7.62635e-06 1.16549e-04 7.62635e-06 7.32149e-04 L = 8.0 σ = 0.5 16 3 2.04123e+03 2.46424e+02 2.04123e+03 7.08188e+01 32 3 2.43393e-02 2.03212e+00 2.43393e-02 6.57090e-01 64 3 1.51689e-04 3.61768e-03 1.51689e-04 2.87658e-03 L = 10.0 σ = 0.5 16 3 -7.36733e+03 -1.07883e+02 32 3 1.38333e-01 9.78043e+00 1.38333e-01 7.43831e-01 64 3 1.03698e-03 4.39372e-02 1.03698e-03 4.01027e-02 Table 9. Spectral accuracy. Values of the numerical residual Q R α (f ∞ h ) ∞ for the 3D Bose-Einstein non-degenerative case with = 3, ρ 0 = 0.2 and σ = 0.5 (z = 0.68499).

where, in the former case,

M q (v) = 1 3 1 1 z 0 h exp |v-u 0 h | 2 2T 0 h -1 , with (z 0 h , T 0 h ) solving          ρ 0 h = (2πT 0 h ) 3/2 3 B 3/2 (z 0 h ) , e 0 h = 3T 0 h 2 B 5/2 (z 0 h ) B 3/2 (z 0 h ) ,
and, in the latter case,

M q (v) = m 0,h δ 0 (v -u 0 h ) + 1 3 1 exp |v-u 0 h | 2 2T 0 h -1 , T 0 h = 2 e 0 h ζ(3/2) 3ζ(5/2)
. L Grid u 0 = 0 u 0 = u b with rescaling without rescaling with rescaling without rescaling L = 4.0 σ = 1.0 16 3 3.74496e-04 2.77192e-02 3.74496e-04 5.36638e-03 32 3 4.12493e-07 5.76615e-06 4.12493e-07 7.85561e-07 64 3 6.86211e-12 3.46635e-07 6.86211e-12 1.63647e-07 L = 6.0 σ = 1.0 16 3 1.49832e-04 7.65264e-02 1.49832e-04 1.11485e-01 32 3 4.51218e-08 1.52468e-05 4.51218e-08 6.22715e-05 64 3 2.65642e-13 1.00237e-09 2.65642e-13 4.21875e-09 L = 8.0 σ = 1.0 16 3 5.11241e-04 4.59837e-01 5.11241e-04 5.35585e-01 32 3 3.22287e-07 2.73180e-04 3.22287e-07 4.08919e-04 64 3 3.36264e-12 2.79570e-09 3.36264e-12 9.71734e-09 L = 10.0 σ = 1.0 16 3 2.78127e-03 2.60318e+00 2.78127e-03 1.18199e+00 32 3 5.03331e-06 4.81342e-03 5.03331e-06 1.81840e-03 64 3 2.03119e-10 1.93959e-07 2.03099e-10 2.62032e-07

Table 10. Spectral accuracy. Values of the numerical residual Q R α (f ∞ h ) ∞ for the 3D Bose-Einstein non-degenerative case with = 3, ρ 0 = 0.2 and σ = 1 (z = 0.30354).

• Fermi-Dirac particles: f h converges to the limit state f ∞ defined as

f ∞ (v) =        M q (v) , if m 0,h = ρ 0 h - K dv dv 4π e 0 h d v dv 2 < 0, χ(v) , if m 0,h = 0, with K 2 = 2, K 3 = 5 3 10 π
, and where

M q (v) = 1 dv 1 1 z 0 h exp |v-u 0 h | 2 2T 0 h + 1 , with (z 0 h , T 0 h ) solving            ρ 0 h = (2πT 0 h ) dv 2 dv F dv 2 (z 0 h ) , e 0 h = d v T 0 h 2 F dv 2 +1 (z 0 h ) F dv 2 (z 0 h )
, and

χ(v) = 1 dv 1 B(u 0 h ,A 0 h ) (v) , A 0 h = 2e 0 h (d v + 2) d v .
4.2.1. 2D Fermi-Dirac relaxation. We consider a 128 × 128 uniform grid in the velocity domain [-8, 8] 2 and a time step ∆t = 0.025. We do not use velocity rescaling by taking (λ, ω, µ) = (0, 1, ( 8 π ) 2 ) and use a Runge-Kutta 2 SSP time integrator. We consider 2 different values for f 0 :

• Fermi-Dirac saturation: we set (4.1)

f 0,1 (v) = ρ 0 π(A 0 ) 2 1 B(u 0 ,A 0 ) (v) , with A 0 = 2 √ e 0 .
Even if the initial distribution is not continuous, we ensure that the entropy is always defined since the maximal value of f 0 is ρ 0 4πe 0 ≈ ( * ) -2 . Hence we expect a relaxation towards a quantum maxwellian M q for any < * , and no time variation when = * (or numerical instabilities leading to a non-physical phenomenon).

r

-2 = (r * ) -2 0.1 7.95775 0.2 1.98944 0.5 3.1832e-01 0.8 1.2434e-01 0.9 9.82438e-02 0.95 8.81745e-02 0.99 8.11932e-02 1 7.95775e-02 Table 11. Values of -2 when the initial distribution is a classical maxwellian.

• Classical maxwellian distribution: we set

(4.2) f 0,2 (v) = ρ 0 2πσ 0 exp - |v -u 0 | 2 2σ 0 ,
with ρ 0 = 1, u 0 = (0, 0) and σ 0 = 1. In that case, one has * ≈ 3.5449. The motivation for choosing such initial datum is that it is a regular distribution whose maximal value is equal to ρ 0 2πσ 0 ≈ 0.15915. Let us remind the reader that an entropic solution must satisfy max(f ) ≤ ( * ) -2 , and that entropy decay is not ensured in the opposite case. Consequently, we must pay attention to the chosen value for for running the simulation with respect to the value of ( * ) -2 ≈ 0.07958 prescribed by the moments of f 0 since taking too large will break the entropy definition criterion above (even if we do not reach the saturation criterion = * ).

Starting with f 0,1 from (4.1), we indeed observe a relaxation of the distribution to a quantum maxwellian for any = r * with r ∈ {0.1, 0.5, 0.8, 0.9, 0.95, 0.99} (see . We also notice that the shape of the distribution at time t = 30 is close to an indicator function of the form χ as we consider values of r that are close to 1 -. We observe that for the particular case r = 0.99, the entropy increases after a finite time. This is not surprising since the limit state as t → +∞ is almost discontinous, so the Fourier transforms applied to such data induces Gibbs numerical phenomenon.

We then run simulations with r ≥ 1 for investigating the behavior of the distribution when Fermi-Dirac saturation occurs. For the specific case r = 1 where we expect a time independent solution, the errors due to the treatment of a non-continous initial datum by Fourier transform make the approximate solution blow up numerically within a finite time (see Figure 4.5). For higher values of r, we also observe a numerical blow up of the approximate solution within a finite time which seems to depend on . This last observation is not surprising since we considered the case where ρ 0 h > 4πe 0 h 2 where we cannot identify a limit state as t → +∞.

We then start with the classical maxwellian f 0,2 from (4.2). In that case, since * ≈ 3.54491 and max(f 0,2 ) ≈ 0.15915, the entropy definition criteria may be broken at initial time. Taking = r * with r ≥ 0 indicates that it is impossible to define the entropy for r ≥ 0.8 (see Table 11).

By taking r ∈ {0.1, 0.5, 0.8}, the entropy is correctly defined at any time since t → max(f (t, •)) is maximal at time t = 0. In such case, we observe the relaxation of the distribution to the associated limit state f ∞ = M q as t → +∞ (see Figure 4.6).

For the cases r ∈ {0.9, 0.95, 0.99}, even if the entropy cannot be correctly defined, we observe that the numerical distribution f h still converges to the steady state f ∞ h = M q defined with the moments (ρ 0 h , u 0 h , e 0 h ) for any r < 1 since the relaxation error

f h (t, •) -f ∞ h 1 decreases as t → +∞ (see Figure 4.7).
Finally, the cases r ≥ 1 are similar to the indicator case presented in Figure 4.3: the entropy decay is not observed numerically, because of the oscillations induced by the singular equilibria. The limiting case r = 1 is presented in Figure 4.8, as well as the time evolution of t → max(f h (t, •)). 1,q (G h ), we restrict our simulations to a 32 3 grid in velocity. Again, we consider two definitions of the initial distribution f 0 :

• Fermi-Dirac saturation: we set with A 0 = 10e 0 3 . As for the 2D case, we expect a relaxation towards a quantum maxwellian M q for any < * with ( * ) -3 ≈ 3ρ 0 4π ( 10e 0 3 )

(4.3) f 0,3 (v) = 3ρ 0 4π (A 0 ) 3 1 B(u 0 ,A 0 ) (v) ,
3 2
and no time variation when = * .

• Classical maxwellian distribution: we set Considering the discontinuous initial distribution (4.3), we observe the relaxation of the distribution to a quantum maxwellian for any = r * with r ∈ {0.1, 0.5, 0.8, 0.9, 0.95}. Indeed, we observe a decay in time for the entropy and the relaxation error f h (t, •) -f ∞ h 1 for these values of (see Figures 4.9 Following the methodology proposed for the 2D Fermi-Dirac case, we also tested the case = 0.99 * where we observe again a numerical blow-up of the distribution instead of the relaxation towards a quantum maxwellian. Again, this blow-up is due to Gibbs phenomenon induced by the treatment of an almost discontinous signal with Fourier transforms. In the present case, this numerical blow-up is highlighted by the loss of the decay of the entropy (see Figure 4.11).

(4.4) f 0,4 (v) = ρ 0 (2πσ 0 ) 3 2 exp - |v -u 0 | 2 2σ 0 ,
To complete this study, we also investigated the values = r * with r ≥ 1 and the conclusions we made for the 2D Fermi-Dirac can be extended to the 3D case. For the specific case r = 1 where Fermi-Dirac saturation occurs, the accumulation of local errors due to Gibbs phenomenon quickly deteriorates the discrete distribution f h so the numerical blow-up appears after a small number of and that, in such case, we cannot identify the time behaviour of the distribution (see Figure 4.13 for r = 1.01).

Let us now consider the case of a classical maxwellian f 0,4 . We observe the expected relaxation to the corresponding steady state for = r * for r ∈ {0.1, 0.5, 0.8, 0.9} (see Figures 4.14 Note that, for the latter case, even if the entropy cannot be defined since max(f 0 ) > -3 , we still have the relaxation of the distribution towards to expected steady state.

Due to the coarse velocity mesh, taking higher values of r leads to a numerical blow-up, even if < * (see Figure 4.16). Finally, taking r ≥ 1 leads to a numerical blow-up of the distribution as for the 2D case with similar conditions on (see Figure 4.17).

3D Bose-Einstein relaxation.

We now investigate the relaxation in time of the discrete distribution f h to the quantum M q or the condensate distribution M q according to the value with respect to * = 1 (ρ 0 h ) 1/3 4π e 0 As for the Fermi-Dirac cases, we consider a 32 3 velocity grid on the domain [-L, L] 3 = [-6, 6] 3 and an initial distribution defined as (4.3) or (4.4). This gives us the thresholds * ≈ 3.97285, and * ≈ 4.81755, respectively. We expect that the computed discrete distribution f h relaxes to M q when < * and to M q when ≥ * (i.e. a blow-up of the distribution localized in u 0 h ). The computational limitations we noticed in the previous paragraph for handling 3D quantum cases are more restrictive here than for the 3D Fermi-Dirac case. Indeed, if → ( * ) -, the 32 3 velocity grid we considered is not sufficient to catch the stiff gradients of the limit state f ∞ . We can see in Figure 4.18 that for = r * with r ∈ {0.1, 0.5, 0.8} we can reach successfully the expected limit state as the entropy (and the numerical relaxation error) decrease.

However, for higher values of r ∈ (0.8, 1), the coarseness of the velocity grid deteriorates the distribution profile after a small number of time iterations, even if < * (see Figure 4.19).

Surprisingly, considering ≥ * , i.e. a Bose-Einstein condensation and the associated singular limit state M q , leads to numerical results that are less deteriorated than for r ∈ {0.9, 0.95} (see Figures 4.20). We observe the expected blow-up in time of the discrete distribution f h but it occurs after a simulation time T * that seems to decreases as long as becomes large. This phenomenon may be due to the fact that, in such case, the stiffness of the distribution is due to the Dirac mass and the singular quantum maxwellian in the definition of M q , whereas in the case where → ( * ) -, the stiffness of f ∞ is only due to the quantum maxwellian M q that is almost singular in u 0 h . All these observations are confirmed when we use the classical maxwellian M c as initial datum (Figures omitted for redundancy

Conclusion and perspectives

In this work, we illustrated the efficiency and accuracy of the Fast Spectral algorithm from [START_REF] Filbet | A numerical scheme for the quantum boltzmann equation efficient in the fluid regim[END_REF] along with the rescaling velocity method from [START_REF] Filbet | A Rescaling Velocity Method for Dissipative Kinetic Equations -Applications to Granular Media[END_REF] for approximating the solutions to the Boltzmann-Nordheim equation. In particular, we showed that this numerical method is able to reproduce very accurately some of the main mathematical features of this equation, including the convergence towards singular steady states (namely the phenomena of Fermi-Dirac saturation and Bose-Einstein condensation). As a consequence, we believe that this method is perfectly suited to study some of the still unknown behavior of this model, in particular the rate of convergence towards the (singular) 3D Bose-Einstein condensates, that was conjectured in [START_REF] Escobedo | Finite time blow-up for the boson Nordheim equation[END_REF]. Unfortunately, the accuracy needed for computing such singular solution in 3D is still unavailable in our approach. This motivates us for studying ways to make our simulation code faster in order to consider finer meshes and complete the For any j ∈ I N , we define g R,p,j and g R,p,j on I N such that ( g R,p,j ) k = α R,p (k) e -2iπ k•j N , ( g R,p,j ) k = α R,p (k) e -2iπ k•j N , so we have where g R,p,j and g R,p,j are obtained with inverse Fourier transform of g R,p,j and g R,p,j respectively.

Q R 1,q (G h , G h , G h

Figure 1 . 1 .

 11 Figure 1.1. Graphs of functions µ -→ F ν (e µ ) for ν = 1 (a) and ν = 3 2 (b), and z → B ν (z) for ν = 1 (c) and ν = 3 2 (d). We notice that, for each case, F ν and B ν are valued in the corresponding I eq set defined in (1.11).

2. 4 .

 4 Expansion of the kernel modes. Let us assume from now on that B can be expressed as follows: B(z, y) = a (|z|) b (|y|) , with some given functions a, b : R + → R. Hence, decomposing polar coordinate z = ρ e and y = ρ e with ρ, ρ ∈ [0, R], e, e ∈ S dv-1 , one has

Figure 4 . 2 .

 42 Figure 4.2. Profiles v x → M q (v x , 0) for 2D Bose-Einstein particles with σ = 0.5 (a) and σ = 1 (b). For both profiles, we considered = 3, ρ 0 = 1 and u 0 = 0.

Figure 4 . 3 .

 43 Figure 4.3. Trends to equilibrium. Time evolution of the entropy (left) and {v y = 0}-profile of the numerical equilibria (right) for the 2D Fermi-Dirac case with initial datum f 0,1 and = r * with r ∈ {0.1, 0.5, 0.8} (top to bottom).

4. 2 . 2 .

 22 3D Fermi-Dirac relaxation. Let us now extend the analysis of the previous paragraph to the case of 3D Fermi-Dirac particles. Due to the large time cost required for computing the term Q R

Figure 4 . 4 .

 44 Figure 4.4. Trends to equilibrium. Time evolution of the entropy (left) and {v y = 0}-profile of the numerical equilibria (right) for the 2D Fermi-Dirac case with initial datum f 0,1 and = r * with r =∈ {0.9, 0.95, 0.99} (top to bottom).

Figure 4 . 5 . 0 (2π σ 0 ) 3 2 ≈ 6 .

 45026 Figure 4.5. Trends to equilibrium. Time evolution of max f h (t, •) (left) and {v y = 0}-profile of the numerical blow up profile (right) for the 2D Fermi-Dirac case with initial datum f 0,1 and = r * with r ∈ {1, 1.01, 1.05} (top to bottom).

  -4.10).

Figure 4 . 6 .

 46 Figure 4.6. Trends to equilibrium. Time evolution of the entropy (left) and {v y = 0}-profile of the numerical equilibria (right) for the 2D Fermi-Dirac case with initial datum f 0,2 and = r * with r ∈ {0.1, 0.5, 0.8} (top to bottom).

Figure 4 . 7 .

 47 Figure 4.7. Trends to equilibrium. Time evolution of f h (t, •) -f ∞ h 1 (left) and {v y = 0}-profile of the numerical equilibria (right) for the 2D Fermi-Dirac case with initial datum f 0,2 and = r * with r ∈ {0.9, 0.95, 0.99} (top to bottom).

  and 4.15).

Figure 4 . 8 .

 48 Figure 4.8. Trends to equilibrium. Time evolution of max f h (t, •) (left) and {v y = 0}-profile of the numerical equilibrium (right) for the 2D Fermi-Dirac case with initial datum f 0,2 and = * .

  ).

Figure 4 . 9 .

 49 Figure 4.9. Trends to equilibrium. Time evolution of the entropy (left) and profile of f h (10, (•, 0, 0)) (right) for the 3D Fermi-Dirac case with initial datum f 0,3 and = r * with r ∈ {0.1, 0.5, 0.8} (top to bottom).

Figure 4 . 10 .

 410 Figure 4.10. Trends to equilibrium. Time evolution of the entropy (left) and profile of f h (10, (•, 0, 0)) (right) for the 3D Fermi-Dirac case with initial datum f 0,3 and = r * with r ∈ {0.9, 0.95} (top to bottom).

Figure 4 . 11 .

 411 Figure 4.11. Trends to equilibrium. Time evolution of the entropy (left) and profile of f h (9.69, (•, 0, 0)) (right) for the 3D Fermi-Dirac case with initial datum f 0,3 and = 0.99 * .

Figure 4 . 14 .P - 1 p 1 p 1 p

 414111 Figure 4.14. Trends to equilibrium. Time evolution of the entropy (left) and profile of f h (10, (•, 0, 0)) (right) for the 3D Fermi-Dirac case with initial datum f 0,4 and = r * with r ∈ {0.1, 0.2, 0.5, 0.8} (top to bottom).We remark that, for any n ∈ I N ,Q R 1,q (G h , G h , G h ) n =

Figure 4 . 15 .

 415 Figure 4.15. Trends to equilibrium. Time evolution of fh (t, •) -f ∞ h 1 (left)and profile of f h (10, (•, 0, 0)) (right) for the 3D Fermi-Dirac case with initial datum f 0,4 and = 0.9 * .

Figure 4 . 16 .,Figure 4 . 17 .P - 1 p

 4164171 Figure 4.16. Trends to equilibrium. Time evolution of max f h (t, •) (left) and {v y = v z = 0}-blowup profile of f h (right) for the 3D Fermi-Dirac case with initial datum f 0,4 and = r * with r ∈ {0.95, 0.99}(top to bottom).

1 p 1 p,Figure 4 . 18 .

 11418 Figure 4.18. Trends to equilibrium. Time evolution of the entropy (left) and profile of f h (10, (•, 0, 0)) (right) for the 3D Bose-Einstein case with initial datum f 0,3 and = r * with r ∈ {0.1, 0.5, 0.8} (top to bottom).

Figure 4 . 19 .

 419 Figure 4.19. Trends to equilibrium. Time evolution of the entropy (left) and {v y = v z = 0}-blowup profile of f h (right) for the 3D Bose-Einstein case with initial datum f 0,3 and = r * with r ∈ {0.9, 0.95} (top to bottom) .

Figure 4 . 20 .

 420 Figure 4.20. Trends to equilibrium. Time evolution of max f h (t, •) (left) and {v y = v z = 0}-blowup profile (right) for the 3D Bose-Einstein case with initial datum f 0,3 and = r * with r ∈ 1, 1.05 (top to bottom .

Table 1 .

 1 Spectral

	u b

accuracy. Values of the numerical residual

Table 2 .

 2 Spectral accuracy. Values of the numerical residual

	u b

Table 4 .

 4 u b with rescaling without rescaling with rescaling without rescaling Spectral accuracy. Values of the numerical residual

		16 2	8.75850e-02	4.77941e+00	8.75850e-02	2.09808e+00
	L = 4.0 σ = 1.0	32 2 64 2 128 2	4.34564e-03 6.97576e-06 8.18675e-10	1.34244e-01 7.27236e-05 8.90710e-06	4.34564e-03 3.31081e-06 1.14691e-09	6.84172e-02 1.13026e-04 2.36536e-05
		256 2	9.06688e-10	8.27488e-06	9.06688e-10	2.52684e-05
		16 2	8.00293e-02	3.30798e+00	9.34962e-02	5.15603e+00
	L = 6.0 σ = 1.0	32 2 64 2 128 2	3.39670e-03 8.44540e-06 2.13457e-11	2.43273e-01 2.84543e-04 5.45586e-09	1.14900e-03 4.22200e-06 1.90767e-11	3.48709e-01 5.26459e-04 2.75652e-08
		256 2	1.24522e-11	2.90246e-09	1.98063e-11	2.91932e-08
		16 2	8.59697e-02	8.09303e+00	8.59697e-02	7.84054e+00
	L = 8.0 σ = 1.0	32 2 64 2 128 2	3.20116e-03 1.89958e-05 1.25941e-09	3.22578e-01 3.67449e-03 8.52464e-08	1.01646e-02 3.76830e-05 1.60801e-09	7.22433e-01 2.61653e-03 2.33918e-07
		256 2	8.48209e-14	4.54954e-12	4.54874e-14	7.25384e-12

Table 8 .

 8 Spectral accuracy. Values of the numerical residual

	1.63056e-05	4.36373e-07	2.17488e-06

Note that, the relaxation of f is always entropic for Bose-Einstein particles but can be non-entropic for Fermi-Dirac particles

T.R. was partially funded by Labex CEMPI (ANR-11-LABX-0007-01) and ANR Project MoHyCon (ANR-17-CE40-0027-01). 1

classification of relaxation phenomena. This will be the main topic of the upcoming followup [START_REF] Mouton | On Deterministic Numerical Methods for the quantum Boltzmann-Nordheim Equation. II. Fast algorithms[END_REF] of our work.

Appendix A. Fast evaluation of the collision kernel

Let us briefly present how we practically compute the collision operator, by focusing on its most computationally expensive part Q R 1,q (G h , G h , G h ).