Numerical study of Bose-Einstein condensation in the Kaniadakis-Quarati model for bosons

Kaniadakis and Quarati (1994) proposed a Fokker-Planck equation with quadratic drift as a PDE model for the dynamics of Bose-Einstein condensation in a homogeneous Bose gas. It is an open question whether this equation has solutions exhibiting condensates in finite time. The main analytical challenge lies in the continuation of exploding solutions beyond their first blow-up time while having a linear diffusion term. We present a thoroughly validated time-implicit numerical scheme capable of simulating solutions for arbitrarily large time, and thus enabling a numerical study of the condensation process in the Kaniadakis-Quarati model. We show strong numerical evidence that above the critical mass rotationally symmetric solutions of the Kaniadakis-Quarati model in 3D form a condensate in finite time and converge in entropy to the unique minimiser of the natural entropy functional at an exponential rate. Moreover, we show that transient condensates may occur for sufficiently concentrated initial data.


Introduction
In this paper we consider the following family of nonlinear Fokker-Planck equations where γ > 0 is a fixed parameter and f = f (t, v) ≥ 0. We are particularly interested in the case γ = 1, the physically most relevant one, in which equation (1.1) is known as the Bose-Einstein-Fokker-Planck equation or Kaniadakis-Quarati model (KQ). It was introduced by Kaniadakis and Quarati [33] as a model for the dynamics of the velocity distribution of a spatially homogeneous Bose gas.
Physical background (γ = 1). The feature in which KQ differs from the linear Fokker-Planck equation consists in the additional factor (1 + f ) in the drift term. This factor, leading to a nonlinear equation, arises from the assumption of indistinguishability of identical quantum particles. Indeed, in contrast to classical mechanics, in a quantum system of identical and indistinguishable particles, the presence of particles in a given energy state influences the probability of further quantum particles joining that state.
Here we are interested in systems of bosons, whose wave function is symmetric. This results in an increase in the transition probability, which, in the continuum model, is encoded in the extra factor (1 + f ). For KQ the choice d = 3 is the physically most interesting space dimension. In this case the problem exhibits a finite critical mass m c above which condensates are expected to emerge in finite time, see below for more details. However, in the literature little is known about the possible formation of condensates in 3D KQ.
Variational structure and steady states. Equation (1.1) has a natural entropy functional, given by where Φ(f ) := 1 γ f 0 log s γ 1+s γ ds and thus Φ ′′ (f ) = 1/h(f ) for h(s) := s(1 + s γ ). Indeed, formally, equation (1.1) can be rewritten as where δH δf (f ) denotes the variational derivative of H at f . Thus, for any sufficiently regular, positive (and hence mass conserving) solution f = f (t, v) of eq. (1.1), one obtains the entropy dissipation identity Notice, however, that due to the presence of the (quantum correction) term s γ in the definition of h(s), equation (1.2) is not a gradient flow of the functional H with respect to the classical Wasserstein metric. The mobility h(s) associated to the nonlinear continuity equation (1.2) is convex leading to well-known issues of ill-defined Wasserstein-like metrics to render rigorous the gradient flow structure [24] in contrast to the Fermi-Dirac case [17,18]. However, this formal gradient flow structure will play an important role for deriving suitable numerical schemes to deal with condensation. We observe that, given a sufficiently regular positive function f , the RHS of equation (1.3) is strictly negative unless ∇ δH δf (f ) = 0. The regular solutions of this equation are henceforth referred to as the steady states associated with problem (1.1). They are explicitly given by Notice that f ∞,θ is smooth and integrable for θ > 0, and the family {f ∞,θ } is strictly ordered and approaches f c := f ∞,0 from below as θ ց 0. Furthermore, letting m c := f c , the map (0, ∞) ∋ θ → m θ := f ∞,θ ∈ (0, m c ) is a bijection, and m c < ∞ if and only if γ > 2 d , i.e. if and only if the problem is L 1 -supercritical. While f ∞,θ is the unique minimiser of H among non-negative integrable functions of mass m = m θ , for m > m c the problem of minimising H under mass constraint does not have a regular solution. Since Φ is sublinear at infinity, the natural extension of the entropy functional to the set of finite non-negative Borel measures M + b , again denoted by H, is given by where f denotes the density of the absolutely continuous part of µ. The extension is convex and lower-semicontinuous with respect to weak-star convergence in M [23,7]. In [7] it is shown via an explicit calculation that the extended functional has a unique minimiser among finite non-negative measures of mass m > m c , which is given by The above comments on the entropy functional and the steady states of equation (1.1) apply to the problem posed on the whole space R d (assuming sufficient decay as |v| → ∞) as well as to the problem on a centred ball B(0, R 1 ) subject to no-flux boundary conditions.
Dynamics of the Kaniadakis-Quarati model. As noted in [16], in the L 1 -subcritical case, d = 1, KQ is globally wellposed in the classical sense for sufficiently regular initial data, and solutions converge to equilibrium at an exponential rate [11]. In the L 1 -critical case, d = 2, solutions are also globally regular and converge to equilibrium [9]-with an exponential rate in the spatially isotropic case f (t, v) = g(t, |v|). The approach in [9] exploits the fact that the 2D KQ in isotropic coordinates can be transformed to a linear Fokker-Planck equation, which leads to explicit solutions also for the nonlinear equation. For 3D KQ Toscani [46] proved via contradiction the existence of solutions blowing up in finite time. Finite-time blow-up in this reference is obtained for any solution of sufficiently large mass m (above a technical threshold far larger than the critical mass), but also for solutions of arbitrarily small mass provided they are initially sufficiently concentrated near the origin.
L 1 -supercritical Fokker-Planck model for bosons in 1D. The one-dimensional case of equation (1.1) with γ > 2 was recently studied in the ref. [16], both on the entire line as well as on a centred interval subject to zero-flux boundary conditions. We would like to point out that the successful numerical experiments reported below, which are based on the equation for the pseudo-inverse cumulative distribution function u(t, ·) of f (t, ·) (cf. Section 2), triggered the rigorous analysis in [16], which is itself based on the equation in mass variables. In [16] the authors obtain global-in-time existence and uniqueness of solutions u in the viscosity sense for initial data corresponding to sufficiently regular positive densities f 0 . These solutions are non-decreasing in the mass variable, here denoted by x. It is further shown that such solutions u = u(t, x) are smooth away from {u = 0} and that, in the original variables, u(t, ·) corresponds to a finite measure µ(t) ∈ M + b of the form where the map t → x p (t) := L 1 ({u(t, ·) = 0}) is continuous and the function f (t, ·) ∈ L 1 + is smooth away from the origin, where it satisfies equation (1.1) in the pointwise sense. Moreover, whenever the density f (t, ·) is unbounded at the origin, its spatial blow-up profile has the form This makes it possible to extend entropy methods globally in time and to deduce convergence to the measure of the same mass which minimises the entropy. In the case m > m c , the minimiser has a positive Dirac mass at the origin, and the solution must eventually have a non-trivial condensate component. On the other hand, if m < m c , the minimiser is smooth, and from the bound (1.5) it can easily be deduced that in this case there exists T ∈ (0, ∞) such that x p (t) = 0 for all t ≥ T (see [16,Corollary 4.5]). From this observation combined with an adaptation of the finite-time blow-up argument in [46] one infers the existence of solutions whose condensate component x p = x p (t) is not identically zero but compactly supported in (0, ∞) (see [30] for details). We refer to this phenomenon as a transient condensate. Below we will see that the L 1 -supercritical case in 1D of the family of nonlinear Fokker-Planck equations (1.1), corresponding to γ > 2, appears to be a good caricature for the dynamical behaviour of the physically interesting case of the 3D KQ model in radial coordinates. Let us finally mention that equation (1.1) in 1D and without the diffusion term was studied in [14] showing that condensates always form in finite time and that their mass is increasing in time so that, once formed, they never dissolve. The results reported here and in [16] show the genuine countereffect of linear diffusion on condensation. Indeed, the presence of diffusion leads to the possibility of non-monotonic behaviour of the condensate part x p (t) and the existence of transient condensates, as proved in one dimension for γ > 2 and conjectured in the three dimensional case for γ = 1.
Purpose and results. The main purpose of this work is to provide strong numerical evidence for the existence of solutions to the 3D KQ forming a Bose-Einstein condensate in finite time. Our numerical results suggest that any rotationally symmetric solution above the critical mass will eventually have a non-trivial condensed component. From our simulations a rather clear picture of the dynamical properties of KQ in 3D in the isotropic case will emerge: the long-time asymptotics will be identified, which the numerical solution converges to in entropy at an exponential rate. Numerical evidence is provided for the possibility of the condensed part failing to be monotonic in time and for even dissolving completely. Before investigating KQ in 3D, we will apply the numerical scheme to the caricature of the L 1 -supercritical case in 1D, i.e. (1.1) with γ > 2, in order to numerically reproduce the analytic results established rigorously in [16], see Section 3.1. Since no explicit (non-stationary) solutions are available in 1D, the 1D scheme (in the L 1 -supercritical case) will be validated by numerically analysing the convergence behaviour under mesh refinement with respect to a reference solution on a very fine mesh. Furthermore, our numerical experiments in 1D will replicate the existence of transient condensates and provide strong evidence that the decay of the entropy is exponential. Concerning the scheme for rotationally symmetric solutions of KQ we perform a careful validation in 2D, where explicit solutions are available.
Numerical scheme. The proposed numerical scheme is based on the variational formulation of equation (1.1) using a mass transportation Lagrangian approach. It corresponds to writing the formal gradient flow with respect to a Wasserstein-like distance in terms of the inverse of the cumulative distribution functions as in [8,20]. We use this variational structure in terms of the inverse cumulative distribution function as a basis of the proposed discretisation. Notice that in 1D the classical Wasserstein distance corresponds to the L 2 distance between the inverses of the distribution functions, and for radial functions in any dimension the classical Wasserstein distance between them can also be expressed in terms of the mass distributions. This variational formulation has the inherent advantage that conservation of mass, positivity and energy decay of solutions follow by construction. Moreover, in the one-dimensional and in the radial setting, the expression in terms of the inverses of distribution functions allows for the efficient computation of the otherwise computationally expensive Wasserstein-like distances. There has been an increased interest in related structure preserving Lagrangian schemes in the last years, see for example [29,8,19,40,22,20,13,15]. The numerical analysis of these schemes is still underdeveloped with partial results in [8,40,10,13,15]. Let us finally point out that free energy decaying numerical schemes in the original variables based on finite volume schemes have been proposed in [12,41,1,2] and references therein. These schemes fail to go beyond the blow-up time since they cannot resolve the presence of Dirac concentrations while accurately following the evolution of the smooth part of the solution.
Comparison with other models for Bose-Einstein condensation. There are many other models in the literature which have been suggested in the context of Bose-Einstein condensation. Of particular interest (due to similar phenomena) is a certain class of kinetic equations generally referred to as quantum Boltzmann equations, which have been derived following the usual derivation of the Boltzmann equation while taking into account the quantum effect by using Bose-Einstein statistics when counting the number of collisions between particles, see e.g. [21].
Significant progress has been made in the analysis of the spatially homogeneous and velocity isotropic Boltzmann-Nordheim equation, also called the Boltzmann equation for Bose-Einstein particles, see [35,36,26,25,3,37,38]. Loosely speaking (and omitting many details), the authors of the cited references are able to establish the existence of generalised solutions which form a Bose-Einstein condensate in finite time. It is not known whether the generalised solutions introduced in these references are unique. Let us remark that for γ = 1 the steady states (1.4) coincide with the classical Bose-Einstein distributions and the functional Φ(f ) dv agrees (up to a sign convention) with the entropy associated to the homogeneous Boltzmann-Nordheim equation for bosons, see [26,32]. In contrast to equation (1.1), the Boltzmann-Nordheim equation formally preserves the kinetic energy |v| 2 2 f dv. On the numerical side, schemes to approximate the Boltzmann-Nordheim equation or quantum Boltzmann equation for bosons have been devised and used to gain insights into these equations, see [39,6,31,28] and the references therein. The formal derivation of the 3D KQ model from the Boltzman-Nordheim equation has not been numerically tackled. Despite the design of numerical schemes for the Boltzmann-Nordheim equation, only few numerical studies attempt to go beyond the first blow-up time (where the velocity distribution ceases to be bounded). In [43,42,34,45] the authors observe that at the first blow-up time the solution has an integrable power law singularity near the lowest energy state and, in general, there will be a non-trivial flux of particles entering that state. The hypothesis of mass conservation then leads to a law for the time evolution of the condensate component, resulting in a coupled system. The methods do not appear to allow to track in a precise way the evolution after blow-up. Our approach is very different as it does not require any distinction between the times where the velocity distribution is bounded and the times where it is unbounded, and enables a detailed study of the dynamics of singular solutions. Let us finally mention that other descriptions have been used both analytically and numerically to study the behaviour beyond condensation in the quantum Boltzmann equation. In some of them the kinetic equation is coupled to a nonlinear Schrödinger equation (cubic, Gross-Pitaevski) modelling the evolution of the condensate, see [44,5,4,27] and the references therein for further details.
The results in the present paper suggest that in the isotropic case the dynamics of condensation in 3D KQ is in several aspects similar to the one of the Boltzmann-Nordheim equation as described rigorously in the references [26,25,3,37,38]. Common features comprise the finite-time blow-up for velocity distributions initially sufficiently concentrated near the origin, the finite-time condensation in the case where the equilibrium has a non-trivial singular part (in form of a Dirac measure at the origin), the spatial blow-up profiles following an integrable power law and the relaxation to equilibrium (under certain conditions). A feature which appears to be unknown for the quantum Boltzmann equations, but which will be studied extensively in this manuscript for 3D KQ, is the existence of transient condensates, i.e. the existence of (generalised) solutions, which are initially smooth, begin to form a condensate after some positive time, and will eventually be smooth again.
Plan of the manuscript. The remaining part of this manuscript is structured as follows: in Section 2 we discuss the numerical scheme for the 1D caricature of the 3D KQ given by the L 1 -supercritical 1D Fokker-Planck equation (1.1) with γ > 2 and its generalisation to the radial case in higher-dimensions with particular focus on the KQ model, γ = 1. Subsection 3.1 shows that the proposed numerical scheme does capture the main behaviour after blow-up in the 1D case: condensation, transient condensates for subcritical initial mass and convergence towards equilibirium. Subsection 3.2 validates the discretisation of the radial case by comparing to the explicit solutions given in [9]. Finally, Subsection 3.3 shows that the caricature given by the L 1 -supercritical Fokker-Planck equation (1.1) in 1D is essentially numerically correct for the 3D KQ model for radial initial data.

The transformed equation 2.1.1 One-dimensional case
Here, we consider the case d = 1 and assume that γ > 2, which represents the L 1supercritical regime. The equation satisfied by the pseudo-inverse u(t, ·) of the cumulative distribution function (cdf) and, upon multiplying by the factor (∂ x u) γ , it can be rewritten as While these new coordinates are generally known to be (numerically) favourable when investigating mass concentration phenomena in 1D, a particular feature of equation (2.2) is that the function u ≡ 0, which at the level of f corresponds to a Dirac delta at the origin, is an actual solution. Since mass conservation is a crucial feature of our Fokker-Planck model for bosons, the natural boundary condition for eq.
It enforces the flux of particles through the boundary to be zero. Formally, at the level of u, this means that the RHS of eq.
Hence, if the solution u is C 1,2 t,x near and up to the boundary, this becomes ∂ t u = 0 or, equivalently, This is the form we use in our numerical scheme. It corresponds to the Dirichlet conditions u(t, 0) = −R 1 , u(t, m) = R 1 for t > 0. As explained in Section 1, given a radius R 1 and a mass m of mass m which minimises the entropy H. At the level of u, we denote this minimiser by u ∞ , we further let H(u) := H(µ), where µ is the measure associated with the generalised inverse of u and will, in places, abbreviate H ∞ := H(u ∞ ) = H(µ ∞ ). The dependence of u ∞ on R 1 and m will be omitted.

Higher-dimensional case under radial symmetry
For isotropic solutions f (t, v) = g(t, |v|), v ∈ R d , we can perform a similar transformation in higher dimensions. In radial form, equation (2.2) reads As a first ansatz one might try to consider the equation for the (pseudo-) inverse R(t, z) of the radial cdfM (t, r) = r 0 g(t, s)s d−1 ds. However, for bounded densities f the function M is of class O(r d ) as r → 0, implying that R(t, ·) is at most 1/d-Hölder near z = 0 and ∂ z R z 1/d−1 → ∞ as z ց 0, whenever d > 1. We therefore consider the normalised version N(t, s) =M (t, s 1/d ) or, equivalently, which satisfies ∂ s N(t, s) = 1 d g(t, s 1/d ), and let S(t, ·) denote the pseudo-inverse of N(t, ·), so that S = R d . From the formal relation N(t, S(t, z)) = z we deduce (omitting the time argument) .
Then, the equation (2.3) for g leads to the following equation for S: Since we want our scheme to be able to deal with condensates, i.e. S(t, ·) ≡ 0 on some subinterval (0, z(t)), we multiply this equation by (∂ z S) γ to obtain Notice that if γ ∈ [1, 2), the viscosity term has a factor which becomes unbounded when S forms a condensate. We therefore consider for a small parameter 0 < ε ≪ 1 the following regularisation We are mostly interested in the KQ model (where γ = 1), and will thus focus on the equation where d = 2, 3. Notice that a positive ε decreases the strength of diffusion significantly when ∂ z S ε. In order to counterbalance this effect, which may potentially lead to numerical artefacts when investigating the expected phenomenon of condensation, we propose an artificial viscosity type regularisation of the form where 0 < δ ≪ 1 is a small parameter. Belowm (resp.m c ) denotes the total mass of the initial datum f 0 (resp. of f c ) on B(0, R 1 ) multiplied by the factor

The fully discrete scheme
The equations (2.2) and (2.5) are discretised fully implicitly in time. The resulting semidiscrete nonlinear system is then discretised using finite differences and solved by the Newton-Raphson method. The proposed schemes are based on earlier works presented in [8,20]. Let τ be the discrete time step and h the spatial discretisation. We denote by u n i = u(nτ, ih) with n ∈ N 0 and i ∈ {1, · · · N − 1} the discrete solution. Then the implicit time discretisation is given by The finite difference approximation in space is chosen in such a way as to preserve the equation's symmetry, viz.
for i = 1, . . . , N − 1, complemented with the boundary conditions u n 0 = u 0 0 = −R 1 and u n N = u 0 N = R 1 . We use a similar discretisation for (2.5), viz. Algorithm. Given u n−1 the discrete approximation u n at the subsequent time point is computed using a Newton-Raphson iteration. The iteration is stopped as soon as the smallness condition F NR (u n , u n−1 , h, τ ) l 2 < 10 −8 is satisfied, where F NR (u n , u n−1 , h, τ ) i is given by the LHS of equation (2.6) multiplied by h γ . For S we proceed similarly.
Remark 2.1. In the simulations exhibiting the numerically somewhat delicate condensation phenomenon, the inverse cdf becomes slightly non-monotonic during the Newton-Raphson iteration, which leads to very small imaginary parts in the above scheme and of the solution at the subsequent time step. In our actual code we therefore rearrange the approximation in each Newton-Raphson iteration to ensure monotonicity. Alternatively, one can replace the first derivatives u x by their absolute values |u x | and discretise and simulate this equation. In practice, the differences between the results using the first and the second option are negligible. A similar statement applies to the higher-dimensional case, where we choose again the option of the monotonic rearrangement.

Numerical experiments
In this section we describe the validation of our scheme, and present and discuss our numerical experiments.

L 1 -supercritical bosonic Fokker-Planck model in 1D: simulations replicating the theory
First, we demonstrate the reliability of the proposed numerical scheme in 1D by reproducing the features proved in [16]. As mentioned in the introduction, these simulations triggered the analytical study of the transformed equation (2.2) since our Newton-Raphson iterations were converging and leading to a clear extension of the solution after the blow-up time. In addition, we use the scheme to predict that the entropy decays at an exponential rate, even after the onset of a condensate.
If not stated otherwise, we choose γ = 2.9 and use a centred Gaussian as initial datum, viz.
for fixed positive constants A and σ. Moreover, we always set R 1 = 1. We remark that for d = 1 and the above choice of γ and R 1 the critical mass m c takes the numerical value m c ≈ 5.37.

Validation in 1D
In the first numerical experiment to validate the 1D scheme (2.6) we compare the solution for a given mesh with a numerical reference solution calculated on a fixed and much finer mesh. We set σ = 0.7, A = 4.5 in (3.1) as well as T = 0.025. For simplicity, the mass variable x ∈ [0, m] is often referred to as the spatial variable. The numerical reference solution is computed on a grid of 12801 (equidistant) spatial mesh points and a total number of 1000 (equidistant) time points. Notice that the values of the parameters A and σ coincide with those in (P1) below and observe that, in the simulations based on (P1), well before the final time T = 0.025 chosen for our validation, a significant amount of mass has accumulated at the origin (cf. Figures 1a and 1c). Therefore, our validation covers the case in which condensation occurs.     Table 2 indicates the L 2 space-time error between computed and reference solution. The results suggest a second order dependence of the error on the spatial increment and a first order dependence on the temporal increment. As long as the solution is not degenerate, this can be explained by the fact that we use an implicit Euler scheme in time (which is first-order accurate), a central finite difference discretisation in space (whose truncation error is of second order) and have chosen a high resolution in time for the test using purely spatial refinement, which makes the temporal error negligible in this test. Notice, however, that the degenerate case requires more care and that, in this work, we do not provide a rigorous numerical analysis of the scheme.

Comparing simulations and theoretical results
In order to numerically confirm the dynamical properties of eq. (1.1) in 1D established in [16], we run our scheme with the following four sets of parameters covering the masssuper resp. -subcritical, the asymmetric case as well as the case of the initial datum being highly concentrated near the origin v = 0:  1.08  The approximate total mass for each of these simulations is indicated in part (a) of the corresponding figure: it is the maximal value of the part of the horizontal axis which is displayed.
The convergence to the minimiser of the entropy can be clearly observed in Figures 1a  and 2a. Beyond, Figures 1b, 2b, 3c and 3d, which show the evolution of the relative entropy H(u(t, ·)) − H(u ∞ ), suggest an exponential decay of the entropy. Notice that in the L 1 -supercritical case no theoretical results have been established regarding the decay rate of the entropy.
The finite-time condensation in the mass-supercritical case is well confirmed by the simulations (P1)&(P4). Recall that at the level of u the size of the condensate corresponds to the measure L 1 ({u(t, ·) = 0}) of the flat part of u(t, ·) at height zero, which we numerically determine by the criterion |u(t, ·)| < 10 −6 . Figure 1c shows the time evolution of the condensed part relative to the (conserved) total mass. It clearly shows the onset of a condensate after some time 0 < t ≪ 0.025. Further figures illustrating the formation of condensates are Fig. 1a, 2a and 2c. Interestingly, in Figure 2c the fraction of mass in the condensate is not monotonic, implying that, even when above the critical mass, a previously formed condensate may partially dissolve. Figures 1d and 2d show the behaviour of f (t, v)/f c (v) for 0 < v ≪ R 1 at the times t = 0.04 and t = 0.1. Notice that in both figures the asymptote agrees well with the predicted one, even though the solution u(t, ·) is not uniformly close to u ∞ .    Transient condensates. In Figure 3 the behaviour of a mass-subcritical, but initially very concentrated solution is compared to the solution emanating from a more spread out datum. In both cases the entropy decays exponentially. Observe that in the case of high concentration, the solution forms a condensate in finite time which eventually vanishes again. We refer to this phenomenon as a transient condensate. Recall that for d = 1 and γ > 2 the existence of transient condensates below the critical mass is known rigorously [30]. The simulations based on (P3) illustrate very explicitly how, after some finite time, the function u(t, ·) begins to forms a flat part at the horizontal axis, which eventually disappears again as the solution converges to the smooth, non-degenerate equilibrium (cf. Figure 3e).

Validation in 2D
In the case d = 2, KQ is L 1 -critical and-as shown in [9]-its isotropic form can be transformed in an explicit way to a linear Fokker-Planck equation, whose solutions are explicit by means of the fundamental solution for this problem in R 2 . Here we want to use these explicit solutions to validate the proposed numerical scheme. Since all simulations are performed on a finite domain with zero flux boundary condition, the solutions to KQ obtained upon this transformation are only approximations of the exact solutions to our problem. However, we obtain a good approximation of the solutions in B(0, R 1 ) ⊂ R 2 with zero flux provided R 1 is chosen sufficiently large. This is due to the fact that the exact solutions in R 2 emanating from the chosen initial data (Gaussians) have exponential decay in |v|. The same is true for their derivative with respect to v, implying that on the boundary ∂B(0, R 1 ) of a centred ball of large enough radius R 1 ≫ 1 the flux is negligible. Hence, the exact solutions on R 2 restricted to B(0, R 1 ) are close to the exact solutions on B(0, R 1 ) with zero flux. Next, we recall the transformation leading to the explicit formula of solutions on the whole space, as observed in [9]: the solutions of the linear Fokker-Planck equation 2) h(0, ·) = h 0 are given by means of the fundamental solution where a(t) = e −2t , b(t) = e 2t − 1, and K b (z) = (2πb) −1 e −|z| 2 /2b . More precisely, (for sufficiently regular data h 0 ) the solution of equation (3.2) takes the form The relation between non-negative, isotropic solutions f of 2D KQ and non-negative, isotropic solutions h of eq. (3.2) is given by      , and from formula (3.3) and relation (3.4) we infer an expression for the solution f , which shows, in particular, that f (T, ·) has exponential decay for any positive time T . In our actual code, we use the inverse cdf of f (T, ·).
Details on the tests. We choose R 1 to be the smallest radius satisfying f c (v) ≤ 10 −4 for |v| ≥ R 1 . This guarantees that for any not too large σ > 0, the function f (t, ·) is small outside B(0, R 1 ). Two different tests are performed using the following common set of parameters: A = 4, σ = 0.9, final time T = 0.04 and size of the coarsest mesh equal to n 0 = 25. Since the solution to the exact problem remains bounded, the tests are performed with ε = δ = 0.
In the first test the dependence of the L 2 distance at time T between exact and computed solution for different spatial resolutions is analysed. More precisely, for j = 0, . . . , N = 5 we compute the error where J j denotes the discrete mesh using a total number of 2 j n 0 + 1 mesh points intersected with the interval [0, m/2], S (j) exact denotes the exact solution restricted to the spatial mesh J j and S (j) the discrete solution computed on the mesh J j using a total number of 400 time steps. Since we expect a polynomial dependence of the error on the spatial increment, we then let rate(j) = log 2 (E j /E j+1 ). The results of the test can be found in Table 3. Theoretically, since in the present case of two space dimensions the original density f remains uniformly bounded in time, which implies that ∂ z S stays away from zero, the spatial discretisation based on central differences should guarantee a quadratic dependence of the truncation error on the spatial increment. The rates displayed in Table 3 are somewhat worse, possibly due to the fact that the mesh size has not been chosen sufficiently large to capture the asymptotic behaviour well enough.
In the second test we analyse the dependence of the L 2 space-time distance between exact and computed solution on the number of spatial and temporal grid points. The procedure is analogous to the first test except that the j-th mesh is obtained by using 2 j n 0 + 1 spatial and 2 j m 0 temporal grid points, where m 0 = 4, and that now the error is given by exact l 2 (I j ×J j ) · 2 −2j , where I j denotes the discrete temporal mesh consisting of 2 j m 0 time points. The results are displayed in Table 4 and suggest a linear rate of convergence. This is in line with the backward Euler scheme used for the time stepping.    For completeness, we also tested the dependence of the computed solution on the regularisation parameters ε and δ, even though this is not necessary for 2D KQ since the density is theoretically known to remain bounded. We obtained a polynomial decrease of the error.
Remark 3.2. The choice of the comparatively fine mesh in (P6) was made in order to obtain a sufficiently good approximation of the evolution of the entropy. See Fig. 5b, which suggests an exponential decay.

Large-time behaviour
Our simulations suggest that 3D KQ has properties which are very similar to the Fokker-Planck model for bosons in 1D in the L 1 -supercritical regime. Figures Fig. 4a, Fig. 5a and Fig. 6a suggest that in the long-time limit the numerical solution S(t, ·) converges to the minimiser of the entropy (at the level of S), which we here denote by S ∞ -omitting its dependence on the given mass m and the radius R 1 . Next, the decay of the relative entropy appears to be exponential in all three cases, see Figures Fig. 4b, Fig. 5b and Fig. 6c

Condensation
In both the mass-supercritical case (P6) and the case of high concentration near the origin (P7) we observe the onset of a flat part at the level of S(t, ·) at height zero after some finite time, see Fig. 5c and 6d. In the original variables this means that mass is gradually absorbed by the origin. Furthermore, Fig. 6d shows that, similarly to the observations in 1D (see Section 3.1), it is possible for mass previously concentrated at velocity zero to escape. In fact, the condensate component may even dissolve completely.
Thus, the fraction of particles in the condensate is, in general, not monotonic in time for the 3D Kaniadakis-Quarati model. At times where the solution has a non-trivial condensate component, we were interested in the spatial behaviour of S(t, ·) close to {S(t, ·) = 0}. Owing to the results on the 1D model, one would expect the function f (t, v)/f c (v) to behave like 1 + O(|v|) as |v| → 0, where f is related to S via formula (2.4). Figure 5d confirms this behaviour at a time, where the solution is still well away from S ∞ . Remark 3.3. In order to produce the transient condensate in Figure 6 it was necessary to choose the parameter δ strictly positive. The same simulation for δ = 0 results in the flat part being trapped at zero once it has formed. As explained in Section 2.1.2 and also in view of our results for the 1D model, this "stickiness" appears to be a numerical artefact resulting from the circumstance that a regularisation based on a positive ε but vanishing δ is imbalanced and favours condensation.

Conclusion
In this work we have proposed a numerical scheme for a class of bosonic Fokker-Planck equations able for the first time to cope with Dirac delta concentrations of (partial) mass at the origin in finite time and to go beyond this blow-up time. This is achieved by changing to mass variables and scaling suitably the equation to obtain an alternative formulation admitting a Dirac delta concentration at the origin as a possible steady state. These new PDEs are solved by implicit schemes, and their approximations by Newton-Raphson type methods are shown to be numerically convergent by mesh refinement, even beyond the blow-up time. We have illustrated different phenomena appearing in the radial 3D KQ model which mimic the phenomena observed and partially proved for the 1D caricature of 3D KQ, see [16]. We have provided strong numerical evidence for the existence of transient condensates for mass-subcritical data and for the convergence to entropy minimising equilibria consisting of a regular part, which is smooth away from the origin, and a condensed part. We have also analysed the numerical convergence rates to the asymptotic states.