# A quantitative DiPerna-Lions theory and exponential mixing rate for passive
scalar transport111Date: September 1, 2024.222Keywords: Passive scalar
transport equation; DiPerna-Lions well-posedness theory, exponential bound on
mixing rate for transport along Sobolev vector field, uniform decay rate of
DiPerna-Lions commutator.333AMS Subject Classification: 35F05, 35Q35, 76B03,
76F25, 76R10.

Lucas Huysmans Department of Pure Mathematics & Mathematical Statistics,
University of Cambridge, Cambridge CB3 0WB. Email address<EMAIL_ADDRESS>Ayman Rimah Said Department of Pure Mathematics & Mathematical Statistics,
University of Cambridge, Cambridge CB3 0WB. Email address<EMAIL_ADDRESS>
( )

###### Abstract

We consider the transport of a passive scalar $f\in\mathbb{R}$ along a
divergence-free velocity vector field $u\in\mathbb{R}^{d}$ on the infinite
space $\mathbb{R}^{d}$. We give a quantitative version of the DiPerna-Lions
well-posedness theory for Sobolev vector fields $u\in L_{t}^{1}W_{x}^{1,p}$
when $1<p<\infty$ by giving a uniform decay rate of the DiPerna-Lions
commutator. We recover a slightly more general form of the known exponential
bound on the mixing rate by a Sobolev vector field $u\in L_{t}^{1}W_{x}^{1,p}$
when $1<p<\infty$.

###### Contents

  1. 1 Introduction
  2. 2 Main results
  3. 3 Proofs
  4. Acknowledgements

##  1 Introduction

In this paper we study the passive scalar transport equation

$\frac{\partial f}{\partial t}+\nabla\cdot(uf)=0,$ (TE)

for a divergence-free velocity vector field
$u:\mathbb{R}^{d}\times[0,T]\to\mathbb{R}^{d}$ with Sobolev regularity
$L_{t}^{1}W_{x}^{1,p}$, and a scalar field
$f:\mathbb{R}^{d}\times[0,T]\to\mathbb{R}$. The well-posedness of the initial
value problem to (TE) is well-studied (see the lecture notes of
[Ambrosio2014]), and is classically approached by the method of characteristic
curves when $u$ is Lipschitz-continuous. The seminal paper
[diperna1989ordinary] of DiPerna and Lions extends well-posedness to weak
solutions under the more general criterion $u\in L_{t}^{1}W_{x}^{1,p}$, see
Theorem 2.4. However, in contrast to the method of characteristic curves, the
well-posedness proof of DiPerna and Lions is non-quantitative. This prevents,
for example, from giving uniform error estimates for approximate solutions to
(TE), such as those given by numerical approximation.

The main contributions of this paper are two-fold. Firstly, we give a
quantitative version of DiPerna and Lions’ well-posedness proof when $u\in
L_{t}^{1}W_{x}^{1,p}$ for $1<p<\infty$, see section 1.3 and Theorem 2.5, which
we show is optimal for $p=2$ and derive a new stability estimate (2).
Secondly, we extend this to give in (4) a slightly more general version of the
$L_{t}^{1}W_{x}^{1,p}$-exponential bound on the geometric mixing scale of
transport, originally due to Crippa and DeLellis in [crippa2008estimates], see
Section 1.4 and Theorem 2.8.

We mention that one may prove analogous results for transport on the torus
$\mathbb{R}^{d}/\mathbb{Z}^{d}$, or along vector fields $u\in
L_{t}^{1}W_{x}^{1,p}$ with bounded divergence $\nabla\cdot u\in
L_{t}^{1}L_{x}^{\infty}$ (as considered in the original paper
[diperna1989ordinary] of DiPerna and Lions).

###  1.1 The DiPerna-Lions commutator

The heart of the well-posedness proof of DiPerna and Lions is to introduce the
so-called DiPerna-Lions commutator,

$r_{\delta}(u,f)=\int_{\mathbb{R}^{d}}f(y)(u(x)-u(y))\cdot\nabla\rho_{\delta}(x-y)\;dy,$

where we suppress any time-dependence, and
$\rho_{\delta}(x)=\delta^{-d}\rho\left(\frac{x}{\delta}\right)\in
C^{\infty}(\mathbb{R}^{d})$ is a choice of spatial mollifier, see Definition
2.2. Then $f*\rho_{\delta}$ solves the classical PDE

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla\right)(f*\rho_{\delta})=r_{\delta}(u,f).$

The key theorem of DiPerna and Lions, from which well-posedness of (TE)
follows, is

######  Theorem 1.2 (DiPerna-Lions [diperna1989ordinary]).

Let $u\in W^{1,p}(\mathbb{R}^{d};\mathbb{R}^{d})$ for some $p\in[1,\infty]$ be
divergence-free, and $f\in L^{p^{*}}(\mathbb{R}^{d})$, then

$\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}\xrightarrow{\delta\to
0}0.$

Our goal is to show that there exists a uniform rate of convergence,
independent of $u$ and $f$, for
$\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$.

###  1.3 Quantitative decay rate

We show the following integral decay rate for
$\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$,

$\left(\int_{0}^{\infty}\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}^{2}\frac{d\delta}{\delta}\right)^{\frac{1}{2}}\leq
C\left\|\nabla
u\right\|_{L^{2}(\mathbb{R}^{d})}\left\|f\right\|_{L^{2}(\mathbb{R}^{d})},$
(1)

see Theorem 2.5 and Corollary 2.6 for the precise statement, including the
modified rate when $\nabla u\in L^{p}(\mathbb{R}^{d})$ for $1<p<\infty$. Note
that the singular nature of the measure $\frac{d\delta}{\delta}$ implies decay
of $\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$ as $\delta\to 0$.

Moreover, we show in Proposition 2.7 that (1) is optimal, and in particular
that a standard uniform convergence rate of the form
$\sup_{u,f}\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}\to 0$ does
exist even for smooth $u$.

We observe that the decay rate (1) is sufficient to give uniform stability
estimates for approximate solutions to (TE). To illustrate this, suppose
$\bar{f}$ solves

$\frac{\partial\bar{f}}{\partial t}+\nabla\cdot(u\bar{f})=\phi,$

for some $\phi$ small in weak-$L^{1}$, and the same initial date $f(x,0)$.
Then

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla\right)(\bar{f}_{\delta}-f_{\delta})=\phi_{\delta}+r_{\delta}(u,\bar{f}-f),$

for all $\delta>0$, where we use the shorthand $f_{\delta}=f*\rho_{\delta}$.
Then by standard energy estimates

$\left\|\bar{f}_{\delta}-f_{\delta}\right\|_{L_{t}^{\infty}L_{x}^{1}}\leq\left\|\phi_{\delta}\right\|_{L_{t}^{1}L_{x}^{1}}+\left\|r_{\delta}(u,\bar{f}-f)\right\|_{L_{t}^{1}L_{x}^{1}},$

for all $\delta^{\prime}>0$. To make use of the decay rate (1), we must
integrate over $\delta$. For any $0<\epsilon<\delta$, we see that

$\int_{\epsilon}^{\delta}\left\|\bar{f}_{\delta^{\prime}}-f_{\delta^{\prime}}\right\|_{L_{t}^{\infty}L_{x}^{1}}\frac{d\delta^{\prime}}{\delta^{\prime}}\leq\int_{\epsilon}^{\delta}\left\|\phi_{\delta^{\prime}}\right\|_{L_{t}^{1}L_{x}^{1}}\frac{d\delta^{\prime}}{\delta^{\prime}}+\sqrt{\log\left(\frac{\delta}{\epsilon}\right)}\left(\int_{\epsilon}^{\delta}\left\|r_{\delta^{\prime}}(u,\bar{f}-f)\right\|_{L_{t}^{1}L_{x}^{1}}^{2}\frac{d\delta^{\prime}}{\delta^{\prime}}\right)^{\frac{1}{2}}.$

Then by (1) and Minkowski’s integral inequality, for a Gaussian mollifier (so
that for convenience,
$\left\|g*\rho_{\delta_{1}}\right\|_{L^{1}(\mathbb{R}^{d})}\leq\left\|g*\rho_{\delta_{2}}\right\|_{L^{1}(\mathbb{R}^{d})}$
for all $g\in L^{1}(\mathbb{R}^{d})$ whenever $\delta_{1}\geq\delta_{2}$), we
obtain the desired stability estimate

$\left\|\bar{f}_{\delta}-f_{\delta}\right\|_{L_{t}^{\infty}L_{x}^{1}}\leq\left\|\phi_{\epsilon}\right\|_{L_{t}^{1}L_{x}^{1}}+\frac{C}{\sqrt{\log\left(\frac{\delta}{\epsilon}\right)}}\left\|\nabla
u\right\|_{L_{t}^{1}L_{x}^{2}}\left\|\bar{f}-f\right\|_{L_{t}^{\infty}L_{x}^{2}},$
(2)

for all $\delta>\epsilon$.

###  1.4 Connection with mixing

‘Mixing’ of a passive scalar refers to the transfer of mass to smaller spatial
scales/higher frequencies. One approach to measure mixing is through decay of
weak norms such as $\left\|f\right\|_{H^{-1}(\mathbb{R}^{d})}$, and is usually
referred to as the functional mixing scale. A stronger notion is the geometric
mixing scale $\delta\in(0,\infty)$, which is defined by conditions of the form
$\left\|f*\rho_{\delta}\right\|_{L^{\infty}(\mathbb{R}^{d})}\leq\kappa$, see
for example [bressan2003lemma]. In particular, the DiPerna-Lions commutator
controls exactly the change in such geometric mixing scales, since

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla\right)(f*\rho_{\delta})=r_{\delta}(u,f).$

Therefore, bounds on the commutator immediately give bounds on the permitted
growth/decay of the mixing scale. For example, one may show through standard
energy estimates as before

$\left\|f_{T}*\rho_{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}\geq\left\|f_{0}*\rho_{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}-\left\|r_{\delta^{\prime}}(u,f)\right\|_{L_{t}^{1}L_{x}^{1}},$

for all $\delta^{\prime}>0$, where we use the shorthand $f_{t}(x)=f(x,t)$.
Then, as a corollary of (1) we have as for the previous computations

$\left\|f_{T}*\rho_{\delta}\right\|_{L^{1}(\mathbb{R}^{d})}\geq\left\|f_{0}*\rho_{\epsilon}\right\|_{L^{1}(\mathbb{R}^{d})}-\frac{C}{\sqrt{\log\left(\frac{\epsilon}{\delta}\right)}}\left\|\nabla
u\right\|_{L_{t}^{1}L_{x}^{2}}\left\|f_{0}\right\|_{L^{2}(\mathbb{R}^{d})},$

for all $\delta<\epsilon$, (where we have again supposed for simplicity $\rho$
is a Gaussian).

If we define the final geometric mixing scale $\delta\in(0,\infty)$ by the
mixing condition

$\left\|f_{T}*\rho_{\delta}\right\|_{L^{1}(\mathbb{R}^{d})}\leq\left\|f_{0}*\rho_{\epsilon}\right\|_{L^{1}(\mathbb{R}^{d})}-\kappa,$

for some $\kappa>0$, this then gives a Gaussian lower bound on the mixing
scale,

$\delta\geq\epsilon\exp\left(-\frac{C^{2}}{\kappa^{2}}\left\|f_{0}\right\|_{L^{2}(\mathbb{R}^{d})}^{2}\left\|\nabla
u\right\|_{L_{t}^{1}L_{x}^{2}}^{2}\right),$

and for this reason we refer to (1) as a ‘Gaussian’ decay rate on the
commutator.

Interestingly, the mixing scale is known to be bounded below by an exponential
in $\left\|\nabla u\right\|_{L_{t}^{1}L_{x}^{p}}$ when $1<p<\infty$, see
[crippa2008estimates, leger2018new, hadvzic2018singular]. Meanwhile, we have
that ‘Gaussian’ decay of the commutator (1) is optimal (see Proposition 2.7).
This discrepancy can be rectified by considering a different kind of decay.

Our second main result is the bound

$\left\|\int_{\delta}^{\epsilon}r_{\delta^{\prime}}(u,f)(x)\;\frac{d\delta^{\prime}}{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}\leq
C_{p}\left\|\nabla
u\right\|_{L^{p}(\mathbb{R}^{d})}\left\|f\right\|_{L^{p^{*}}(\mathbb{R}^{d})},$
(3)

for all $1<p<\infty$, where $\frac{1}{p}+\frac{1}{p^{*}}=1$, and any
$\delta<\epsilon$, see Theorem 2.8 for the precise statement. That is, by
reversing the order of integration over $L^{1}(\mathbb{R}^{d})$ and
$\frac{d\delta}{\delta}$, we recover the ‘exponential’ decay rate of the
commutator as suggested by the exponential mixing bounds [crippa2008estimates,
leger2018new, hadvzic2018singular].

Indeed, as a corollary of (3), we may reproduce these exponential mixing
bounds. For any $\delta<\epsilon$ we have transport of

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla\right)\left(\int_{\delta}^{\epsilon}f*\rho_{\delta^{\prime}}\;\frac{d\delta^{\prime}}{\delta^{\prime}}\right)=\int_{\delta}^{\epsilon}r_{\delta^{\prime}}(u,f)\;\frac{d\delta^{\prime}}{\delta^{\prime}},$

and so by standard energy estimates,

$\left\|\int_{\delta}^{\epsilon}f_{T}*\rho_{\delta^{\prime}}\;\frac{d\delta^{\prime}}{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}\geq\left\|\int_{\delta}^{\epsilon}f_{0}*\rho_{\delta^{\prime}}\;\frac{d\delta^{\prime}}{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}-\left\|\int_{\delta}^{\epsilon}r_{\delta^{\prime}}(u,f)\;\frac{d\delta^{\prime}}{\delta^{\prime}}\right\|_{L_{t}^{1}L_{x}^{1}},$

where we again use the shorthand $f_{t}(x)=f(x,t)$. Then as a corollary of
(3), and if we take for simplicity $\rho$ to be a Gaussian, we obtain the
following lower bound

$\left\|f_{T}*\rho_{\delta}\right\|_{L^{1}(\mathbb{R}^{d})}\\\ \begin{aligned}
&\qquad\geq\frac{1}{\log\left(\frac{\epsilon}{\delta}\right)}\left(\left\|\int_{\delta}^{\epsilon}f_{0}\;\frac{d\delta^{\prime}}{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}-\left\|\int_{\delta}^{\epsilon}\left(f_{0}-f_{0}*\rho_{\delta^{\prime}}\right)\;\frac{d\delta^{\prime}}{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}-C_{p}\left\|\nabla
u\right\|_{L_{t}^{1}L_{x}^{p}}\left\|f_{0}\right\|_{L^{p^{*}}(\mathbb{R}^{d})}\right)\\\
&\qquad\geq\left(\left\|f_{0}\right\|_{L^{1}(\mathbb{R}^{d})}-\sup_{\delta^{\prime}<\epsilon}\left\|f_{0}-f_{0}*\rho_{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}\right)-\frac{C_{p}}{\log\left(\frac{\epsilon}{\delta}\right)}\left\|\nabla
u\right\|_{L_{t}^{1}L_{x}^{p}}\left\|f_{0}\right\|_{L^{p^{*}}(\mathbb{R}^{d})},\end{aligned}$

for all $\delta<\epsilon$.

Then, under the mixing condition

$\left\|f_{T}*\rho_{\delta}\right\|_{L^{1}(\mathbb{R}^{d})}\leq\left(\left\|f_{0}\right\|_{L^{1}(\mathbb{R}^{d})}-\sup_{\delta^{\prime}<\epsilon}\left\|f_{0}-f_{0}*\rho_{\delta^{\prime}}\right\|_{L^{1}(\mathbb{R}^{d})}\right)-\kappa,$

for some $\kappa>0$, we deduce the desired exponential bound on the geometric
mixing scale

$\delta\geq\epsilon\exp\left(-\frac{C_{p}}{\kappa}\left\|f_{0}\right\|_{L^{p^{*}}(\mathbb{R}^{d})}\left\|\nabla
u\right\|_{L_{t}^{1}L_{x}^{p}}\right).$ (4)

This generalises the exponential mixing bound of Crippa and DeLellis in
[crippa2008estimates] to all passive scalars $f_{0}\in
L^{1}(\mathbb{R}^{d})\cap L^{p^{*}}(\mathbb{R}^{d})$.

We mention that, from the continuous embedding
$L^{\infty}(\mathbb{R}^{d}/\mathbb{Z}^{d})\hookrightarrow
L^{1}(\mathbb{R}^{d}/\mathbb{Z}^{d})$, the analogous result on the torus gives
an exponential bound on the geometric mixing scale defined instead by
$\left\|f_{T}*\rho_{\delta}\right\|_{L^{\infty}(\mathbb{R}^{d})}\leq\kappa$,
as in Bressan’s conjecture [bressan2003lemma].

######  Remark 1.5.

When sharing a version of this manuscript [eliaprivate] we were made aware of
an unpublished result similar to (3) by Brué, Colombo and De Philippis.

###  1.6 Further comments

One may ask if these results extend to $u\in L_{t}^{1}W_{x}^{1,1}$. This is
particularly important in regards to the mixing conjecture of Bressan in
[bressan2003lemma]. Harmonic estimates of the kind (3) are shown to fail for
$\left\|\nabla u\right\|_{L^{1}(\mathbb{R}^{d})}$ in [hadvzic2018singular],
but are shown in [seeger2019multilinear] to be valid if $\left\|\nabla
u\right\|_{L^{1}(\mathbb{R}^{d})}$ is replaced by the Hardy space
$\left\|\nabla u\right\|_{h^{1}(\mathbb{R}^{d})}$.

On the other hand, the quantitative DiPerna-Lions estimate (1) relies instead
on stability of $\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$ under
convergence of $\nabla u,f$ in weak-$L^{p}(\mathbb{R}^{d})$, and weak
compactness of the $L^{p}$-unit ball. To illustrate this, one has the
stability estimate

$\left\|r_{\delta}(u,f)-r_{\delta}(\bar{u},f)\right\|_{L^{1}(\mathbb{R}^{d})}\leq\frac{1}{\delta}\left\|u-\bar{u}\right\|_{L^{p}(\mathbb{R}^{d})}\left\|f\right\|_{L^{p^{*}}(\mathbb{R}^{d})},$

for any $\bar{u}\in W^{1,p}(\mathbb{R}^{d})$, where also
$\frac{1}{p}+\frac{1}{p^{*}}=1$, and we have assumed for simplicity
$\left\|\nabla\rho\right\|_{L^{1}(\mathbb{R}^{d})}=1$. We know from Theorem
1.2 (by DiPerna and Lions) that for fixed $\bar{u}\in
W^{1,p}(\mathbb{R}^{d})$, $f\in L^{p^{*}}(\mathbb{R}^{d})$, then
$\left\|r_{\delta}(\bar{u},f)\right\|_{L^{1}(\mathbb{R}^{d})}\xrightarrow{\delta\to
0}0$. By the above stability estimate, for sufficiently small $\delta>0$ there
exists some weak ball
$\left\|u-\bar{u}\right\|_{L^{p}(\mathbb{R}^{d})}\leq\epsilon$ around
$\bar{u}\in W^{1,p}(\mathbb{R}^{d})$ such that
$\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$ is also small,
depending only on $\bar{u}$ and $\epsilon$. In this way we may produce a weak
open cover of $W^{1,p}(\mathbb{R}^{d})$, and for each open set a $\delta>0$
such that $\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$ is small
for all $u$ in that open set.

If some $S\subset W^{1,p}$ is weakly compact, then we may find a finite
subcover covering $S$, and so obtain a uniform $\delta_{\mathrm{min}}>0$ for
which $\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$ is small for
some $\delta>\delta_{\mathrm{min}}$ depending on $u\in S$. That is we obtain a
uniform decay rate for
$\inf_{\delta^{\prime}>\delta}\left\|r_{\delta^{\prime}}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}$,
of which the integral decay rate (1) is an example. It will turn out that the
optimal choice of $\bar{u}\in W^{1,p}(\mathbb{R}^{d})$, i.e. choice of weak
ball, corresponds exactly to real interpolation via the Lions-Peetre
$K$-function method, see [lions1964classe] and the proof of Theorem 2.5.

For the case $p=1$, one has that the unit ball $S\subset
W^{1,1}(\mathbb{R}^{d})$ fails to be weakly compact, and so the required
finite subcover cannot be found. On the other hand, one instead has weak
compactness of the $BV(\mathbb{R}^{d})$-unit ball. It should therefore be
possible to employ the above strategy to give a quantitative proof of
Ambrosio’s $BV$ well-posedness theory in [ambrosio2004transport]. The
resulting bound on the geometric mixing scale found by the authors via this
method was of the form
$\delta\gtrsim\left(\exp^{C\left\|u\right\|_{L_{t}^{1}BV_{x}}}(1)\right)^{-1}$,
where $\exp^{n}(x)$ refers to repeated exponentiation (tetration),
$\exp(\exp(...(\exp(x))))$.

##  2 Main results

We start with the following standard definitions.

######  Definition 2.1 (Weak solution to the transport equation (TE)).

Consider a vector field $u\in
L^{1}\left([0,T];L^{1}(\mathbb{R}^{d};\mathbb{R}^{d})\right)$ with
$\nabla\cdot u=0$ in the distributional sense.

We say $f\in L^{1}\left([0,T];L^{1}(\mathbb{R}^{d};\mathbb{R})\right)$ with
$uf\in L^{1}L^{1}$ is a weak solution to the transport equation along $u$

$\frac{\partial f}{\partial t}+\nabla\cdot(uf)=0,$ (TE)

with initial data $f_{0}\in L^{1}(\mathbb{R}^{d})$ if, for any $\phi\in
C_{c}^{\infty}(\mathbb{R}^{d}\times[0,T))$,

$\int_{\mathbb{R}^{d}\times[0,T)}f\left(\frac{\partial\phi}{\partial
t}+u\cdot\nabla\phi\right)\;dxdt=-\int_{\mathbb{R}^{d}}f_{0}\phi_{0}\;dx,$

where $\phi_{0}(x)=\phi(x,0)$.

######  Definition 2.2 (DiPerna-Lions commutator).

Let $u\in L^{1}(\mathbb{R}^{d};\mathbb{R}^{d})$, $f\in
L^{1}(\mathbb{R}^{d};\mathbb{R})$, with $uf\in
L^{1}(\mathbb{R}^{d};\mathbb{R}^{d})$. Let $\rho\in L^{1}(\mathbb{R}^{d})$
with

$\int_{\mathbb{R}^{d}}\rho\;dx=1\text{ and define }\rho_{\delta}\in
L^{1}(\mathbb{R}^{d}),\
\rho_{\delta}(x)=\delta^{-d}\rho\left(\frac{x}{\delta}\right).$

Then we define the DiPerna-Lions commutator

$r_{\delta}(u,f)(x)=\int_{\mathbb{R}^{d}}f(y)(u(x)-u(y))\cdot\nabla\rho_{\delta}(x-y)\;dy.$

######  Remark 2.3.

The motivation behind the DiPerna-Lions commutator is the following. Let $u\in
L^{1}([0,T];L^{1}(\mathbb{R}^{d};\mathbb{R}^{d}))$ be divergence-free, and
$f_{0}\in L^{1}(\mathbb{R}^{d})$. Let $f\in
L^{1}\left([0,T];L^{1}(\mathbb{R}^{d};\mathbb{R})\right)$ with $uf\in
L_{t}^{1}L_{x}^{1}$ be a solution to (TE) along $u$ with initial datum
$f_{0}$. Then $f*\rho_{\delta}$ solves the classical PDE

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla\right)(f*\rho_{\delta})=r_{\delta}(u,f).$

As an immediate corollary of Theorem 1.2 on the decay of $r_{\delta}(u,f)$ we
have

######  Corollary 2.4 (DiPerna-Lions [diperna1989ordinary]).

Let $u\in L^{1}([0,T];W^{1,p}(\mathbb{R}^{d};\mathbb{R}^{d}))$ be divergence-
free for some $p\in[1,\infty]$, and $f_{0}\in L^{p^{*}}(\mathbb{R}^{d})$ where
$\frac{1}{p}+\frac{1}{p^{*}}=1$. Then there exists a unique $f\in
L_{t}^{\infty}L_{x}^{p^{*}}$ such that $f$ is a solution to (TE) along $u$
with initial datum $f_{0}$.

One of the main results of our paper is a quantitative version of Theorem 1.2.
More precisely, we define the homogenous Besov space
$\dot{B}^{0}_{p,q}(\mathbb{R}^{d})$ as the Schwartz distributions
$\mathcal{S}^{\prime}(\mathbb{R}^{d})$ for which the following homogenous
Besov norm is finite,

$\left\|f\right\|_{\dot{B}^{0}_{p,q}(\mathbb{R}^{d})}=\left(\sum_{n=-\infty}^{\infty}\left\|\psi_{n}*f\right\|_{L^{p}(\mathbb{R}^{d})}^{q}\right)^{\frac{1}{q}},$

where $\psi_{n}\in C^{\infty}(\mathbb{R}^{d})$ is a Littlewood-Paley
decomposition, given by the Fourier transform

$\hat{\psi}_{n}(\xi)=\varphi(2^{-n}\xi),$

for a fixed choice of $\varphi\in C_{c}^{\infty}(\mathbb{R}^{d})$ satisfying
$\mathrm{supp}\;\varphi=\\{\xi:\frac{1}{2}\leq|\xi|\leq 2\\}$,
$\varphi(\xi)>0$ if $\frac{1}{2}<|\xi|<2$, and
$\sum_{n=-\infty}^{\infty}\varphi(2^{-n}\xi)=1$ for $\xi\neq 0$, as in
[bergh2012interpolation].

Recall the continuous embedding
$L^{p}(\mathbb{R}^{d})\hookrightarrow\dot{B}_{p,\mathrm{max}({p,2})}^{0}(\mathbb{R}^{d})$
for all $1<p<\infty$, see [bahouri2011fourier].

We prove the following uniform decay rate of the DiPerna-Lions commutator.

######  Theorem 2.5.

Let $p,p^{\prime},q,r\in[1,\infty]$ with
$\frac{1}{p}+\frac{1}{p^{\prime}}=\frac{1}{r}$.

Let $\nabla u\in
L^{p}(\mathbb{R}^{d};\mathbb{R}^{d})\cap\dot{B}^{0}_{p,q}(\mathbb{R}^{d};\mathbb{R}^{d})$
be divergence-free, and $f\in
L^{p^{\prime}}(\mathbb{R}^{d})\cap\dot{B}^{0}_{p^{\prime},q}(\mathbb{R}^{d})$,
then

$\left(\int_{0}^{\infty}\left\|r_{\delta}(u,f)\right\|_{L^{r}}^{q}\;\frac{d\delta}{\delta}\right)^{\frac{1}{q}}\leq
C\left(\left\|\nabla
u\right\|_{\dot{B}^{0}_{p,q}}\left\|f\right\|_{L^{p^{\prime}}}+\left\|\nabla
u\right\|_{L^{p}}\left\|f\right\|_{\dot{B}^{0}_{p^{\prime},q}}\right),$

for some constant $C>0$ depending only on the dimension $d$, and the following
norm of the mollifier
$\left\|(1+|x|^{2})\nabla\rho\right\|_{L^{1}}+\left\||x|\nabla^{2}\rho\right\|_{L^{1}}$.

By the continuous embedding
$L^{p}(\mathbb{R}^{d})\hookrightarrow\dot{B}_{p,\mathrm{max}({p,2})}^{0}(\mathbb{R}^{d})$
for all $1<p<\infty$, and Minkowski’s integral inequality we have the
following simple corollary.

######  Corollary 2.6.

Let $u\in L^{1}([0,T];W^{1,p}(\mathbb{R}^{d};\mathbb{R}^{d}))$ be divergence-
free. Let $f\in L^{\infty}([0,T];L^{p^{\prime}}(\mathbb{R}^{d}))$ with
$\frac{1}{r}=\frac{1}{p}+\frac{1}{p^{\prime}}\leq 1$. If
$p,p^{\prime}<\infty$, then

$\left(\int_{0}^{\infty}\left\|r_{\delta}(u,f)\right\|_{L_{t}^{1}L_{x}^{r}}^{\mathrm{max}(p,p^{\prime})}\;\frac{d\delta}{\delta}\right)^{\frac{1}{{\mathrm{max}(p,p^{\prime})}}}\leq
C\left\|\nabla
u\right\|_{L_{t}^{1}L_{x}^{p}}\left\|f\right\|_{L_{t}^{\infty}L_{x}^{p^{\prime}}},$

for some constant $C>0$ depending only on $p,p^{\prime}$, the dimension $d$,
and the following norm of the mollifier
$\left\|(1+|x|^{2})\nabla\rho\right\|_{L^{1}}+\left\||x|\nabla^{2}\rho\right\|_{L^{1}}$.

Next, we give an example showing that the results of the previous theorem are
optimal in the case $p=p^{\prime}=2$.

######  Proposition 2.7.

Let $u\in C^{\infty}(\mathbb{R}^{2};\mathbb{R}^{2})$ be the linear shear

$u(x,t)=(x_{2},0)\in\mathbb{R}^{2}.$

For any compactly supported mollifier $\rho\in
C_{c}^{\infty}(\mathbb{R}^{2})$, then both

$\liminf_{\delta\to
0}\left(\sup_{|f|\leq\phi}\left\|r_{\delta}(u,f)\right\|_{L^{1}(\mathbb{R}^{d})}\right)>0,\text{
and
}\inf_{q\in[1,2)}\left(\int_{0}^{1}\left\|r_{\delta}(u,\bar{f})\right\|_{L^{1}(\mathbb{R}^{d})}^{q}\frac{d\delta}{\delta}\right)^{\frac{1}{q}}=\infty,$

where $\bar{f}\in L^{p}(\mathbb{R}^{2})$ for all $p\in[1,\infty)$ and is
supported on $|x|\leq 1$, and for any cutoff $\phi\in
C_{c}^{\infty}(\mathbb{R}^{2})$ with $\phi(x)=1$ whenever $|x|\leq 1$.

Finally, the next theorem gives the control on the DiPerna-Lions commutator
that transfers to exponential mixing bounds.

######  Theorem 2.8.

Let $1<p,p^{\prime}\leq\infty$, and $1\leq r<\infty$ with
$\frac{1}{r}=\frac{1}{p}+\frac{1}{p^{\prime}}$.

Let $u\in W^{1,p}(\mathbb{R}^{d})$ be divergence-free, and $f\in
L^{p^{\prime}}(\mathbb{R}^{d})$, then for any $0<a<b$,

$\left\|\int_{a}^{b}r_{\delta}(u,f)(x)\;\frac{d\delta}{\delta}\right\|_{L^{r}(\mathbb{R}^{d})}\leq
C\left\|\nabla u\right\|_{L^{p}}\left\|f\right\|_{L^{p^{\prime}}},$

for some constant $C$ depending only on $p,p^{\prime}$, the dimension $d$, and
some high-order norm of the mollifier $\rho$ (and in particular not on $a,b$).

In particular, this includes any mollifier $\rho$ a Schwartz function.

##  3 Proofs

###### Proof of Theorem 2.5.

Throughout $\lesssim$ will denote less than or equal to up to a constant
depending only on the dimension $d$, and the norm of the mollifier
$\left\|(1+|x|^{2})\nabla\rho\right\|_{L^{1}}+\left\||x|\nabla^{2}\rho\right\|_{L^{1}}$,
and in particular not on $\delta\in(0,\infty)$.

We first rewrite

$\displaystyle r_{\delta}(u,f)(x)$
$\displaystyle=\int_{\mathbb{R}^{d}}f(y)(u(x)-u(y))\cdot\nabla\rho_{\delta}(x-y)\;dy,$
$\displaystyle=\int_{\mathbb{R}^{d}}f(x-\delta h)\left(\frac{u(x)-u(x-\delta
h)}{\delta}\right)\cdot\nabla\rho(h)\;dh,$

and so by testing against $\phi(x)\in L^{r^{*}}(\mathbb{R}^{d})$ we see
immediately that

$\left\|r_{\delta}(u,f)\right\|_{L^{r}}\lesssim\frac{1}{\delta}\left\|f\right\|_{L^{p^{\prime}}}\left\|u\right\|_{L^{p}}.$
(5)

Next, we consider the following approximation

$\displaystyle r_{\delta}(u,f)(x)-\sum^{d}_{i,j=1}\frac{\partial
u_{i}}{\partial x_{j}}(x)\int_{\mathbb{R}^{d}}f(x-\delta
h)h_{j}\frac{\partial\rho}{\partial h_{i}}\;dh$
$\displaystyle=\int_{\mathbb{R}^{d}}f(x-\delta
h)\sum^{d}_{i=1}\left(\frac{u_{i}(x)-u_{i}(x-\delta
h)}{\delta}-\sum^{d}_{j=1}h_{j}\frac{\partial u_{i}}{\partial
x_{j}}(x)\right)\frac{\partial\rho}{\partial h_{i}}\;dh$
$\displaystyle=\int_{0}^{1}\int_{\mathbb{R}^{d}}f(x-\delta
h)\sum^{d}_{i,j=1}\left(h_{j}\frac{\partial u_{i}}{\partial x_{j}}(x-t\delta
h)-h_{j}\frac{\partial u_{i}}{\partial
x_{j}}(x)\right)\frac{\partial\rho}{\partial h_{i}}\;dh\;dt$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\int_{\mathbb{R}^{d}}f(x-\delta
h)\sum^{d}_{i,j,k=1}\left(-t\delta h_{j}h_{k}\frac{\partial^{2}u_{i}}{\partial
x_{j}\partial x_{k}}(x-st\delta h)\right)\frac{\partial\rho}{\partial
h_{i}}\;dh\;dt\;ds,$

and so by testing against $\phi(x)\in L^{r^{*}}(\mathbb{R}^{d})$ we have

$\left\|r_{\delta}(u,f)-\sum^{d}_{i,j=1}\frac{\partial u_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}f(x-\delta h)h_{j}\frac{\partial\rho}{\partial
h_{i}}\;dh\right\|_{L^{r}}\lesssim\delta\left\|f\right\|_{L^{p^{\prime}}}\left\|\nabla^{2}u\right\|_{L^{p}}.$
(6)

Now notice, since $\nabla\cdot u=0$,

$\displaystyle\sum^{d}_{i,j=1}\frac{\partial u_{i}}{\partial
x_{j}}(x)\int_{\mathbb{R}^{d}}f(x-\delta h)h_{j}\frac{\partial\rho}{\partial
h_{i}}\;dh$ $\displaystyle=\sum^{d}_{i,j=1}\frac{\partial u_{i}}{\partial
x_{j}}(x)\int_{\mathbb{R}^{d}}f(x-\delta h)\frac{\partial(h_{j}\rho)}{\partial
h_{i}}\;dh$ $\displaystyle=\sum^{d}_{i,j=1}\frac{\partial u_{i}}{\partial
x_{j}}(x)\int_{\mathbb{R}^{d}}\delta\frac{\partial f}{\partial x_{i}}(x-\delta
h)h_{j}\rho\;dh,$

and so by testing against $\phi(x)\in L^{r^{*}}(\mathbb{R}^{d})$ we have the
bound

$\left\|\sum^{d}_{i,j=1}\frac{\partial u_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}f(x-\delta h)h_{j}\frac{\partial\rho}{\partial
h_{i}}\;dh\right\|_{L^{r}}\lesssim\delta\left\|\nabla
u\right\|_{L^{p}}\left\|\nabla f\right\|_{L^{p^{\prime}}}.$ (7)

Finally, assume $f=-\Delta g$ for some $g\in
W^{1,p^{\prime}}(\mathbb{R}^{d})$, then

$\displaystyle\sum^{d}_{i,j=1}\frac{\partial u_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}f(x-\delta h)\frac{\partial(h_{j}\rho)}{\partial
h_{i}}\;dh$ $\displaystyle=\sum^{d}_{i,j,k=1}\frac{\partial u_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}\frac{\partial^{2}g}{\partial x_{k}\partial
x_{k}}(x-\delta h)\frac{\partial(h_{j}\rho)}{\partial h_{i}}\;dh$
$\displaystyle=\sum^{d}_{i,j,k=1}\frac{\partial u_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}-\frac{1}{\delta}\frac{\partial}{\partial
h_{k}}\frac{\partial g}{\partial x_{k}}(x-\delta
h)\frac{\partial(h_{j}\rho)}{\partial h_{i}}\;dh$
$\displaystyle=\sum^{d}_{i,j,k=1}\frac{\partial u_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}\frac{1}{\delta}\frac{\partial g}{\partial
x_{k}}(x-\delta h)\frac{\partial^{2}(h_{j}\rho)}{\partial h_{k}\partial
h_{i}}\;dh,$

and so

$\left\|\frac{\partial u_{i}}{\partial x_{j}}\int_{\mathbb{R}^{d}}f(x-\delta
h)h_{j}\frac{\partial\rho}{\partial
h_{i}}\;dh\right\|_{L^{r}}\lesssim\frac{1}{\delta}\left\|\nabla
u\right\|_{L^{p}}\left\|\nabla(-\Delta)^{-1}f\right\|_{L^{p^{\prime}}}.$ (8)

Putting eqs. (5) to (8) together gives, for any $\bar{u}\in W^{2,p}$,
$\bar{f}\in W^{1,p^{\prime}}$ such that $(-\Delta)^{-1}(f-\bar{f})\in
W^{1,p^{\prime}}$,

$\displaystyle\left\|r_{\delta}(u,f)\right\|_{L^{r}}$
$\displaystyle\leq\begin{aligned}
&\left\|r_{\delta}(u-\bar{u},f)\right\|_{L^{r}}+\left\|r_{\delta}(\bar{u},f)-\sum^{d}_{i,j=1}\frac{\partial\bar{u}_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}f(x-\delta h)h_{j}\frac{\partial\rho}{\partial
h_{i}}\;dh\right\|_{L^{r}}\\\
&+\left\|\sum^{d}_{i,j=1}\frac{\partial\bar{u}_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}(f-\bar{f})(x-\delta
h)h_{j}\frac{\partial\rho}{\partial
h_{i}}\;dh\right\|_{L^{r}}+\left\|\sum^{d}_{i,j=1}\frac{\partial\bar{u}_{i}}{\partial
x_{j}}\int_{\mathbb{R}^{d}}\bar{f}(x-\delta
h)h_{j}\frac{\partial\rho}{\partial h_{i}}\;dh\right\|_{L^{r}}\end{aligned}$
$\displaystyle\lesssim\begin{aligned}
&\left\|f\right\|_{L^{p^{\prime}}}\left(\frac{1}{\delta}\left\|u-\bar{u}\right\|_{L^{p}}+\delta\left\|\nabla^{2}\bar{u}\right\|_{L^{p}}\right)\\\
&+\left\|\nabla\bar{u}\right\|_{L^{p}}\left(\frac{1}{\delta}\left\|\nabla(-\Delta)^{-1}(f-\bar{f})\right\|_{L^{p^{\prime}}}+\delta\left\|\nabla\bar{f}\right\|_{L^{p^{\prime}}}\right).\end{aligned}$
(9)

We now need some standard real interpolation inequalities. For this we recall
the definition of $\dot{B}^{s}_{p,q}$ as the Schwartz distributions
$\mathcal{S}^{\prime}(\mathbb{R}^{d})$ such that the following norm is finite,

$\left\|f\right\|_{\dot{B}^{s}_{p,q}}=\left(\sum_{n=-\infty}^{\infty}2^{ns}\left\|\psi_{n}*f\right\|_{L^{p}}^{q}\right)^{\frac{1}{q}}.$

By Theorem $6.4.5$ page 152 of [bergh2012interpolation] the homogeneous Besov
spaces $\dot{B}^{1}_{p,q},\dot{B}^{0}_{p^{\prime},q}$ are the real
interpolation spaces between the spaces

$(\dot{B}^{0}_{p,1},\dot{B}^{2}_{p,1})_{\frac{1}{2},q}\text{ and
}(\dot{B}^{-1}_{p^{\prime},1},\dot{B}^{1}_{p^{\prime},1})_{\frac{1}{2},q},$

respectively. Thus it is standard that using the Lions-Peetre $K$-function
method [lions1964classe], we have the following characterisation of the norms
$\left\|\cdot\right\|_{\dot{B}^{1}_{p,q}}$ and
$\left\|\cdot\right\|_{\dot{B}^{0}_{p^{\prime},q}}$. Defining

$\displaystyle\left|f\right|_{\dot{B}^{1}_{p,q}}=\left(\int_{0}^{+\infty}\frac{\left(\inf_{\bar{f}\in\dot{B}^{2}_{p,1},f-\bar{f}\in\dot{B}^{0}_{p,1}}\left(\left\|f-\bar{f}\right\|_{\dot{B}^{0}_{p,1}}+t\left\|\bar{f}\right\|_{\dot{B}^{2}_{p,1}}\right)\right)^{q}}{t^{\frac{q}{2}}}\frac{dt}{t}\right)^{\frac{1}{q}},$
$\displaystyle\left|f\right|_{\dot{B}^{0}_{p^{\prime},q}}=\left(\int_{0}^{+\infty}\frac{\left(\inf_{\bar{f}\in\dot{B}^{1}_{p^{\prime},1},f-\bar{f}\in\dot{B}^{-1}_{p^{\prime},1}}\left(\left\|f-\bar{f}\right\|_{\dot{B}^{-1}_{p^{\prime},1}}+t\left\|\bar{f}\right\|_{\dot{B}^{1}_{p^{\prime},1}}\right)\right)^{q}}{t^{\frac{q}{2}}}\frac{dt}{t}\right)^{\frac{1}{q}},$

there exists a universal constant $C$, independent of $p,q\in[1,\infty]$, such
that

$\frac{1}{C}\left\|f\right\|_{\dot{B}^{1}_{p,q}}\leq\left|f\right|_{\dot{B}^{1}_{p,q}}\leq
C\left\|f\right\|_{\dot{B}^{1}_{p,q}}\text{ and
}\frac{1}{C}\left\|f\right\|_{\dot{B}^{0}_{p^{\prime},q}}\leq\left|f\right|_{\dot{B}^{0}_{p^{\prime},q}}\leq
C\left\|f\right\|_{\dot{B}^{0}_{p^{\prime},q}}.$

Thus there exists some choice of such $\bar{u}_{\delta},\bar{f}_{\delta}$ for
each $\delta\in(0,\infty)$ so that

$\displaystyle\left(\int_{0}^{\infty}\left(\frac{1}{\delta}\left\|u-\bar{u}_{\delta}\right\|_{\dot{B}^{0}_{p,1}}+\delta\left\|\bar{u}_{\delta}\right\|_{\dot{B}^{2}_{p,1}}\right)^{q}\frac{d\delta}{\delta}\right)^{\frac{1}{q}}\lesssim\left\|u\right\|_{\dot{B}^{1}_{p,q}},$
$\displaystyle\left(\int_{0}^{\infty}\left(\frac{1}{\delta}\left\|f-\bar{f}_{\delta}\right\|_{\dot{B}^{-1}_{p^{\prime},1}}+\delta\left\|\bar{f}_{\delta}\right\|_{\dot{B}^{1}_{p^{\prime},1}}\right)^{q}\frac{d\delta}{\delta}\right)^{\frac{1}{q}}\lesssim\left\|f\right\|_{\dot{B}^{0}_{p^{\prime},q}},$

for all $q\in[1,\infty]$. Moreover, $\bar{u}_{\delta}\in\dot{B}^{2}_{p,1}$ may
always be chosen of the form

$\bar{u}_{\delta}=\sum_{n=-\infty}^{K_{\delta}}\psi_{n}*u,$

for some $K_{\delta}\in\mathbb{Z}$, see [bergh2012interpolation], and so in
addition we have the bound
$\left\|\nabla\bar{u}_{\delta}\right\|_{L^{p}}\lesssim\left\|\nabla
u\right\|_{L^{p}}$ independent of $\delta$. In light of this, we may apply (9)
for each $\delta\in(0,\infty)$ to give the result as follows

$\displaystyle\left(\int_{0}^{\infty}\left\|r_{\delta}(u,f)\right\|_{L^{r}}^{q}\;\frac{d\delta}{\delta}\right)^{\frac{1}{q}}$
$\displaystyle\lesssim\begin{aligned}
&\left(\int_{0}^{\infty}\left(\left\|f\right\|_{L^{p^{\prime}}}\left(\frac{1}{\delta}\left\|u-\bar{u}_{\delta}\right\|_{L^{p}}+\delta\left\|\nabla^{2}\bar{u}_{\delta}\right\|_{L^{p}}\right)\right)^{q}\;\frac{d\delta}{\delta}\right)^{\frac{1}{q}}\\\
&+\left(\int_{0}^{\infty}\left(\left\|\nabla\bar{u}_{\delta}\right\|_{L^{p}}\left(\frac{1}{\delta}\left\|\nabla(-\Delta)^{-1}(f-\bar{f}_{\delta})\right\|_{L^{p^{\prime}}}+\delta\left\|\nabla\bar{f}_{\delta}\right\|_{L^{p^{\prime}}}\right)\right)^{q}\;\frac{d\delta}{\delta}\right)^{\frac{1}{q}}\end{aligned}$
$\displaystyle\lesssim\begin{aligned}
&\left\|f\right\|_{L^{p^{\prime}}}\left(\int_{0}^{\infty}\left(\frac{1}{\delta}\left\|u-\bar{u}_{\delta}\right\|_{\dot{B}^{0}_{p,1}}+\delta\left\|\bar{u}_{\delta}\right\|_{\dot{B}^{2}_{p,1}}\right)^{q}\;\frac{d\delta}{\delta}\right)^{\frac{1}{q}}\\\
&+\left\|\nabla
u\right\|_{L^{p}}\left(\int_{0}^{\infty}\left(\frac{1}{\delta}\left\|f-\bar{f}_{\delta}\right\|_{\dot{B}^{-1}_{p^{\prime},1}}+\delta\left\|\bar{f}_{\delta}\right\|_{\dot{B}^{1}_{p^{\prime},1}}\right)^{q}\;\frac{d\delta}{\delta}\right)^{\frac{1}{q}}\end{aligned}$
$\displaystyle\lesssim\left\|f\right\|_{L^{p^{\prime}}}\left\|u\right\|_{\dot{B}^{1}_{p,q}}+\left\|\nabla
u\right\|_{L^{p}}\left\|f\right\|_{\dot{B}^{0}_{p^{\prime},q}},$

where we have used by Theorem $6.2.4$ page $142$ of [bergh2012interpolation]
that we have the continuous embedding say $\dot{B}^{0}_{p,1}\hookrightarrow
L^{p}$ and we conclude the proof by
$\left\|u\right\|_{\dot{B}^{1}_{p,q}}\lesssim\left\|\nabla
u\right\|_{\dot{B}^{0}_{p,q}}$. ∎

###### Proof of Proposition 2.7.

For any $f\in L^{1}(\mathbb{R}^{2})$,

$\displaystyle r_{\delta}(u,f)(x)$
$\displaystyle=\int_{\mathbb{R}^{2}}f(y)(x_{2}-y_{2})\frac{\partial\rho_{\delta}}{\partial
x_{1}}(x-y)\;dy$
$\displaystyle=\int_{\mathbb{R}^{2}}f(x-h)\frac{\partial(h_{2}\rho_{\delta})}{\partial
h_{1}}(h)\;dh,$ $\displaystyle=f*K_{\delta},$

where $K_{\delta}(x)=\delta^{-2}K(\delta^{-1}x)$ and
$K(x)=\frac{\partial(x_{2}\rho)}{\partial x_{1}}(x)$.

In particular, the Fourier transforms of $K_{\delta}$ and $K$ are given by

$\displaystyle\hat{K}_{\delta}(\xi)=\hat{K}(\delta\xi),$
$\displaystyle\hat{K}(\xi)=-\xi_{1}\frac{\partial\hat{\rho}}{\partial\xi_{2}}.$

For now consider any $f\in\mathcal{S}^{\prime}(\mathbb{R}^{2})$ a Schwartz
distribution. Then $f*K_{\delta}$ will coincide with $r_{\delta}(u,f)$ if in
addition $f\in L^{1}(\mathbb{R}^{2})$.

To prove the first statement of the theorem we take the harmonic
$f(x)=e^{i\delta^{-1}\bar{\xi}\cdot x}$ for some $\bar{\xi}\in\mathbb{R}^{2}$
such that $|\hat{K}(\bar{\xi})|>0$, which exists since $\rho\in
L^{1}(\mathbb{R}^{2})$ with $\int_{\mathbb{R}^{2}}\rho\;dx=1$. Then
$f*K_{\delta}=e^{i\delta^{-1}\bar{\xi}\cdot x}\hat{K}(\bar{\xi})$.

We now cutoff $f$ so that it lies in $L^{1}(\mathbb{R}^{d})$. Let the
mollifier $\rho$ be supported on $|x|\leq R$, then $\rho_{\delta}$ is
supported on $|x|\leq\delta R$. Taking a compactly supported cutoff $\phi(x)$
with $\phi(x)=1$ on $|x|\leq 1$, then $r_{\delta}(u,\phi f)=f*K_{\delta}$
whenever $|x|\leq 1-\delta R$. Therefore, for all $\delta\leq\frac{1}{2R}$,
$\left\|r_{\delta}(u,\phi f)\right\|_{L^{1}}\geq|\hat{K}(\bar{\xi})|$ as
required.

Note that currently $\phi f$ is complex valued, but by taking either the real,
or imaginary part, a similar non-zero bound must hold.

To prove the second statement of the theorem, pick
$\bar{\xi}\in\mathbb{R}^{2}$ with $|\bar{\xi}|=1$, and in addition some
$0<c_{1}<c_{2}<2c_{1}$, such that
$\inf_{\delta\in[c_{1},c_{2}]}|\hat{K}(\delta\bar{\xi})|=\epsilon>0$, noting
that $\hat{K}(\xi)$ is continuous if we assume in addition $x_{2}\rho(x)\in
L^{1}(\mathbb{R}^{2})$.

Let $f\in\mathcal{S}^{\prime}(\mathbb{R}^{2})$ be of the form
$f(x)=\sum_{n=1}^{\infty}a_{n}e^{i\delta_{n}^{-1}\bar{\xi}\cdot x}$ for some
$a_{n}\in[0,1]$ and $\delta_{n}\in(0,\infty)$.

Now,

$\displaystyle\left|e^{i\delta_{n}^{-1}\bar{\xi}\cdot x}*K_{\delta}\right|$
$\displaystyle=\left|\hat{K}\left(\frac{\delta}{\delta_{n}}\bar{\xi}\right)\right|$
$\displaystyle\leq\min\left\\{\frac{\delta}{\delta_{n}}\left\|\frac{\partial\hat{\rho}}{\partial\xi_{2}}\right\|_{L^{\infty}},\frac{\delta_{n}}{\delta}\left\||\xi|^{2}\frac{\partial\hat{\rho}}{\partial\xi_{2}}\right\|_{L^{\infty}}\right\\}$
$\displaystyle\leq
C\min\left\\{\frac{\delta}{\delta_{n}},\frac{\delta_{n}}{\delta}\right\\}.$

Now let $\delta_{n}=2^{-n^{2}}$. Then for
$\frac{\delta}{\delta_{n}}\in[c_{1},c_{2}]\subset[c_{1},2c_{1}]$,

$\displaystyle|f*K_{\delta}|$ $\displaystyle\geq
a_{n}\left|\hat{K}\left(\frac{\delta}{\delta_{n}}\bar{\xi}\right)\right|-C\sum_{1\leq
m\leq
n-1}|a_{m}|\frac{c_{1}2^{-n^{2}}}{2^{-m^{2}}}-C\sum_{m=n+1}^{\infty}|a_{m}|\frac{2^{-m^{2}}}{c_{2}2^{-n^{2}}}$
$\displaystyle\geq
a_{n}\left|\hat{K}\left(\frac{\delta}{\delta_{n}}\bar{\xi}\right)\right|-c_{1}C2^{2-2n}-c_{2}^{-1}C2^{-2n}$
$\displaystyle\geq a_{n}\epsilon-C(c_{1}+c_{2}^{-1})4^{1-n}.$

Therefore, for any $q\in[1,2)$,

$\displaystyle\int_{\delta_{n}c_{1}}^{2\delta_{n}c_{1}}\left(\int_{\left[-\frac{1}{2},\frac{1}{2}\right]^{2}}|f*K_{\delta}|(x)\;dx\right)^{q}\frac{d\delta}{\delta}$
$\displaystyle\geq\left(\frac{1}{2}a_{n}^{q}\epsilon^{q}-(C(c_{1}+c_{2}^{-1})4^{1-n})^{q}\right)\log\left(\frac{c_{2}}{c_{1}}\right),$

and so, for $N\in\mathbb{N}$ such that $2\delta_{N}c_{1}\leq 1$,

$\displaystyle\int_{0}^{1}\left(\int_{\left[-\frac{1}{2},\frac{1}{2}\right]^{2}}|f*K_{\delta}|(x)\;dx\right)^{q}\frac{d\delta}{\delta}\geq\frac{1}{2}\log\left(\frac{c_{2}}{c_{1}}\right)\epsilon^{q}\sum_{n=N}^{\infty}a_{n}^{q}-\frac{C^{q}(c_{1}+c_{2}^{-1})^{q}}{1-4^{-q}}\log\left(\frac{c_{2}}{c_{1}}\right).$

Now, choose $\\{a_{n}\\}_{n=1}^{\infty}$ so that
$\sum_{n=N}^{\infty}a_{n}^{q}=\infty$ but
$\sum_{n=N}^{\infty}a_{n}^{2}<\infty$. Then
$f\in\dot{B}^{0}_{p,2}([-\frac{1}{2},\frac{1}{2}]^{2})$ for all
$p\in[1,\infty)$ (noting that $f$ is periodic with period 1, and each
$e^{i\delta_{n}^{-1}\bar{\xi}\cdot x}$ lies in a different Littlewood-Paley
block). By the embedding $\dot{B}^{0}_{p,2}\hookrightarrow L^{p}$ whenever
$p\in[2,\infty)$ we have that $f\in L^{p}([-\frac{1}{2},\frac{1}{2}]^{2})$ for
all $p\in[1,\infty)$.

We now cutoff $f$ so that it lies in $L^{1}(\mathbb{R}^{d})$. Let the
mollifier $\rho$ be supported on $|x|\leq R$, then $\rho_{\delta}$ is
supported on $|x|\leq\delta R$. Taking a compactly supported cutoff $\phi(x)$
with $\phi(x)=1$ on $|x|\leq\frac{3}{4}$ and $\phi(x)=0$ on $|x|\geq 1$, then
$\phi f\in L^{p}(\mathbb{R}^{2})$ for all $p\in[1,\infty)$, and
$r_{\delta}(u,\phi f)(x)=f*K_{\delta}$ provided that
$|x|\leq\frac{3}{4}-\delta R$.

Therefore,

$\int_{0}^{1}\left\|r_{\delta}(u,f\phi)\right\|_{L^{1}}^{q}\frac{d\delta}{\delta}=\infty.$

Note that $\phi f$ above is complex valued, but then for either the real, or
imaginary part, the above integral must also be infinite. ∎

###### Proof of Theorem 2.8.

Throughout $\lesssim$ will denote less than or equal to up to a constant
depending only on $p,p^{\prime}$, the dimension $d$, and the mollifier $\rho$
(and in particular not on $a,b$).

Now, as before we write

$\displaystyle r_{\delta}(u,f)(x)$
$\displaystyle=\int_{\mathbb{R}^{d}}f(x-h)(u(x)-u(x-h))\cdot\nabla\rho_{\delta}(h)\;dy$
$\displaystyle=\int_{\mathbb{R}^{d}}f(x-h)\sum^{d}_{i,j=1}\left(\int_{0}^{1}h_{j}\frac{\partial
u_{i}}{\partial x_{j}}(x-th)\;dt\right)\frac{\partial\rho_{\delta}}{\partial
h_{i}}(h)\;dh.$

Recalling that
$\rho_{\delta}(h)=\delta^{-d}\rho\left(\frac{h}{\delta}\right)$, for any
$0<a<b$ we have

$\int_{a}^{b}r_{\delta}(u,f)(x)\;\frac{d\delta}{\delta}=\int_{\mathbb{R}^{d}}f(x-h)\sum^{d}_{i,j=1}\left(\int_{0}^{1}h_{j}\frac{\partial
u_{i}}{\partial
x_{j}}(x-th)\;dt\right)\left(\int_{a}^{b}\frac{\partial\rho}{\partial
h_{i}}\left(\frac{h}{\delta}\right)\;\frac{d\delta}{\delta^{d+2}}\right)\;dh.$

Using that $\nabla\cdot u=0$ we may rewrite this as

$\displaystyle\int_{a}^{b}r_{\delta}(u,f)(x)\;\frac{d\delta}{\delta}$
$\displaystyle=\int_{\mathbb{R}^{d}}f(x-h)\sum^{d}_{i,j=1}\left(\int_{0}^{1}\frac{\partial
u_{i}}{\partial
x_{j}}(x-th)\;dt\right)\left(\int_{a}^{b}\frac{\partial(h_{j}\rho)}{\partial
h_{i}}\left(\frac{h}{\delta}\right)\;\frac{d\delta}{\delta^{d+1}}\right)\;dh$
$\displaystyle=\int_{\mathbb{R}^{d}}f(x-h)\sum^{d}_{i,j=1}\left(\int_{0}^{1}\frac{\partial
u_{i}}{\partial x_{j}}(x-th)\;dt\right)K_{i,j}(h)\;dh,$ (10)

where

$K_{i,j}(x)=\int_{a}^{b}\frac{\partial(x_{j}\rho)}{\partial
x_{i}}\left(\frac{x}{\delta}\right)\;\frac{d\delta}{\delta^{d+1}}.$

In particular, we have the following Calderon-Zygmund estimates on $K_{i,j}$

$\displaystyle\begin{aligned}
|K_{i,j}(x)|&\leq\left\|(1+|x|^{d+1})\frac{\partial(x_{j}\rho)}{\partial
x_{i}}\right\|_{L^{\infty}}\left(\int_{|x|}^{\infty}\delta^{-d-1}d\delta+\int_{0}^{|x|}|x|^{-d-1}d\delta\right)\\\
&\lesssim|x|^{-d},\end{aligned}$ $\displaystyle\begin{aligned}
\left|\frac{\partial K_{i,j}}{\partial
x_{k}}(x)\right|&=\left|\int_{a}^{b}\frac{\partial^{2}(x_{j}\rho)}{\partial
x_{k}\partial
x_{i}}\left(\frac{x}{\delta}\right)\;\frac{d\delta}{\delta^{d+2}}\right|\\\
&\leq\left\|(1+|x|^{d+2})\frac{\partial^{2}(x_{j}\rho)}{\partial x_{k}\partial
x_{i}}\right\|_{L^{\infty}}\left(\int_{|x|}^{\infty}\delta^{-d-2}d\delta+\int_{0}^{|x|}|x|^{-d-2}d\delta\right)\\\
&\lesssim|x|^{-d-1},\end{aligned}$ $\displaystyle\begin{aligned}
|\hat{K}_{i,j}(\xi)|&=\left|\int_{\mathbb{R}^{d}}e^{i\xi\cdot
x}\int_{a}^{b}\frac{\partial(x_{j}\rho)}{\partial
x_{i}}\left(\frac{x}{\delta}\right)\;\frac{d\delta}{\delta^{d+1}}\;dx\right|\\\
&=\left|\int_{a}^{b}\left(\int_{\mathbb{R}^{d}}e^{i\delta\xi\cdot
h}\frac{\partial(x_{j}\rho)}{\partial
x_{i}}(h)\;dh\right)\frac{d\delta}{\delta}\right|\\\
&=\left|\int_{a}^{b}-i\delta\xi_{i}\frac{\partial\hat{\rho}}{\partial\xi_{j}}(\delta\xi)\;\frac{d\delta}{\delta}\right|\\\
&\leq\int_{0}^{\infty}|\xi_{i}|\left|\frac{\partial\hat{\rho}}{\partial\xi_{j}}\right|(\delta\xi)\;d\delta\\\
&\leq\left\|(1+|\xi_{i}|^{2})\frac{\partial\hat{\rho}}{\partial\xi_{j}}\right\|_{L^{\infty}}\left(\int_{0}^{|\xi_{i}|^{-1}}|\xi_{i}|\;d\delta+\int_{|\xi_{i}|^{-1}}^{\infty}|\xi_{i}|^{-1}\delta^{-2}\;d\delta\right)\\\
&\lesssim 1,\end{aligned}$

where in particular none of the bounds depend on $a,b$.

The result then follows by applying Theorem 1.1 in [seeger2019multilinear] to
(10). ∎

## Acknowledgements

Lucas Huysmans acknowledges support from the UK Engineering and Physical
Sciences Research Council (EPSRC) under grant numbers EP/T517847/1 and
EP/V52024X/1. Ayman R. Said acknowledges funding from the UKRI grant SWAT.

## References


#  Thermal equilibrium states and instability of self-gravitating particles in
an asymptotically AdS spacetime

Hiroki<EMAIL_ADDRESS>Division of Particle
and Astrophysical Science, Graduate School of Science, Nagoya University,
Nagoya 464-8602, Japan  Chul-Moon<EMAIL_ADDRESS>Division of Particle and Astrophysical Science, Graduate School of Science,
Nagoya University, Nagoya 464-8602, Japan

###### Abstract

We investigate the existence and the stability of spherically symmetric
thermal equilibrium states of the self-gravitating many-particle system that
satisfies the Einstein-Vlasov equations with a negative cosmological constant.
While a thermal equilibrium state of the self-gravitating particle system
cannot have a finite mass without an artificial wall in the asymptotically
flat case, in the asymptotically AdS case, the total mass can be finite due to
the AdS potential barrier without any artificial wall. In this case,the AdS
radius characterizes the typical size of the system. The two independent
parameters parametrize the equilibrium states. Taking the total rest mass as
the unit and fixing the AdS radius, we obtain the one-parameter family of
equilibria which describes a curve in the parameter space spanned by the
gravothermal energy and the temperature. Then we investigate the instability
of the system based on the turning point method for each value of the AdS
radius. We find that the curve typically has a double spiral structure as in
the asymptotically flat case with an artificial wall. The equilibrium
solutions exist only for the finite parameter region of the gravothermal
energy between the onsets of the spirals.

##  1 Introduction

Self-gravitating systems have distinct thermodynamical properties and have
attracted much attention. One of the most interesting examples describing such
properties is the gravothermal catastrophe (instability) associated with the
negative specific heat. Antonov investigated the self-gravitating many-
particle system surrounded by a spherical adiabatic wall of the radius $R$ by
introducing a statistical description [1]. Ref. [1] showed that the Maxwell-
Jüttner distribution gives the distribution function for an equilibrium state,
and the system has no global maximum of the entropy while it may have a local
maximum. They also showed that the second variation of the entropy cannot be
negative if the ratio between the density at the center and one at the edge of
the system is smaller than 709. Lynden-Bell and Wood applied the Poincaré’s
stability criterion[2] to the self-gravitating system[3]. Poincaré’s criterion
says that a series of equilibria can destabilize only when it passes through a
turning point. The number of unstable modes changes through a turning point of
the series of equilibria [4, 5, 6, 7]. By using the turning point method,
Lynden-Bell and Wood[3] showed that no equilibrium solutions exist if the
gravothermal energy $E$ is lower than the value
$E_{\mathrm{min}}=-0.335N^{2}/R$ with $N$ being the number of particles, and
the point with this value corresponds to the critical point found by
Antonov[1]. In other words, if the radius $R$ of the artificial wall which
confines the system is not sufficiently small, the system is unstable. While
the above statement is based on the microcanonical ensemble, i.e., the
particle number and the total energy are fixed in the variation, the
thermodynamical stabilities in other processes were studied in Refs. [3, 8, 9,
10].

In the case of Newtonian gravity, the system can be thermodynamically stable
as long as the radius of the artificial wall is sufficiently small. On the
other hand, in general relativity, it is expected that the system becomes
unstable against gravitational collapse if the confined region is too small
compared with the total energy of the system. This expectation was shown to be
correct in Refs. [11, 12, 13, 14, 15, 16] for massive particle systems (see
also [17, 18] for radiation). They showed that a curve for the series of
equilibria has a double spiral structure. One of the spirals appears in the
low energy region, and the other appears in the high energy region. The former
one corresponds to the Newtonian gravothermal catastrophe, and the latter one,
which causes gravitational collapse, is peculiar to general relativity. In
this paper, we analyze the relativistic many-particle system satisfying the
Einstein-Vlasov equations without introducing any artificial wall but with a
negative cosmological constant.

Apart from the classical gravothermal catastrophe, the asymptotically AdS
gravitational system would be attractive in the context of the AdS/CFT
correspondence [19, 20, 21] and the AdS instability [22]. Because of the
confined structure and the nonlinearity of the field equations, arbitrarily
small perturbations may form black holes in an asymptotically AdS spacetime.
In Ref. [23, 24], they showed that arbitrarily small data that form black
holes exist for the spherically symmetric Einstein-massless Vlasov system. the
final states may depend on the symmetry of the spacetime and initial
conditions of perturbations [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35], and
the final fate has not been well understood yet[36, 37]. The self-gravitating
particle system considered in this paper might provide a macroscopic model of
excitations in the asymptotically AdS spacetime. Therefore the thermodynamical
stability of the system might be helpful to guess the final fate of the AdS
instability.

This paper is organized as follows. In Sec. 2, we derive basic equations and
expressions for physical quantities of an equilibrium state of the self-
gravitating many-particle system. The parameter dependence of the spatial
profile of the system is shown in Sec. 3.1. The instabilities of the
equilibria are discussed based on the turning point method for a typical case
in Sec. 3.2, and we show the parameter dependence of the configuration of the
curve describing the series of equilibria in Sec. 3.3. Sec. 4 is devoted to a
summary and conclusion.

Throughout this paper, we use the geometrized units in which both the speed of
light and Newton’s gravitational constant are unity, $c=G=1$.

##  2 Equations and physical quantities of an equilibrium state

Let us consider the relativistic system of self-gravitating massive particles,
each of which has the same rest mass. For notational simplicity, hereafter, we
use the units in which the rest mass of a particle is unity. Then the number
of particles is equivalent to the total rest mass of the system. We focus on a
spherically symmetric spacetime, whose metric is given by

$\displaystyle\differential{s^{2}}=-\mathrm{e}^{2\nu(r)}\differential{t^{2}}+\mathrm{e}^{2\mu(r)}\differential{r^{2}}+r^{2}\differential{\Omega_{2}^{2}},$
(1)

where
$\differential{\Omega_{2}^{2}}=\differential{\theta^{2}}+\sin^{2}\theta\differential{\phi^{2}}$
is a line element on the sphere. Considering the one-particle distribution
function $f(x^{\mu},p^{i})$, we can write the energy-momentum tensor of the
self-gravitating system as

$T_{\mu\nu}=\int dV_{p}p_{\mu}p_{\nu}f(x^{\mu},p^{i}),$ (2)

where $p^{\mu}$ is the particle four-momentum and $dV_{p}$ is the invariant
volume element in the momentum space. The invariant volume element in the
momentum space can be written as

$dV_{p}=2\sqrt{-g}\ d^{4}p\ \delta(p^{2}+1)\ \theta(\varepsilon)=\sqrt{-g}\
d^{3}p/\varepsilon,$ (3)

where $\delta$ and $\theta$ are the delta function and the Heaviside’s step
function, respectively, and $\varepsilon:=-p_{t}$. We note that, in the last
expression, $\varepsilon$ is a dependent variable of $p^{r}$, $p^{\theta}$ and
$p^{\phi}$ through the on-shell condition.

A thermal equilibrium state is realized by the distribution function
maximizing the Boltzmann-Gibbs entropy:

$S:=\int_{\Sigma}\differential{\Sigma_{\mu}}\ s^{\mu},\quad
s^{\mu}:=-\int\differential{V_{p}}p^{\mu}f(\log f-1),$ (4)

where $\Sigma$ is a $t={\mathrm{const.}}$ hypersurface and $s^{\mu}$ is the
entropy current. We introduce the quasi-local mass $M(r)$ as

$\displaystyle M(r):=\int_{0}^{r}\differential{s}4\pi s^{2}\rho(s)$ (5)

with $\rho=-T^{t}_{t}$. We also define the number of particles $N$ within the
radius $r$ as

$\displaystyle N(r)=\int_{0}^{r}\differential{s}4\pi s^{2}n(s),$ (6)

where $n(r):=\mathrm{e}^{\nu+\mu}n^{t}$ is the number density with the current
$n^{\mu}$ given by

$\displaystyle n^{\mu}:=\int dV_{p}p^{\mu}f(x^{\mu},p^{i}).$ (7)

Considering isolated systems, we fix the total mass $M:=\lim_{r\to\infty}M(r)$
and the number of particles $N:=\lim_{r\to\infty}N(r)$. Then, introducing the
positive numbers $\alpha$ and $\beta$ as the Lagrange multipliers, we can
obtain the distribution function for the thermal equilibrium state from the
condition called the Gibbs relation:

$\variation{S}+\alpha\variation{N}-\beta\variation{M}=0.$ (8)

As is shown in Refs. [38, 17], this condition leads to the Maxwell-Jüttner
distribution function $f=\exp(\alpha-\beta\varepsilon)$.

Since the Maxwell-Jüttner distribution is isotropic in the momentum space, the
energy momentum tensor is written in the isotropic form
$T^{\mu}_{\nu}=\mathrm{diag}(-\rho,p,p,p)$ where $p:=T_{r}^{r}$ is the
pressure. We derive the expression for each component of the energy-momentum
tensor according to the method used in the context of the Einstein-Vlasov
system[39, 40, 41]. Introducing the new integration variables $(J,\psi)$ as

$\displaystyle p^{\theta}=:\frac{\sqrt{J}}{r^{2}}\cos\psi,\quad
p^{\phi}\sin\theta=:\frac{\sqrt{J}}{r^{2}}\sin\psi,$ (9)

we can rewrite the volume element as
$dV_{p}=\frac{\mathrm{e}^{\nu+\mu}}{2r^{2}\varepsilon}dp^{r}\wedge dJ\wedge
d\psi$. The relevant regions of $p^{r}$ and $J$ are $-\infty<p^{r}<\infty$ and
$0\leq J<\infty$. Considering $\varepsilon$ as an independent variable instead
of $p^{r}$, we obtain the relevant region
$\varepsilon_{0}:=\mathrm{e}^{\nu}\leq\varepsilon$ and $0\leq J\leq J_{\rm
max}:=r^{2}(\mathrm{e}^{-2\nu}\varepsilon^{2}-1)$.

By performing the transformation $(p^{r},J)\to(\varepsilon,s)$ with
$s:=J/J_{\mathrm{max}}$, we obtain

$\displaystyle
dV_{p}=\frac{\varepsilon_{0}^{-2}}{\sqrt{1-s}}\pqty{\varepsilon^{2}-\varepsilon_{0}^{2}}^{\frac{1}{2}}d\varepsilon\wedge
ds\wedge d\psi,$ (10)

where we multiplied the factor two to take into account the degeneracy
associated with the sign of $p^{r}$. The energy density and the pressure for
the Maxwell-Jüttner distribution function are given by the integration of
$T_{t}^{t}$ and $T_{r}^{r}$ over the momentum space as

$\displaystyle\rho(r)$
$\displaystyle=4\pi\mathrm{e}^{\alpha}\bqty{\frac{K_{1}(\beta\varepsilon_{0})}{\beta\varepsilon_{0}}+\frac{3K_{2}(\beta\varepsilon_{0})}{(\beta\varepsilon_{0})^{2}}},$
(11a) $\displaystyle p(r)$
$\displaystyle=\frac{4\pi\mathrm{e}^{\alpha}K_{2}(\beta\varepsilon_{0})}{(\beta\varepsilon_{0})^{2}},$
(11b)

where $K_{1}$ and $K_{2}$ are the modified Bessel functions of the second
kind. The number density for the Maxwell-Jüttner distribution is also
calculated as

$\displaystyle
n(r)=4\pi\mathrm{e}^{\alpha+\mu}\frac{K_{2}(\beta\varepsilon_{0})}{\beta\varepsilon_{0}}\,.$
(12)

The Einstein’s equations with a negative cosmological constant
$G_{\mu\nu}+\Lambda g_{\mu\nu}=8\pi T_{\mu\nu}$ reduce to the following
equations:

$\displaystyle\mathrm{e}^{-2\mu}=1-\frac{2M(r)}{r}+\frac{r^{2}}{L^{2}},$ (13a)
$\displaystyle\nu^{\prime}=\frac{\frac{M(r)}{r^{2}}+\frac{r}{L^{2}}+4\pi
rp(r)}{1-\frac{2M(r)}{r}+\frac{r^{2}}{L^{2}}},$ (13b)

with the regularity condition at the center: $\mu(0)=0$, where
$L:=\sqrt{-3/\Lambda}$ is the AdS radius.

Let us introduce $\rho_{c}$, $p_{c}$ and $\nu_{c}$ as the values of $\rho$,
$p$ and $\nu$ at the origin. We also define the following variables at the
origin: $w_{c}=p_{c}/\rho_{c}$, $\beta_{c}:=\beta\mathrm{e}^{\nu_{c}}$,
$n_{c}=n(0)=\beta_{c}w_{c}\rho_{c}$, and $\ell:=(4\pi\rho_{c})^{-1/2}$.
Introducing $x=r/\ell$, $y:=\nu-\nu_{c},\ \tilde{\rho}=\rho/\rho_{c}$ and
$\tilde{p}=p/p_{c}$, we obtain the mass and the particle number as

$\displaystyle M(r)=\ell\int_{0}^{x}dzz^{2}\tilde{\rho}(z),\quad
N(r)=\ell\int_{0}^{x}dzz^{2}\tilde{n}(z).$ (14)

Then by defining $\tilde{M}(x)=M(r)/l$ and $\tilde{N}(x)=M(r)/l$, we can
rewrite Eq.(13b) as follows:

$\displaystyle\derivative{y}{x}=\frac{\frac{\tilde{M}(x)}{x^{2}}+\frac{x}{\lambda^{2}}+w_{c}x\tilde{p}(x)}{1-\frac{2\tilde{M}(x)}{x}+\frac{x^{2}}{\lambda^{2}}}\,.$
(15)

where $\lambda:=L/\ell$ is the normalized AdS radius. We can get a self-
consistent equilibrium solution by solving Eq.(15) together with
$d\tilde{M}/dx=x^{2}\tilde{\rho}$. The boundary conditions are given by
$y(0)=y^{\prime}(0)=\tilde{M}(0)=0$ from the definitions of $y(x)$ and
$\tilde{M}(x)$. Giving the parameter set ($\beta_{c}$,$\lambda$) and solving
Eq.(15) from the center to infinity, we obtain a unique self-consistent
solution. In this sense, the equilibrium solution is described by the two
independent parameters $\beta_{c}$ and $\lambda$, therefore it forms a two
parameter family of solutions.

Normalizing the physical quantities with respect to the particle number
$N\pqty{r}$, we obtain

$\displaystyle\tilde{E}(x):=\frac{\tilde{M}(x)-\tilde{N}(x)}{\tilde{N}(x)}$
$\displaystyle=\frac{\int_{0}^{x}\differential{z}z^{2}\tilde{\rho}(z)}{\beta_{c}w_{c}\int_{0}^{x}\differential{z}z^{2}\tilde{n}(z)}-1,$
(16a) $\displaystyle\tilde{L}(x):=\frac{\lambda}{\tilde{N}(x)}$
$\displaystyle=\frac{\lambda}{\beta_{c}w_{c}\int_{0}^{x}\differential{z}z^{2}\tilde{n}(z)},$
(16b) $\displaystyle\tilde{r}(x):=\frac{x}{\tilde{N}(x)}$
$\displaystyle=\frac{x}{\beta_{c}w_{c}\int_{0}^{x}\differential{z}z^{2}\tilde{n}(z)}\,,$
(16c)

where $\tilde{n}:=n/n_{c}$ is the normalized particle density. We note that we
can regard $\tilde{E}$ and $\tilde{L}$ at infinity as the independent
parameters to specify an equilibrium solution instead of $\beta_{c}$ and
$\lambda$ although there may be a degeneracy in the space of
$(\tilde{E}(\infty),\tilde{L}(\infty))$.

##  3 Results

###  3.1 Parameter dependence of the profile

Fig. 1 and Fig. 2 show the density profile and $M(r)/r$ as functions of $r$
for each parameter set $(\beta_{c},\lambda)$, respectively.

$x=r/\ell$

$\tilde{\rho}=\rho/\rho_{c}$

$\log\beta_{c}$

(a) $\beta_{c}$ dependence for $\lambda\to\infty$.

$x=r/\ell$

$\tilde{\rho}=\rho/\rho_{c}$

$\log\lambda$

(b) $\lambda$ dependence for $\beta_{c}=1$.

Fig. 1:  The energy density profile as a function of $x$.

$\log x=\log\frac{r}{\ell}$

$\log\frac{M(r)}{r}$

$\log\beta_{c}$

(a) $\beta_{c}$ dependence for $\lambda\to\infty$.

$\log x=\log\frac{r}{\ell}$

$\log\frac{M(r)}{r}$

$\log\lambda$

(b) $\lambda$ dependence for $\beta_{c}=1$.

Fig. 2:  The value of $M\pqty{r}/r$ as a function of $x$ in log scale. Fig.
2(a) implies that the system limits to the case of the Newtonian gravity in
the low temperature limit $\beta_{c}\to\infty$. Fig. 2(b) implies that an
effective wall exists around $x=\lambda$.

As we can see in Fig. 1 the system becomes more compact as $\lambda$
decreases. We also find $\tilde{\rho}(x)\ll 1$ for $x\gg\lambda$ in Fig. 1,
and $M\pqty{r}/r\propto 1/x$ for $x\gg\lambda$ for a finite value of $\lambda$
while $M\pqty{r}/r$ asymptotes to some finite value for the case
$\lambda=\infty$, namely $\Lambda=0$. Therefore the AdS radius $\lambda$ works
like the radius of the wall which confines the particles, and the total mass
$M$ and the particle number $N$ diverge for $\Lambda=0$ without an artificial
wall. The Newtonian limit can be realized in the low temperature limit, namely
$\beta_{c}\rightarrow\infty$, and we find the maximum value of $M\pqty{r}/r$
decreases with increasing $\beta_{c}$ as is expected in the Newtonian limit.
Note that the solutions satisfy the Buchdahl inequality $2M(r)/r<8/9$, i.e.
$\log{M/r}<log(4/9)\sim-0.81093$. This result is consistent with the results
in Ref .[42].

###  3.2 Cold and Hot Turning points

Before considering our system, let us briefly review the turning point method
and the gravothermal catastrophe. The self-gravitating system inside the
adiabatic wall of the radius $R$ has the instability called gravothermal
catastrophe, which can be classified into two types[43]. One is the “strong
instability” indicating that no equilibrium solutions exist, and the other is
the “weak instability” indicating that the entropy of the equilibrium is not a
local maximum, but just an extremum. In the asymptotically flat Newtonian
case, the series draws a one-dimensional curve parametrized by
$\sigma:=\rho_{N}(0)/\rho_{N}(R)$, i.e.
$(\bar{E},\bar{\beta})=(\bar{E}(\sigma),\bar{\beta}(\sigma))$, where the
variables $\bar{E}$ and $\bar{\beta}$ are conjugate with each other in the
sense that the variation of the entropy $\bar{S}$ is given by

$\delta\bar{S}=\bar{\beta}\delta\bar{E}.$ (17)

The stability of the system changes at the so-called turning points satisfying
$\derivative*{\bar{E}}{\sigma}=0$ [2]. The turning point corresponds to the
point at which one of the eigenvalues of the second variation of the entropy
vanishes. Therefore, the weak instability sets in at the turning point.
Futhermore, the strong instability also sets in at the turning point because
the gravothermal energy is bounded by the value at the point. It should be
noted that, when we apply the turning point method, the thermal equilibrium
condition, which we impose on the system, would be essential[4, 5]. Actually,
in the Einstein-Vlasov system without thermal equilibrium condition, there is
an example indicating the violation of the turning point principle[44].

In the relativistic system [15, 16], as is explicitly shown in Sec. 2 for the
asymptotically AdS case, the series of equilibria draws a two-dimensional
surface. In order to apply the turning point method [2] to the system, let us
consider the one-parameter family of equilibria for each value of $\tilde{L}$,
namely the isolated system with a fixed cosmological constant. Since we take
the number of particles as a unit and $\variation{N}=0$ in Eq. (8), the
variation of the entropy is given by Eq. (17) with
$(\bar{E},\bar{\beta})=(LE/N^{2},N\beta/L)\evaluated{}_{x\to\infty}$ and
$\bar{S}=S/N$.

An equilibrium state is uniquely identified by a parameter set
$(\beta_{c},\lambda)$. Then, for a given value of $\gamma:=\log\tilde{L}$, the
sequence is described by a curve on the plane spanned by $\bar{E}$ and
$\bar{\beta}$ as is shown in Fig. 3(a). We also define the sharpness parameter
$\sigma:=-\log\tilde{\rho}(\lambda)$ in analogy with the Newtonian case. Fig.
3 shows the series of equilibria for $\gamma=2.8$, and the sharpness parameter
is also plotted as a function of $\bar{E}$.

$\bar{E}=\frac{LE}{N^{2}}$

$\bar{\beta}=\frac{N\beta}{L}$

(a) $\bar{\beta}=\bar{\beta}(\bar{E})$.

$\bar{E}=\frac{LE}{N^{2}}$

$\sigma=-\log\frac{\rho(L)}{\rho_{c}}$

(b) $\sigma=\sigma(\bar{E})$.

Fig. 3:  Fig. 3(a) shows a typical curve for a series of equilibria with
compactness parameter $\gamma=2.8$. Fig. 3(b) shows the sharpness parameter
$\sigma$ as a function of $\bar{E}$. The curve has a double spiral structure
and we can see that $\sigma$ gets larger from the center to both sides and the
curve forms spiral structure after $\sigma$ exceeds a certain value.

A curve for a series of equilibria typically has a double spiral structure.
This structure in Fig. 3(a) is similar to a curve in a relativistic particle
system confined by an artificial wall[15, 16]. Having this similarity,
following Refs. [15, 16], we call the right (high energy) spiral a “hot
spiral” and the left (low energy) one a “cold spiral”. We note that, in Fig.
3(b), a curve upward to the right indicates the negative specific heat.
Therefore the specific heat, which is the second variation of the entropy
$\bar{S}$ with respect to the gravothermal energy $\bar{E}$, vanishes at the
turning point. It means that, along the line towards a spiral, another
unstable mode comes in at each turning point on the spiral. The sharpness
parameter $\sigma$ takes the minimum value at a certain point between the two
spirals, and increases towards the end of the spirals as is shown in Fig. 3.
Therefore, towards both ends of the curve, the density profile gets sharper
and the system gets more unstable.

The two spirals are characterized by $M/L$ which describes the compactness of
the total mass $M$. It is related to the gravothermal energy as
$M/L=N/L\pqty{N\bar{E}/L+1}=\mathrm{e}^{-\gamma}\pqty{1+\mathrm{e}^{-\gamma}\bar{E}}$.
Since the value $\gamma$ is fixed for each curve, the system is more compact
for a larger value of $\bar{E}$. Therefore the cold spiral is located in the
region where the system is less compact. In this case, particles can take
energy away from the central region, and the temperature in the central region
gets higher because of the negative specific heat. Then the central region
gets more compact, and the instability associated with the fragmentation is
induced. Therefore the cold spiral corresponds to the well-known Newtonian
gravothermal catastrophe. In contrast with the cold spiral, the hot spiral is
formed because the system is too compact to be stable. That is, the system
collapses if too many particles and too much energy are forced to cram into a
small region. Such a behavior can be interpreted by considering the effective
potential for a particle motion in a Schwarzschild-AdS spacetime. When the AdS
radius is sufficiently large, the effective potential has a local minimum and
a stable orbit exists. On the other hand, when the AdS radius is smaller than
a certain value, the potential minimum and the stable orbit disappear, and the
particles inevitably fall into the center of the system.

###  3.3 Parameter dependence of the spiral structure

Let us consider the $\gamma$-dependence of the configuration of the curve. If
we take $\gamma\to\infty$, namely in the dilute limit, the hot spiral moves to
the right infinitely, which corresponds to the Newtonian limit. As in the case
of the relativistic system with an artificial wall, two spirals approach to
each other as $\gamma$ gets smaller, and they merge at the “merge point”
$\gamma=\gamma_{\mathrm{m}}\simeq 2.5300$. As $\gamma$ gets even smaller, two
spirals connect to each other and form a loop at the “loop point”
$\gamma=\gamma_{\mathrm{\ell}}\simeq 2.4468$. If we continue to make $\gamma$
smaller, the loop shrinks and vanishes at the “vanishing point”
$\gamma=\gamma_{\mathrm{v}}\simeq 2.2286$. All these behaviors are shown in
Fig. 4 by plotting curves for some specific values of $\gamma$, and
qualitatively similar to the case of the relativistic system with an
artificial wall reported in Refs. [15, 16].

$\bar{E}=\frac{LE}{N^{2}}$

$\bar{\beta}=\frac{N\beta}{L}$

(a) $\gamma=4.00$.

$\bar{E}=\frac{LE}{N^{2}}$

$\bar{\beta}=\frac{N\beta}{L}$

(b) $\gamma=2.50$.

$\bar{E}=\frac{LE}{N^{2}}$

$\bar{\beta}=\frac{N\beta}{L}$

(c) $\gamma=2.44$.

$\bar{E}=\frac{LE}{N^{2}}$

$\bar{\beta}=\frac{N\beta}{L}$

(d) $\gamma=2.25$.

Fig. 4:  The series of equilibria for various $\gamma=\log\frac{L}{N}$. There
are three critical points
$\gamma_{\mathrm{m}},\gamma_{\mathrm{\ell}},\gamma_{\mathrm{v}}$ and there are
no equilibrium solutions for $\gamma<\gamma_{\mathrm{v}}$.

##  4 Conclusion

We have studied the existence and the stability of the thermal equilibrium
states of a self-gravitating many-particle system confined by an AdS potential
by applying the turning point method to the series of one-parameter families
of the equilibria. We found that all properties, which are explained below,
are similar to the case of the relativistic system confined with an artificial
wall without a cosmological constant but confined with an artificial wall [15,
16]. In our analyses, we regard the AdS radius as the size of an effective
wall.

The equilibria can be parametrized by two parameters. Fixing the rest mass
compactness $\gamma$, the set of equilibria reduces to a one-dimensional curve
in the space spanned by the normalized gravothermal energy $\bar{E}$ and the
inverse temperature $\bar{\beta}$. Each curve typically has a double spiral
structure indicating two independent instabilities. One of them, called a cold
spiral, appears at the region of the lower gravothermal energy $\bar{E}$. The
other spiral called a hot spiral appears in the region of the higher
gravothermal energy. Therefore, the equilibrium solutions exist only for the
finite parameter region of the gravothermal energy between the onsets of the
two spirals, and it indicates the strong instability. The turning points are
also onsets of weak instabilities. The cold spiral corresponds to the well-
known Newtonian gravothermal catastrophe while the hot spiral corresponds to
the instability associated with the strong gravity of the relativistic compact
object. Although the origins of the instabilities are different, both of them
are associated with the negative specific heat.

The configuration of the curve in the $(\bar{E},\bar{\beta})$ space depends on
the value of $\gamma$. In the limit $\gamma\to\infty$, the hot spiral moves
away to infinity, and there are three critical values at which the topology of
the curve changes. At the merge point $\gamma_{\mathrm{m}}\simeq 2.5300$, two
spirals merge with each other. At the loop point $\gamma_{\mathrm{\ell}}\simeq
2.4468$, the curve forms a loop and the two endpoints disappear. At the
vanishing point $\gamma_{\mathrm{v}}\simeq 2.2286$, the loop vanishes and no
equilibrium solutions exist for $\gamma<\gamma_{\mathrm{v}}$. Thus the system
cannot be a thermal equilibrium state if too many particles are confined in a
small region.

It should be noted that we have not shown that the solutions between two
spirals are stable. In the Newtonian system with an artificial adiabatic wall,
the stability of the system can be analytically proven by showing all
eigenvalues of the second variation of the entropy are negative if the
equilibrium state is not in the cold spiral. The same analysis cannot be
directly applied to our case, and the clarification of the stability is left
as a future work.

## Acknowledgments

We would like to thank T. Harada for useful discussions. This work was
supported by JSPS KAKENHI Grant Numbers JP19H01895 (C.Y.), JP20H05850 (C.Y.),
JP20H05853 (C.Y.).

## References

  * [1] V. A. ANTONOV.  Most probable phase distribution in spherical star systems and conditions for its existence.  Vest. Leningrad Univ., 7:135, 1962.
  * [2] H. PoincarÃ©.  Sur l’Ã©quilibre d’une masse fluide animÃ©e d’un mouvement de rotation.  Acta Math., 7:259–380, 1885.
  * [3] D. Lynden-Bell and Roger Wood.  The gravo-thermal catastrophe in isothermal spheres and the onset of red-giant structure for stellar systems.  mnras, 138:495, Jan 1968.
  * [4] J. Katz.  On the number of unstable modes of an equilibrium.  Mon. Not. R. Astron. Soc., 183:765–770, June 1978.
  * [5] J. Katz.  On the Number of Unstable Modes of an Equilibrium - Part Two.  Mon. Not. R. Astron. Soc., 189:817, December 1979.
  * [6] R. Sorkin.  Criterion for the onset of instability at a turning point.  Astrophys. J., 249:1, 10 1981.
  * [7] R. Sorkin.  A Stability Criterion for Many Parameter Equilibrium Families.  Astrophys. J., 257:847, June 1982.
  * [8] M. Lecar and J. Katz.  The stability of the grand microcanonical ensemble for bounded isothermal spheres.  Astrophys. J., 243:983–986, February 1981.
  * [9] Chavanis, P. H.  Gravitational instability of finite isothermal spheres.  Astron. Astrophys., 381(1):340–356, 2002.
  * [10] Chavanis, P. H.  Gravitational instability of isothermal and polytropic spheres.  Astron. Astrophys., 401(1):15–42, 2003.
  * [11] Zacharias Roupas.  Thermodynamical instabilities of perfect fluid spheres in General Relativity.  Class. Quant. Grav., 30(11):115018, 2013.  [Erratum: Class.Quant.Grav. 32, 119501 (2015)].
  * [12] Zacharias Roupas.  Relativistic Gravothermal Instabilities.  Class. Quant. Grav., 32(13):135023, 2015.
  * [13] Zacharias Roupas.  Corrigendum: Thermodynamical instabilities of perfect fluid spheres in general relativity (2013 class. quantum grav.30115018).  Classical and Quantum Gravity, 32(11):119501, may 2015.
  * [14] Zacharias Roupas.  Gravitational instability caused by the weight of heat.  Symmetry, 11(12):1435, 2019.
  * [15] Giuseppe Alberti and Pierre-Henri Chavanis.  Caloric curves of classical self-gravitating systems in general relativity.  Phys. Rev. E, 101(5):052105, 2020.
  * [16] Pierre-Henri Chavanis.  Statistical mechanics of self-gravitating systems in general relativity: Ii. the classical boltzmann gas.  The European Physical Journal Plus, 135(3):310, 2020.
  * [17] Rafael D. Sorkin, Robert M. Wald, and Zhang Zhen Jiu.  Entropy of self-gravitating radiation.  General Relativity and Gravitation, 13(12):1127–1146, 1981.
  * [18] Pierre-Henri Chavanis.  Relativistic stars with a linear equation of state: analogy with classical isothermal spheres and black holes.  Astron. Astrophys., 483:673, 2008.
  * [19] Juan Martin Maldacena.  The Large N limit of superconformal field theories and supergravity.  Int. J. Theor. Phys., 38:1113–1133, 1999.
  * [20] S.S. Gubser, I.R. Klebanov, and A.M. Polyakov.  Gauge theory correlators from non-critical string theory.  Physics Letters B, 428(1-2), May 1998.
  * [21] Edward Witten.  Anti de sitter space and holography, 1998.
  * [22] Piotr BizoÅ and Andrzej Rostworowski.  Weakly turbulent instability of antiâde sitter spacetime.  Physical Review Letters, 107(3), Jul 2011.
  * [23] Georgios Moschidis.  The Einstein–null dust system in spherical symmetry with an inner mirror: structure of the maximal development and Cauchy stability.  4 2017.
  * [24] Georgios Moschidis.  A proof of the instability of AdS for the Einstein–null dust system with an inner mirror.  Anal. Part. Diff. Eq., 13(6):1671–1754, 2020.
  * [25] Oscar J. C. Dias, Gary T. Horowitz, Don Marolf, and Jorge E. Santos.  On the Nonlinear Stability of Asymptotically Anti-de Sitter Solutions.  Class. Quant. Grav., 29:235019, 2012.
  * [26] Maciej Maliborski and Andrzej Rostworowski.  Time-Periodic Solutions in an Einstein AdS–Massless-Scalar-Field System.  Phys. Rev. Lett., 111:051102, 2013.
  * [27] Alex Buchel, Steven L. Liebling, and Luis Lehner.  Boson stars in AdS spacetime.  Phys. Rev. D, 87(12):123006, 2013.
  * [28] Venkat Balasubramanian, Alex Buchel, Stephen R. Green, Luis Lehner, and Steven L. Liebling.  Holographic Thermalization, Stability of Anti–de Sitter Space, and the Fermi-Pasta-Ulam Paradox.  Phys. Rev. Lett., 113(7):071601, 2014.
  * [29] Piotr Bizoń and Andrzej Rostworowski.  Comment on “Holographic Thermalization, Stability of Anti–de Sitter Space, and the Fermi-Pasta-Ulam Paradox”.  Phys. Rev. Lett., 115(4):049101, 2015.
  * [30] Venkat Balasubramanian, Alex Buchel, Stephen R. Green, Luis Lehner, and Steven L. Liebling.  Reply to Comment on “Holographic Thermalization, Stability of Anti–de Sitter Space, and the Fermi-Pasta-Ulam Paradox”.  Phys. Rev. Lett., 115(4):049102, 2015.
  * [31] Fotios Dimitrakopoulos and I-Sheng Yang.  Conditionally extended validity of perturbation theory: Persistence of AdS stability islands.  Phys. Rev. D, 92(8):083013, 2015.
  * [32] Stephen R. Green, Antoine Maillard, Luis Lehner, and Steven L. Liebling.  Islands of stability and recurrence times in AdS.  Phys. Rev. D, 92(8):084001, 2015.
  * [33] Ben Craps, Oleg Evnin, and Joris Vanhoof.  Renormalization group, secular term resummation and AdS (in)stability.  JHEP, 10:048, 2014.
  * [34] Ben Craps, Oleg Evnin, and Joris Vanhoof.  Renormalization, averaging, conservation laws and AdS (in)stability.  JHEP, 01:108, 2015.
  * [35] Piotr Bizoń, Maciej Maliborski, and Andrzej Rostworowski.  Resonant Dynamics and the Instability of Anti–de Sitter Spacetime.  Phys. Rev. Lett., 115(8):081103, 2015.
  * [36] Matthew W. Choptuik, Jorge E. Santos, and Benson Way.  Charting islands of stability with multioscillators in antiâde sitter space.  Physical Review Letters, 121(2), Jul 2018.
  * [37] Ramon Masachs and Benson Way.  New islands of stability with double-trace deformations.  Physical Review D, 100(10), Nov 2019.
  * [38] J. R. Ipser.  A binding-energy criterion for the dynamical stability of spherical stellar systems in general relativity.  Astrophys. J., 238:1101–1110, 1980.
  * [39] Gerhard Rein and Alan D. Rendall.  Compact support of spherically symmetric equilibria in nonrelativistic and relativistic galactic dynamics.  Math. Proc. Cambridge Phil. Soc., 128:363–380, 2000.
  * [40] Hakan Andreasson.  The Einstein-Vlasov System/Kinetic Theory.  Living Rev. Rel., 14:4, 2011.
  * [41] Håkan Andréasson, David Fajman, and Maximilian Thaller.  Static solutions to the Einstein-Vlasov system with non-vanishing cosmological constant.  SIAM J. Math. Anal., 47(4):2657–2688, 2014.
  * [42] Piyabut Burikham, Krai Cheamsawat, Tiberiu Harko, and Matthew J. Lake.  The minimum mass of a spherically symmetric object in $D$-dimensions, and its implications for the mass hierarchy problem.  Eur. Phys. J. C, 75(9):442, 2015.
  * [43] T. Padmanabhan.  Statistical mechanics of gravitating systems.  Physics Reports, 188(5):285 – 362, 1990.
  * [44] Sebastian Günther, Jacob Körner, Timo Lebeda, Bastian Pötzl, Gerhard Rein, Christopher Straub, and Jörg Weber.  A numerical stability analysis for the Einstein-Vlasov system.  Class. Quant. Grav., 38(3):035003, 2021.