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

TL;DR

This paper numerically investigates the Kaniadakis-Quarati model for bosons, demonstrating finite-time condensate formation, convergence to equilibrium, and universal blow-up profiles, thus providing insights into the PDE's solution behavior.

Contribution

It introduces a validated numerical scheme to simulate the model's solutions over large times, revealing finite-time condensation and universal blow-up characteristics.

Findings

Condensates form in finite time for supercritical mass

Solutions converge exponentially to entropy minimizers

Universal power-law blow-up profile near the origin

Abstract

Kaniadakis and Quarati (1994) proposed a Fokker--Planck equation with quadratic drift as a PDE model for the dynamics of bosons in the spatially homogeneous setting. 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. Our simulations further indicate that the spatial blow-up profile near the origin follows a universal power law and that transient condensates can occur for sufficiently concentrated initial data.

Tables

Table 1. Table 1 : Convergence to reference solution at time T = 0.025 𝑇 0.025 T=0.025 .

timesteps	meshsize	$L_x^2$ error	rate
1000	50	7.3825e-3	-
1000	100	2.1290e-3	1.7939
1000	200	5.6056e-4	1.9253
1000	400	1.4222e-4	1.9788
1000	800	3.5598e-5	1.9982
1000	1600	8.8061e-6	2.0152
1000	3200	2.0991e-6	2.0687

Table 2. Table 3 : Convergence to reference solutions using CN and (P3) .

timesteps	meshsize	$L_{t,x}^2$ error	rate
10	50	5.2392e-3	-
20	100	1.1085 e-3	2.2408
40	200	2.4257 e-4	2.1921
80	400	5.6873e-05	2.0926
160	800	1.3983e-05	2.0241

Table 3. Table 4 : Convergence to exact solution at the final time T = 0.04 𝑇 0.04 T=0.04 .

number of	mesh size	$L^2$ error	rate
time points		(at time $T$)
4000	25	6.2783e-3	-
4000	50	2.2323e-3	1.4919
4000	100	7.9661e-4	1.4866
4000	200	2.6080e-4	1.6109
4000	400	7.7921e-5	1.7428
4000	800	1.9283e-5	2.0147

Equations

$${\partial }_{t}f$$

$$f(0,\cdot )$$

$$\mathcal{H}(f):=\int \left(\frac{|v{|}^2}{2}f+\Phi (f)\right)\mathrm{d}v,$$

$${\partial }_{t}f=\nabla \cdot \left(h(f)\nabla \frac{\delta \mathcal{H}}{\delta f}(f)\right),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{H}(f)=-\int h(f){\left|\nabla \frac{\delta \mathcal{H}}{\delta f}(f)\right|}^2\mathrm{d}v.$$

$$f_{\infty ,\theta }(v)={\left({\mathrm{e}}^{\gamma (\frac{|v{|}^2}{2}+\theta )}-1\right)}^{-1/\gamma },\theta \ge 0.$$

$$\tilde{\mathcal{H}}:\mu \mapsto \int \left(\frac{|v{|}^2}{2}\mu (\mathrm{d}v)+\Phi (f)\mathrm{d}v\right),$$

$$f_c\cdot {\mathcal{L}}^d+(m-m_c){\delta }_0.$$

$$u(x)=\inf \left\{r:{\int }_{\{r^{\prime }\le r\}}f(r^{\prime })\mathrm{d}{r}^{\prime }\ge x\right\},x\in (0,\Vert f{\Vert }_{L^1}),$$

$$\mu (t)=f(t,\cdot )\cdot {\mathcal{L}}^1+x_p(t){\delta }_0,$$

$$f(t,v)=f_c(v)\cdot (1+O(|v|)=c_{\gamma }|v{|}^{-\frac{2}{\gamma }}\left(1+O(|v|)\right)\text{as }|v|\to 0,c_{\gamma }={\left(2/\gamma \right)}^{\frac{1}{\gamma }}.$$

$$M(t,v)={\int }_{-R_1}^{v}f(t,w)\mathrm{d}w,v\in [-R_1,R_1],$$

$${\partial }_{t}u=({\partial }_{x}u{)}^{-2}{\partial }_x^{2}u-u(1+(u_x{)}^{-\gamma }),$$

$${\partial }_{x}u=\frac{1}{f(u)}.$$

$$({\partial }_{x}u{)}^{\gamma }{\partial }_{t}u-\frac{1}{\gamma -1}{\partial }_x\left(({\partial }_{x}u{)}^{\gamma -1}\right)+u(({\partial }_{x}u{)}^{\gamma }+1)=0.$$

$$u=u_0\text{on }(0,\infty )\times \{0,m\}.$$

$$H(u)={\int }_{(0,m)}\left(\frac{|u{|}^2}{2}+\Psi (u_x)\right)\mathrm{d}x,$$

$$\Psi (s):=s\Phi (1/s)\text{ is convex with }{\Psi }^{\prime \prime }(s)=s^{-3}{\Phi }^{\prime \prime }(1/s)=\frac{1}{s^{3}h(1/s)}.$$

$${\partial }_{t}g=r^{1-d}{\partial }_r\left(r^{d-1}{\partial }_{r}g+r^{d}g(1+g^{\gamma })\right),t,r>0.$$

$$N(t,s)=\frac{1}{d}{\int }_0^{s}g(t,{\sigma }^{1/d})\mathrm{d}\sigma ,$$

$${\partial }_{z}S=\frac{d}{g(R)}.$$

$$\frac{1}{d}{\partial }_{t}S-d\frac{S^{2-2/d}}{({\partial }_{z}S{)}^2}{\partial }_z^{2}S+S(1+d^{\gamma }({\partial }_{z}S{)}^{-\gamma })=0.$$

$$({\partial }_{z}S{)}^{\gamma }\frac{1}{d}{\partial }_{t}S-d\cdot S^{2-2/d}({\partial }_{z}S{)}^{\gamma -2}{\partial }_z^{2}S+S(({\partial }_{z}S{)}^{\gamma }+d^{\gamma })=0.$$

$$({\partial }_{z}S{)}^{\gamma }\frac{1}{d}{\partial }_{t}S-d\cdot S^{2-2/d}({\partial }_{z}S+\epsilon {)}^{\gamma -2}{\partial }_z^{2}S+S(({\partial }_{z}S{)}^{\gamma }+d^{\gamma })=0$$

$$\{\begin{array}{ll}({\partial }_{z}S{)}^{\gamma }\frac{1}{d}{\partial }_{t}S-\frac{d}{\gamma -1}\cdot S^{2-2/d}\frac{\mathrm{d}}{\mathrm{d}z}({\partial }_{z}S+\epsilon {)}^{\gamma -1}+S(({\partial }_{z}S{)}^{\gamma }+d^{\gamma })=0,& \text{ if }\gamma >1,\\ ({\partial }_{z}S{)}^{\gamma }\frac{1}{d}{\partial }_{t}S-d\cdot S^{2-2/d}\frac{\mathrm{d}}{\mathrm{d}z}\log ({\partial }_{z}S+\epsilon )+S(({\partial }_{z}S{)}^{\gamma }+d^{\gamma })=0,& \text{ if }\gamma =1.\end{array}$$

$$d^{-1}{\partial }_{z}S{\partial }_{t}S-d{S}^{2-2/d}\frac{\mathrm{d}}{\mathrm{d}z}\log ({\partial }_{z}S+\epsilon )+S({\partial }_{z}S+d)=0,$$

$$d^{-1}{\partial }_{z}S{\partial }_{t}S-d(S+\delta {)}^{2-2/d}\frac{\mathrm{d}}{\mathrm{d}z}\log ({\partial }_{z}S+\epsilon )+S({\partial }_{z}S+d)=0,$$

$$\det \nabla u(t,x)\cdot f(t,u)=1.$$

$$\left[(\det \nabla u{)}^2{\Psi }^{\prime \prime }(\det \nabla u)\right]{\partial }_t{u}^i-{\partial }_{x_k}\left({\Psi }^{\prime }(\det \nabla u)(\mathrm{cof}(\nabla u){)}_k^i\right)+u^i=0$$

$$H_{\mathrm{ani},d}(u)={\int }_U\left(\frac{1}{2}|u{|}^2+\Psi (\det \nabla u)\right)\mathrm{d}x.$$

$${\left({\partial }_{x}u\right)}^{\gamma }\frac{u-u^n}{\tau }-\frac{1}{\gamma -1}{\partial }_x\left(({\partial }_{x}u{)}^{\gamma -1}\right)+u(({\partial }_{x}u{)}^{\gamma }+1)=0$$

Full text

Keywords: nonlinear Fokker–Planck equation; Bose–Einstein condensation; entropy decay; implicit-in-time discretization; Lagrangian scheme.

1 Introduction

In this paper we consider the following family of nonlinear Fokker–Planck equations

[TABLE]

where $\gamma>0$ is a fixed parameter and $f=f(t,v)\geq 0$. We are particularly interested in the case $\gamma=1$, in which equation (1.1) is known as the Kaniadakis–Quarati model for bosons (KQ). It was introduced by Kaniadakis and Quarati [36] as a model for quantum particles following Bose–Einstein statistics, obtained by adapting accordingly the transition probability rates in the kinetic model.

Physical background ($\gamma=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 with respect to permutations of particles. This results in an increase in the transition probability, which is encoded, in the continuum model, 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

[TABLE]

where $\Phi(f):=\frac{1}{\gamma}\int_{0}^{f}\log\left(\frac{s^{\gamma}}{1+s^{\gamma}}\right)\mathrm{d}s$ and thus $\Phi^{\prime\prime}(f)=1/h(f)$ for $h(s):=s(1+s^{\gamma})$. Indeed, formally, equation (1.1) can be rewritten as

[TABLE]

where $\frac{\delta\mathcal{H}}{\delta f}$ denotes the variational derivative of $\mathcal{H}$. 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

[TABLE]

Notice, however, that due to the presence of the (quantum correction) term $s^{\gamma}$ in the definition of $h(s)$, equation (1.2) is not a gradient flow of the functional $\mathcal{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].

We observe that, given a sufficiently regular positive function $f$, the RHS of equation (1.3) is strictly negative unless $\nabla\frac{\delta\mathcal{H}}{\delta 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

[TABLE]

Notice that $f_{\infty,\theta}$ is smooth and integrable for $\theta>0$, and the family $\{f_{\infty,\theta}\}$ is strictly ordered and approaches $f_{c}:=f_{\infty,0}$ from below as $\theta\searrow 0$. Furthermore, letting $m_{c}:=\int f_{c}$, the map $(0,\infty)\ni\theta\mapsto m_{\theta}:=\int f_{\infty,\theta}\in(0,m_{c})$ is a bijection, and $m_{c}<\infty$ if and only if $\gamma>\frac{2}{d},$ i.e. if and only if the problem is $L^{1}$-supercritical. While $f_{\infty,\theta}$ is the unique minimiser of $\mathcal{H}$ among non-negative integrable functions of mass $m=m_{\theta}$, for $m>m_{c}$ the problem of minimising $\mathcal{H}$ under mass constraint does not have a regular solution. Since $\Phi$ is sublinear at infinity, the natural extension $\mathcal{\widetilde{H}}$ of the entropy functional to the set of finite non-negative Borel measures $\mathcal{M}^{+}_{b}$ is given by

[TABLE]

where $f$ denotes the density of the absolutely continuous part of $\mu$. The extension is convex and lower-semicontinuous with respect to weak-star convergence in $\mathcal{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

[TABLE]

The above comments on the entropy functional and the steady states of equation (1.1) apply to the problem posed on the whole space $\mathbb{R}^{d}$ (assuming sufficient decay as $|v|\to\infty$) 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 [15], 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 [21]. 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 $2$D KQ in isotropic coordinates can be transformed to a linear Fokker–Planck equation, which leads to explicit solutions also for the nonlinear equation. For $3$D KQ Toscani [49] 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. Formal results on the dynamics of isotropic solutions to 3D KQ based on matched asymptotic expansions have been obtained in [47]. Our numerical simulations will qualitatively confirm some of the main findings in [47], which suggests that the dynamics depicted in this reference give a good hint at the typical behaviour of solutions. Our numerical experiments will, however, also indicate that the dynamics may, in general, display a richer variety of phenomena. The formal considerations in [47] rely on the assumption of a sufficiently spread out initial datum $f_{0}$. We would also like to emphasize that in contrast to [47] our scheme allows for a numerical study beyond the first blow-up time.

$L^{1}$-supercritical Fokker–Planck model for bosons in 1D.

The one-dimensional case of equation (1.1) with $\gamma>2$ was recently studied in the ref. [15], 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 in the present manuscript, which are based on the equation for the pseudo-inverse cumulative distribution function

[TABLE]

of the original density $f=f(t,\cdot)$ (cf. equation (2.1) in Section 2 below) triggered the rigorous analysis in [15], which is itself based on this reformulation. Let us briefly review those results of [15] which are of relevance for the present paper: 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}$ of finite mass $m$. 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 the push-forward measure $u(t,\cdot)_{\#}\mathcal{L}^{1}_{|(0,m)}=:\mu(t)\in\mathcal{M}^{+}_{b}$, generalising the problem in the original variables, has the the form

[TABLE]

where the map $t\mapsto x_{p}(t):=\mathcal{L}^{1}(\{u(t,\cdot)=0\})$ is continuous and the function $f(t,\cdot)\in L^{1}_{+}$ is smooth away from the origin, where it satisfies equation (1.1) in the pointwise sense. Moreover, whenever the density $f(t,\cdot)$ is unbounded at the origin, its spatial blow-up profile has the form

[TABLE]

See [33] for a refinement of (1.6). The above framework 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.6) it can easily be deduced that in this case there exists $T\in(0,\infty)$ such that $x_{p}(t)=0$ for all $t\geq T$ (see [15, Cor. 3.16]). From this observation combined with an adaptation of the finite-time blow-up argument in [49] one infers the existence of solutions whose condensate component $x_{p}=x_{p}(t)$ is not identically zero but compactly supported in $(0,\infty)$ (see [15, Cor. 3.18]). 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 $\gamma>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 analysed in [13] 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 [15] show the genuine countereffect of linear diffusion on condensation leading to transient condensates and non-monotonic behaviour of the condensate part $x_{p}(t)$, proved in one dimension for $\gamma>2$ and conjectured in the three dimensional case for $\gamma=1$.

Main numerical findings.

The main purpose of this work is to provide strong numerical evidence for the existence of solutions to 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 condensate component. From our simulations a rather clear picture of the dynamical properties of KQ in $3$D 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 $3$D, we will apply the numerical scheme to the caricature of the $L^{1}$-supercritical case in 1D, i.e. (1.1) with $\gamma>2$, in order to numerically reproduce the analytical results established rigorously in [15], see Section 3.1. Since non-stationary explicit solutions are not available in $1$D, the $1$D 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. Concerning the scheme for rotationally symmetric solutions of KQ we perform a 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 is motivated by the approach in [8, 20], where the gradient flow with respect to the Wasserstein distance is expressed in terms of the inverse of the cumulative distribution functions. Inherent in this approach is the conservation of mass property, which follows by construction. We would like to emphasize that concerning the Kaniadakis–Quarati model in 3D studied in the present work, far less is known rigorously as compared to the equations simulated in [8, 20] (porous medium equation, critical Keller–Segel) which have been exhaustively studied in the literature. This is partially explained by the fact that the variational structure for this problem cannot directly be exploited by resorting to established tools from optimal transportation theory. In fact, the potential difficulty in our situation lies in the circumstance that we do not have the Wasserstein gradient flow structure in a rigorous sense. We will, however, see that this precise structure is not required and our proposed scheme will be shown to preserve in particular the entropy decay property (rigorously in 1D and 2D for the semidiscrete case, see Section 2.2). Our numerical scheme is able to go beyond the first blow-up time and allows for exploring the qualitative behaviour after blow-up: blow-up profile, transient condensates and entropy decay. These good numerical properties, consistent in 1D with the existing theory, reassure us in our numerical findings in Section 3.3 concerning the 3D isotropic case. The fact that our numerical experiments clearly support the conjecture that the qualitative behaviour of condensates in 1D proven in [15] is expected in the most realistic case of 3D radially symmetric initial data can be regarded as the main contribution of this paper.

There has been an increased interest in related structure preserving Lagrangian schemes in the last years, see for example [32, 8, 19, 42, 22, 20, 12, 14]. The numerical analysis of these schemes is still underdeveloped with partial results in [8, 42, 10, 12, 14]. Let us finally point out that free energy decaying numerical schemes in the original variables based on finite volume schemes have been proposed in [11, 43, 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, in contrast to classical Boltzmann equations, are derived using Bose–Einstein statistics. Let us note that for $\gamma=1$ the steady states (1.4) coincide with the classical Bose–Einstein distributions and the functional $\int\Phi(f)\,\mathrm{d}v$ agrees (up to a sign convention) with the entropy associated to the homogeneous Boltzmann–Nordheim equation for bosons, see [28, 35]. In contrast to equation (1.1), the Boltzmann–Nordheim equation formally preserves the kinetic energy $\int\frac{|v|^{2}}{2}f\,\mathrm{d}v$. In the last two decades, significant progress has been made in the analysis of the Boltzmann–Nordheim equation in the homogeneous and velocity isotropic case [25, 26, 38, 28, 27, 3, 39, 40]. To roughly summarise the main results, the authors of the cited references are able to establish the existence of generalised mass- and energy-conserving solutions, which form a Bose–Einstein condensate in finite time and converge, in some sense and under certain conditions, to the entropy minimiser in the large-time limit. The results in the present paper suggest that in the isotropic case the dynamics of condensation in 3D KQ is in some aspects similar to the one of the Boltzmann–Nordheim equation as described rigorously in the references [28, 27, 3, 39, 40]. We note that, regarding the nature of singularities, in the Boltzmann–Nordheim equation many questions are still open.

Numerical schemes to approximate the Boltzmann–Nordheim equation or quantum Boltzmann equation for bosons have also been devised and used to understand their qualitative properties, see [41, 6, 34, 31] and the references therein. However, only few numerical studies attempt to go beyond the first blow-up time (where the velocity distribution ceases to be bounded). In [45, 44, 37, 48] 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 distinguishing 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 [46, 5, 4, 29] and the references therein for further details.

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 $\gamma>2$ and its generalisation to the radial case in higher-dimensions with particular focus on the KQ model, $\gamma=1$. We also briefly discuss the anisotropic case. Section 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 equilibrium. Section 3.2 validates the discretisation of the radial case by comparing to the explicit solutions given in [9]. In Section 3.3 we present the simulations of 3D KQ, which allow us to conclude 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.

2 Numerical method

Since we want our scheme to be able to deal with Dirac masses at the origin, our simulations are not based on the formulation (1.1); instead we follow and generalise the ansatz in the ref. [15] considering the equation satisfied by the (pseudo-) inverse cumulative distribution function (cdf) of $f(t,\cdot)$. In higher dimensions $d>1$, assuming rotational symmetry, we will consider the inverse of the radial cdf of $f(t,\cdot)$ (i.e. of the partial mass function)—appropriately normalised. As in the first part of [15], we consider our equations posed on a bounded domain, more precisely on the centred ball $B(0,R_{1})$ of radius $R_{1}>0$ with zero-flux boundary conditions.

2.1 Change of variables

2.1.1 One-dimensional case

Here, we consider the case $d=1$ and assume that $\gamma>2$, which represents the $L^{1}$-supercritical regime. The total mass of the initial datum $f_{0}$ is denoted by $m$. Then, the equation satisfied by the pseudo-inverse $u(t,\cdot)$ of the cumulative distribution function (cdf)

[TABLE]

of $f(t,\cdot)$ (cf. (1.5)) formally states

[TABLE]

where $x\in(0,m)$ denotes the mass variable. Formally, $f$ is related to $u$ by the identity

[TABLE]

Upon multiplying eq. (2.1) by the factor $(\partial_{x}u)^{\gamma}$, it can be rewritten as

[TABLE]

While these new coordinates are generally known to be (numerically) favourable when investigating mass concentration phenomena in $1$D, a particular feature of equation (2.2) is that the function $u\equiv 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, the natural boundary condition for eq. (1.1) states $\partial_{r}f+rf(1+f^{\gamma})=0$ on $(0,\infty)\times\{-R_{1},R_{1}\}$. 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. (2.1) is zero on $(0,\infty)\times\{0,m\}$. Hence, if the solution $u$ is $C^{1,2}_{t,x}$ near and up to the boundary, this becomes $\partial_{t}u=0$ or, equivalently,

[TABLE]

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=\|f_{0}\|_{L^{1}(-R_{1},R_{1})}$ there exists a unique $\mu_{\infty}\in\mathcal{M}^{+}_{b}([-R_{1},R_{1}])$ of mass $m$ which minimises the entropy $\mathcal{\widetilde{H}}$. At the level of $u$, we denote this minimiser by $u_{\infty}$. We further let $H(u):=\mathcal{H}(f)$ resp. $\mathcal{\widetilde{H}}(\mu)$, where $\mu=u_{\#}\mathcal{L}^{1}$ is the push-forward measure of the Lebesgue measure on $[0,m]$ under the map $u$ and will, in places, abbreviate $H_{\infty}:=H(u_{\infty})=\mathcal{\widetilde{H}}(\mu_{\infty})$. The dependence of $u_{\infty}$ on $R_{1}$ and $m$ will be omitted. For later reference, let us observe that $H(u)$ is formally given by

[TABLE]

where the function

[TABLE]

2.1.2 Higher dimensions – isotropic case

For isotropic solutions $f(t,v)=g(t,|v|),$ $v\in\mathbb{R}^{d}$, we can perform a similar transformation in higher dimensions. In radial form, equation (2.2) reads

[TABLE]

As a first ansatz one might try to consider the equation for the (pseudo-) inverse $R(t,z)$ of the radial cdf $\bar{M}(t,r)=\int_{0}^{r}g(t,s)s^{d-1}\,\mathrm{d}s$. However, for bounded densities $f$ the function $\bar{M}$ is of class $O(r^{d})$ as $r\to 0$, implying that $R(t,\cdot)$ is at most $1/d$-Hölder near $z=0$ and $\partial_{z}R\gtrsim z^{1/d-1}\to\infty$ as $z\searrow 0$, whenever $d>1$. We therefore consider the normalised version $N(t,s)=\bar{M}(t,s^{1/d})$ or, equivalently,

[TABLE]

which satisfies $\partial_{s}N(t,s)=\frac{1}{d}g(t,s^{1/d})$, and let $S(t,\cdot)$ denote the pseudo-inverse of $N(t,\cdot)$, so that $S=R^{d}$. From the formal relation $N(t,S(t,z))=z$ we deduce (omitting the time argument)

[TABLE]

Then, the equation (2.5) for $g$ leads to the following equation for $S$:

[TABLE]

Since we want our scheme to be able to deal with condensates, i.e. $S(t,\cdot)\equiv 0$ on some subinterval $(0,z(t))$, we multiply this equation by $(\partial_{z}S)^{\gamma}$ to obtain

[TABLE]

Notice that if $\gamma\in[1,2)$, the viscosity term has a factor which becomes unbounded when $S$ forms a condensate. We therefore consider for a small parameter $0<\varepsilon\ll 1$ the following regularisation

[TABLE]

or, equivalently,

[TABLE]

We are mostly interested in the KQ model (where $\gamma=1$), and will thus focus on the equation

[TABLE]

where $d=2,3$. Notice that a positive $\varepsilon$ decreases the strength of diffusion significantly when $\partial_{z}S\lesssim\varepsilon$. 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

[TABLE]

where $0<\delta\ll 1$ is a small parameter. Below $\bar{m}$ (resp. $\bar{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 $\frac{1}{|\partial B(0,1)|}$. Then, as in the $1$D case, the appropriate boundary conditions for equation (2.8) are $S(t,0)=0$ and $S(t,\bar{m})=R_{1}^{d}$.

As in Section 2.1.1 we denote by $S_{\infty}=S_{\infty}(R_{1},\bar{m})$ the pseudo-inverse normalised radial cdf of the unique (isotropic) minimising measure in $\mathcal{M}^{+}_{b}(\overline{B}(0,R_{1}))$ corresponding to the choice $(R_{1},m)$ of parameters, and generally let $H_{d}(S):=\mathcal{\widetilde{H}}(\mu)$, where $\mu$ is the unique isotropic measure in $\mathcal{M}^{+}_{b}(\overline{B}(0,R_{1}))$ satisfying $\mu(\overline{B}(0,r))=\nu([0,r^{d}])\cdot|\partial B(0,1)|$ and $\nu$ denotes the measure associated with the generalised inverse of $S$. We also abbreviate $H_{\infty}:=H(S_{\infty})$ and $H(t):=H_{d}(S(t))$.

2.1.3 Higher dimensions – anisotropic case

Let us briefly discuss that one can perform a related change of variables in higher dimensions without radial symmetry. With this aim, one needs to consider vector-valued transformations $u(t,\cdot):U\to V$, $U,V\subset\mathbb{R}^{d}$, which are formally related to the original density $f$ via

[TABLE]

Similarly to [30, 19], one finds that the system governing the evolution of $u=(u^{1},\dots,u^{d})^{T}$ can formally be written as

[TABLE]

for $i=1,\dots,d$, where $\Psi$ is defined as in (2.4). The entropy $H_{\mathrm{ani},d}(u)$ in the new variables takes the form

[TABLE]

Observe that in the vectorial case $H_{\mathrm{ani},d}(u)$ is no longer convex but merely polyconvex in $\nabla u$. This route could potentially allow to numerically analyse concentrations without radial symmetry in higher dimensions, as it is the case in 2D for aggregation and Keller–Segel type problems close to the blow-up time [20]. Even if this method deserves further exploration, we focus here on the isotropic case to capture the direct generalisation of the 1D behaviour in the 3D realistic setting.

2.2 The semidiscrete scheme

The scalar equations (2.2) and (2.8) are discretised fully implicitly in time. We let $\tau$ be the discrete time step and denote by $\{u^{n}\}_{n\in\mathbb{N}}$ the time-discrete solution of the implicit Euler discretisation of equation (2.2). More precisely, given a non-decreasing function $u^{n}$ satisfying $u^{n}(0)=-R_{1}$ and $u^{n}(m)=R_{1}$, the problem for $u=u^{n+1}$ reads

[TABLE]

subject to the Dirichlet boundary conditions $u^{n+1}(0)=-R_{1},u^{n+1}(m)=R_{1}$.

Let us here make a short digression to explain the main difference and potential difficulty of the present problem with respect to the Wasserstein gradient flows treated in [8, 20]. Those works are based on the idea that the Wasserstein gradient flow of the entropy/free energy in the original variables is equivalent to an $L^{2}$ gradient flow for the problem in the $u$-variables. Loosely speaking, the semidiscrete $L^{2}$ gradient flow for $H(u)$ reads as follows: given $\tilde{u}^{n}$ formally define $\tilde{u}^{n+1}$ as a solution of the problem

[TABLE]

The associated Euler–Lagrange equations, $\frac{\tilde{u}-\tilde{u}^{n}}{\tau}=-\nabla_{L^{2}}H(\tilde{u}),$ read

[TABLE]

To compare this with our problem, we write eq. (2.9) in the more concise equivalent form

[TABLE]

which suggests that in some sense a gradient flow structure is kept. At least, as will be shown below, we keep an important property in the semidiscrete numerical scheme, namely the monotonicity of the entropy. Recall that in 1D the entropy $H(u)$ in the $u$-variables (see (2.3)) is convex in the classical sense, and it is well-known that the implicit Euler scheme applied to a gradient flow of a convex functional satisfies the semidiscrete entropy inequality $H(\tilde{u}^{n+1})\leq H(\tilde{u}^{n})$ for all $n$. In our situation, thanks to the convexity of the integrand of $H$, the entropy decay along the sequence $\{u^{n}\}$ can be recovered by a simple estimate:

[TABLE]

Here, we used the fact that in the above integration by parts the boundary terms vanish since, by construction, $u=u^{n}$ on $\partial(0,m)$. This shows the entropy decay property of the semidiscrete scheme (2.9): $H(u^{n+1})\leq H(u^{n})$ for all $n$. We note that similar properties with a similar strategy of proof are found for finite volume schemes of gradient flows [2].

Remark 2.1* (Higher dimensions, isotropic case).*

In higher dimensions the entropy $H_{d}(S)$, introduced in Section 2.1.2, takes the form (see also (2.6))

[TABLE]

where $\Psi_{d}(s)=\Psi(\tfrac{s}{d})$ is again convex. If $d=2$, thanks to convexity, the implicit Euler discretisation of eq. (2.7) can be shown to keep the entropy decay by arguing as in the 1D case. In higher dimensions, $d>2$, this argument breaks down due to the kinetic part of the entropy failing to be a convex function of $S$. Notice, however, that the convexity in the highest order term, $\partial_{z}S$, is maintained.

2.3 The fully discrete scheme

The semidiscrete nonlinear system (2.9) is discretised using finite differences and solved by the Newton–Raphson method. In the one dimensional case, the finite difference approximation in space is chosen in such a way as to preserve the equation’s symmetry, viz.

[TABLE]

for $i=1,\dots,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 full discretisation for (2.8), viz.

[TABLE]

for $i=1,\dots,N$, where the boundary conditions are given by $S_{0}^{n}=S_{0}^{0}=0$ and $S_{N}^{n}=S_{N}^{0}=R_{1}^{d}$.

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_{\mathrm{NR}}(u^{n},u^{n-1},h,\tau)\|_{l^{2}}<10^{-8}$ is satisfied, where $F_{\mathrm{NR}}(u^{n},u^{n-1},h,\tau)_{i}$ is given by the LHS of equation (2.11) multiplied by $h^{\gamma}$. For $S$ we proceed similarly.

Remark 2.2*.*

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.

3 Numerical experiments

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

3.1 $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 [15]. 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 $\gamma=2.9$ and use a centred Gaussian as initial datum, viz.

[TABLE]

for fixed positive constants $A$ and $\sigma$. Moreover, we always set $R_{1}=1$. We remark that for $d=1$ and the above choice of $\gamma$ and $R_{1}$ the critical mass $m_{c}$ takes the numerical value $m_{c}\approx 5.37$.

3.1.1 Validation in 1D

We begin with validating the $1$D scheme (2.11) by comparing the solution for a given mesh with a numerical reference solution calculated on a fixed and much finer mesh. We set $\sigma=0.7$, $A=4.5$ in (3.1) as well as $T=0.025$. For simplicity, the mass variable $x\in[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 $\sigma$ 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 1 and 1). Therefore, our validation covers the case in which condensation occurs.

Table 2 displays the discrete $L^{2}_{x}$ error of the solution on the coarser mesh with respect to the reference solution, evaluated at the final time $T$, while 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.

Remark 3.1*.*

Higher-order implicit time discretisations could be considered to obtain better accuracy. Table 3 displays the convergence rates upon refinement of the space-time mesh using a Crank–Nicolson-type (CN) time discretisation for (2.2) with parameters (P3) and clearly confirms the second order accuracy of CN. However, we would like to point out that the initial datum determined by (P3) is mass-subcritical, and therefore the 2nd order accuracy is obtained for smooth solutions. Our simulations beyond blow-up indicate that the Newton solver for the implicit Euler scheme has better stability properties to cope with condensates than the CN scheme though.

3.1.2 Comparing simulations and theoretical results

In order to numerically confirm the dynamical properties of eq. (1.1) in $1$D, we run our scheme with the following four sets of parameters covering the mass-super resp. -subcritical, the asymmetric case as well as the case of the initial datum being highly concentrated near the origin $v=0$:

1. (P1)

$m>m_{c}:$ $\sigma=0.7$, $A=4.5$, $T=0.4$, $\tau=0.001$, $n=2001$ ($n:=$ number of spatial grid points). 2. (P2)

Asymmetric & $m>m_{c}:$ translated Gaussian $f_{0}(v)=A\mathrm{e}^{-|v-v_{0}|^{2}/(2\sigma^{2})}+0.1$ chosen as initial datum using the parameters $v_{0}=-1$, $\sigma=0.7$ and $A=4.5$. Moreover, $T=0.4$, $\tau=0.001$, $n=2001$. The shift by $+0.1$ ensures that the cdf of $f_{0}$ is numerically still well invertible close to $v=R_{1}$. 3. (P3)

$m<m_{c}:$ $\sigma=0.7$, $A=1.5$, $T=0.4$, $\tau=0.001$, $n=2001$. 4. (P4)

Concentrated & $m<m_{c}:$ $\sigma=0.1$, $A=1.5$, $T=0.4$, $\tau=10^{-6}$, $n=10001$.

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.

Entropy decay.

The convergence to the minimiser of the entropy can be clearly observed in Figures 1 and 2. Beyond, Figures 1, 2, 3 and 3, which show the evolution of the relative entropy $H(u(t,\cdot))-H(u_{\infty})$, indicate an exponential decay of the entropy. The red slopes in Figures 1, 2, 3 and 3 indicate the approximate slopes of the graphs averaged over the intervals where they are plotted. The computed slopes imply quantitative decays rates for the entropy of the form $e^{-\alpha t}$ with the following numerical values for $\alpha$: $\alpha\approx 23.7$ for (P1), $\alpha\approx 23.8$ for (P3), $\alpha\approx 23.1$ for (P4), and $\alpha\approx 23.0$ for (P2).

Remark 3.2*.*

In the mass-subcritical case $m<m_{c}$ there exists $T=T(u_{0})<\infty$ such that the mapping $u(t,\cdot)$ has no critical point for $t>T$, so that the density $f(t,\cdot)$ of its inverse is smooth (see [15, Corollary 4.5]). In this case one can exploit the fact that the entropy functional of the bosonic Fokker–Planck equation in 1D coincides with that of a nonlinear diffusion equation with linear drift to which the theory developed in [16] applies in order to deduce exponential decay of the entropy with rate $\alpha=2$ for $t\geq T$, i.e.

[TABLE]

This idea was used before in [21] for 1D KQ. The rate of convergence in the general case is still open.

Finite-time condensation for $m>m_{c}$.

The finite-time condensation in the mass-supercritical case is well confirmed by the simulations (P1)&(P2). Recall that the condensate corresponds to the zero level set of $u(t,\cdot)$, which we numerically determine by the criterion $|u(t,\cdot)|<10^{-6}$. Figure 1 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\ll 0.025$. Further figures depicting the formation of condensates are Fig. 1, 2 and 2. Interestingly, in Figure 2 the fraction of mass in the condensate is not monotonic, illustrating that, even when above the critical mass, a previously formed condensate may partially dissolve.

Blow-up profile.

Figures 1 and 2 show the behaviour of $f(t,v)-f_{c}(v)$ for $0<v\ll R_{1}$ at the times $t=0.04$ and $t=0.1$. The figures indicate an error of the form

[TABLE]

for suitable constants $c_{+}(t),c_{-}(t)\in\mathbb{R}$, which, for asymmetric solutions, need not necessarily coincide. The asymptotic behaviour in (3.2) not only confirms the leading order spatial profile obtained rigorously in [15] (see (1.6)), but indicates that the error with respect to $f_{c}$ may typically be of first order in $|v|$ and thus smaller than the order $1-2/\gamma$ ensured by formula (1.6). (A rigorous derivation of the improved error control can be found in [33].) Let us also mention that in both figures the solution $u(t,\cdot)$ is not uniformly close to $u_{\infty}$, so that the asymptotic behaviour of the density near the origin at the chosen times is not due to the fact that the long-time limit of the density equals $f_{c}$.

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 $\gamma>2$ the existence of transient condensates is known rigorously [15]. The simulations based on (P4) illustrate very explicitly how, after some finite time, the function $u(t,\cdot)$ 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 3).

3.2 Validating KQ by means of explicit solutions 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 $\mathbb{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})\subset\mathbb{R}^{2}$ with zero flux provided $R_{1}$ is chosen sufficiently large. This is due to the fact that the exact solutions in $\mathbb{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 $\partial B(0,R_{1})$ of a centred ball of large enough radius $R_{1}\gg 1$ the flux is negligible. Hence, the exact solutions on $\mathbb{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

[TABLE]

are given by means of the fundamental solution

[TABLE]

where $a(t)=\mathrm{e}^{-2t},\;\;b(t)=\mathrm{e}^{2t}-1,\text{ and }K_{b}(z)=(2\pi b)^{-1}\mathrm{e}^{-|z|^{2}/2b}.$ More precisely, (for sufficiently regular data $h_{0}$) the solution of equation (3.3) takes the form

[TABLE]

The relation between non-negative, isotropic solutions $f$ of $2$D KQ and non-negative, isotropic solutions $h$ of eq. (3.3) is given by

[TABLE]

where

[TABLE]

We initialise our tests again with a centred Gaussian of the form

[TABLE]

for fixed positive constants $A$ and $\sigma$. Then the initial datum $h_{0}$ corresponding to $f_{0}$ via the transformation (3.5) is given by

[TABLE]

and from formula (3.4) and relation (3.5) we infer an expression for the solution $f$, which shows, in particular, that $f(T,\cdot)$ has exponential decay for any positive time $T$. In our actual code, we use the inverse cdf of $f(T,\cdot)$.

Details on the tests.

We choose $R_{1}$ to be the smallest radius satisfying $f_{c}(v)\leq 10^{-4}$ for $|v|\geq R_{1}$. This guarantees that for any not too large $\sigma>0$, the function $f(t,\cdot)$ is small outside $B(0,R_{1})$.

Two different tests are performed using the following common set of parameters: $A=4$, $\sigma=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 $\varepsilon=\delta=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,\dots,N=5$ we compute the error

[TABLE]

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)}_{\mathrm{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 $\mathrm{rate}(j)=\log_{2}(E_{j}/E_{j+1})$. The results of the test can be found in Table 5. Theoretically, since in the present case of two space dimensions the original density $f$ remains uniformly bounded in time, which implies that $\partial_{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 5 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

[TABLE]

where $I_{j}$ denotes the discrete temporal mesh consisting of $2^{j}m_{0}$ time points. The results are displayed in Table 5 and suggest a linear rate of convergence. This is in line with the backward Euler scheme used for the time stepping.

Remark 3.3* (Validation of regularisation).*

For completeness, we also tested the dependence of the computed solution on the regularisation parameters $\varepsilon$ and $\delta$, even though this is not necessary for $2$D KQ since the density is theoretically known to remain bounded. We obtained a polynomial decrease of the error.

3.3 Simulations of 3D KQ in radial coordinates

Here, we simulate equation (2.8) with $d=3$ for suitable choices of $\varepsilon,\delta$, $0<\varepsilon,\delta\ll 1$, where we choose $R_{1}=1$. We recall our notation $\bar{m}_{c}=\frac{1}{|\partial B(0,1)|}\int_{B(0,R_{1})}f_{c}(v)\,\mathrm{d}v$, where now $|\partial B(0,1)|=4\pi$ denotes the area of the 2-sphere, and remark that the numerical value of $\bar{m}_{c}$ is approximately given by $\bar{m}_{c}\approx 1.84$. We perform three simulations with a mass-supercritical, a mass-subcritical and a highly concentrated initial datum, respectively. More precisely, choosing as initial data again Gaussians of the form $f_{0}(v)=A\mathrm{e}^{-|v|^{2}/(2\sigma)}$, we run our scheme with the following three sets of parameters:

1. (P5)

$m<m_{c}:$ $\sigma=0.3$, $A=3$, $T=0.2$, $\tau=0.001$, $n=2001$, $\varepsilon=0$, $\delta=0$. 2. (P6)

$m>m_{c}:$ $\sigma=0.9$, $A=10$, $T=0.25$, $\tau=5\cdot 10^{-6}$, $n=50001$, $\varepsilon=10^{-12}$, $\delta=0$. 3. (P7)

$m<m_{c}:$ $\sigma=0.15$, $A=50$, $T=0.25$, $\tau=5\cdot 10^{-5}$, $n=2001$, $\varepsilon=10^{-10}$, $\delta=10^{-10}$.

The quantity $\bar{m}:=m/|\partial B(0,1)|$ associated with the above choice of parameters takes the value $\bar{m}\approx 0.335$ for (P5), $\bar{m}\approx 2.59$ for (P6), and $\bar{m}\approx 1.41$ for (P7) (see Figures 4, 5 and 6).

The size of the condensate divided by $|\partial B(0,1)|$, i.e. $\bar{x}_{p}(t):=\mathcal{L}^{1}(\{S(t,\cdot)=0\})$, is numerically determined by replacing the condition $S(t,\cdot)=0$ with the smallness criterion $S(t,\cdot)<10^{-10}$.

Remark 3.4*.*

The choice of the comparatively fine mesh in (P6) was made in order to ensure a sufficiently good approximation of the evolution of the entropy. See Fig. 5, which suggests an exponential decay.

Long-time behaviour.

Our simulations suggest that $3$D KQ has properties which are very similar to the Fokker–Planck model for bosons in $1$D in the $L^{1}$-supercritical regime. Figures 4, 5 and 6 suggest that in the long-time limit the numerical solution $S(t,\cdot)$ approximates the minimiser of the entropy (at the level of $S$), which we here denote111For simplicity, in our notation $S_{\infty}$ for the entropy minimiser we omit its dependence on the given mass $m$ and the radius $R_{1}$. by $S_{\infty}$.

Entropy.

The decay of the relative entropy appears to be exponential in all three cases (P5)–(P7), see Figures 4, 5 and 6. In each of these plots the red slope indicates the approximate slope of the graph averaged over the interval where it is plotted. Numerically, the relative entropy $H(t)-H_{\infty}$ appears to decay to zero like $e^{-\alpha t}$, where $\alpha\approx 35.3$ for (P5), $\alpha\approx 21.1$ for (P6), and $\alpha\approx 21.7$ for (P7).

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,\cdot)$ at height zero after some finite time, see Fig. 5 and 6. In the original variables this means that mass is gradually absorbed by the origin. Furthermore, Fig. 6 shows that, similarly to the observations in $1$D (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, at least numerically, the fraction of particles in the condensate is, in general, not monotonic in time for the 3D Kaniadakis–Quarati model.

Remark 3.5*.*

In order to produce the transient condensate in Figure 6 it was necessary to choose the parameter $\delta$ appearing in equation (2.8) (and its discrete counterpart) strictly positive. The same simulation for $\delta=0$ results in the flat part being trapped at height zero once it has formed. As explained in Section 2.1.2 and also in view of our results for the $1$D model, this “stickiness” appears to be a numerical artefact resulting from the circumstance that a regularisation based on a positive $\varepsilon$ but vanishing $\delta$ is imbalanced and favours condensation.

Blow-up profile.

At times where the solution has a non-trivial condensate component, we were interested in the spatial behaviour of $S(t,\cdot)$ close to $\{S(t,\cdot)=0\}$. Owing to the results on the $1$D model, one may expect the function $f(t,\cdot)$ to behave to leading order like the limiting steady state $f_{c}$, i.e. like $2|v|^{-2}$. Furthermore, the formal expansions in [47, Section III.C] suggest that for isotropic solutions of 3D KQ the error by which $f(t,\cdot)$ deviates from $f_{c}$ has the form

[TABLE]

for some constant $c(t)\in\mathbb{R}$. Our experiments corroborate formula (3.6). Indeed, Figures 5 and 7 displaying the quantity $f(t,v)/f_{c}(v)$ at times where $f(t,\cdot)$ is unbounded at the origin show that numerically it behaves like $1+\tilde{c}(t)|v|+o(|v|)$ as $|v|\to 0$.

4 Conclusion

In this work we propose a numerical scheme for nonlinear Fokker–Planck equations for bosons 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 considering appropriately normalised pseudo-inverse distributions 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. The physical entropy associated to these problems is shown to be decreasing for the semidiscrete schemes in 1D and 2D. We illustrate different phenomena appearing in the 3D radial KQ model mimicking the phenomena observed and partially proved for the 1D caricature of the KQ model in the $L^{1}$-supercritical case, see [15].

Acknowledgements

JAC was partially supported by the EPSRC grant number EP/P031587/1. KH was supported by MASDOC DTC at the University of Warwick, which is funded by the EPSRC grant number EP/HO23364/1. MTW acknowledges partial support by the EPSRC grant number EP/P01240X/1.