$\Phi(\textup{d}x)=\phi(x)\textup{d}x+c_{0}\delta_{\rho(M_{\gamma+\lambda})}(x)$
with density $\phi\in L^{1}(\mathbb{R}_{*})$ and with $c_{0}\in\mathbb{R}$. To
exclude that the self-similar solution $\Phi$ has a Dirac in
$\rho(M_{\gamma+\lambda})$, we notice that

$\displaystyle\int_{\delta}^{\xi-\delta}K(\xi-\eta,\eta)c_{0}\delta_{\rho(M_{\gamma+\lambda})}(\eta)c_{0}\delta_{\rho(M_{\gamma+\lambda})}(\xi-\eta)\textup{d}\eta$
$\displaystyle=c_{0}^{2}\int_{\delta}^{\xi-\delta}K(\xi-\eta,\eta)\delta_{\rho(M_{\gamma+\lambda})}(\eta)\delta_{\rho(M_{\gamma+\lambda})}(\xi-\eta)\textup{d}\eta$

is non-zero only if $\xi=2\rho(M_{\gamma+\lambda})$. Hence the coagulation
operator applied to measure $c_{0}\delta_{\rho(M_{\gamma+\lambda})}$ produces
a Dirac in $2\rho(M_{\gamma+\lambda})$, contradicting the fact that
$\Phi_{|(\rho(M_{\gamma+\lambda}),\infty)}$ is absolutely continuous with
respect to the Lebesgue measure. As a consequence, we deduce that $c_{0}=0$. ∎

######  Proposition 5.15 (Moment bounds).

Assume $K$ to be a homogeneous symmetric coagulation kernel, of homogeneity
$\gamma$, satisfying (1.5), (1.6), with $\gamma,\lambda$ such that (1.7) holds
and such that $\gamma+2\lambda\geq 1.$ Let
$\Phi\in\mathcal{M}_{+}(\mathbb{R}_{*})$ be a self-similar profile as in
Definition 3.1 and assume additionally that $\Phi((0,\delta))=0$, where
$\delta>0$. Then

$\int_{\mathbb{R}_{*}}x^{\mu}\Phi(\textup{d}x)=\int_{[\delta,\infty)}x^{\mu}\Phi(\textup{d}x)<\infty,\quad\forall\mu\in\mathbb{R}.$

###### Proof.

For $\mu\leq\gamma+\lambda$ the statement follows by the fact that
$\Phi((0,\delta))=0$ and that by definition of self-similar profile we know
that

$\int_{[\delta,\infty)}x^{\gamma+\lambda}\Phi(\textup{d}x)<\infty.$

We want to prove the statement for $\mu>\gamma+\lambda.$ As a first step, we
prove that there exists a $\overline{\delta}>0$ such that

$\int_{\mathbb{R}_{*}}x^{\gamma+\lambda+\overline{\delta}}\phi(x)\textup{d}x<\infty,$
(5.45)

where $\phi$ is the density of $\Phi$. As a second step we will prove that

$\int_{\mathbb{R}_{*}}x^{\gamma+\lambda+n\overline{\delta}}\phi(x)\textup{d}x<\infty$
(5.46)

for $n\geq 1$ implies that

$\int_{\mathbb{R}_{*}}x^{\gamma+\lambda+(n+1)\overline{\delta}}\phi(x)\textup{d}x<\infty.$
(5.47)

The desired conclusion then follows by induction.

Step 1
Consider an integrable function $\varphi$ and integrate both sides of equation
(5.43) against the function $\psi(x)=\int_{x}^{\infty}\varphi(y)\textup{d}y$
to deduce that

$\displaystyle\int_{0}^{\infty}\psi(z)z\phi(z)\textup{d}z-\int_{0}^{\infty}\phi(z)\psi^{\prime}(z)\left(\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2}{1-\gamma}z^{2}\right)\textup{d}z-\frac{1}{M_{\gamma+\lambda}}\int_{0}^{\infty}\psi(z)z^{\gamma+\lambda}\phi(z)\textup{d}z$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{\infty}xK(x,y)\left(\psi(x+y)-\psi(x)\right)\phi(x)\phi(y)\textup{d}y\textup{d}x.$
(5.48)

We remark that the above formulation requires that $\psi$ has to decay fast
enough so that all the integrals in (5.4) are finite.

Select $\overline{\delta}>0$ such that
$\max\\{0,1-2(\gamma+\lambda)\\}<\overline{\delta}<1-\gamma-\lambda$ and
consider $\psi_{R}(x)=x^{\gamma+\lambda+\overline{\delta}-1}$ if $x\leq R$
while $\psi_{R}(x)=R^{\overline{\delta}}x^{\gamma+\lambda-1}$ if $x>R$. Using
(5.4) and the fact that $\gamma+2\lambda\geq 1$ and that
$\psi_{R}(x+y)-\psi_{R}(x)\leq 0$ we deduce that

$\displaystyle
0<\frac{\gamma+2\lambda-1+2\overline{\delta}}{1-\gamma}\int_{0}^{R}z^{\gamma+\lambda+\overline{\delta}}\phi(z)\textup{d}z\leq\frac{1}{M_{\gamma+\lambda}}\int_{0}^{\infty}z^{2(\gamma+\lambda)+\overline{\delta}-1}\phi(z)\textup{d}z<\infty,$

where the moment $M_{2(\gamma+\lambda)+\overline{\delta}-1}$ is bounded
because, due to the choice of $\overline{\delta}$, we have that
$2(\gamma+\lambda)+\overline{\delta}-1<\gamma+\lambda$. Passing to the limit
as $R$ tends to infinity we deduce that (5.45) holds.

Step 2
We assume that inequality (5.46) holds and we want to prove (5.47). Taking
$\psi_{R}(x)=x^{\gamma+\lambda+(n+1)\overline{\delta}-1}$ when $x\leq R$ and
$\psi_{R}(x)=R^{\overline{\delta}}x^{\gamma+\lambda+n\overline{\delta}-1}$
when $x>R$ in (5.4) we deduce that

$\displaystyle
0<\frac{\gamma+2\lambda-1+2\overline{\delta}(n+1)}{1-\gamma}\int_{0}^{R}z^{\gamma+\lambda+(n+1)\overline{\delta}}\phi(z)\textup{d}z\leq\frac{1}{M_{\gamma+\lambda}}\int_{0}^{\infty}z^{2(\gamma+\lambda)+(n+1)\overline{\delta}-1}\phi(z)\textup{d}z$
$\displaystyle+\int_{\delta}^{\infty}\int_{\delta}^{\infty}xK(x,y)\left[\psi_{R}(x+y)-\psi_{R}(x)\right]\phi(x)\phi(y)\textup{d}y\textup{d}x+\int_{0}^{\infty}\phi(z)\psi_{R}^{\prime}(z)z^{\gamma+\lambda+1}\phi(z)\textup{d}z.$

Thanks to the choice of $\overline{\delta}$, we have that
$2(\gamma+\lambda)+(n+1)\overline{\delta}-1\leq\gamma+\lambda+n\overline{\delta}$
and the desired conclusion follows for $n$ such that
$\gamma+\lambda+(n+1)\overline{\delta}-1<1$ by the fact that in that case
$\psi_{R}$ is decreasing and hence the coagulation term is negative as well as
the term
$\int_{0}^{\infty}\phi(z)\psi_{R}^{\prime}(z)z^{\gamma+\lambda+1}\phi(z)$. We
examine now the case in which $\gamma+\lambda+(n+1)\overline{\delta}-1\geq 1$.
In this case we have that

$\displaystyle\int_{\delta}^{\infty}\int_{\delta}^{\infty}xK(x,y)\left[\psi_{R}(x+y)-\psi_{R}(x)\right]\phi(x)\phi(y)\textup{d}y\textup{d}x$
$\displaystyle\leq
c\int_{\delta}^{\infty}\int_{\delta}^{\infty}xK(x,y)\left[(x+y)^{\gamma+\lambda+(n+1)\overline{\delta}-1}-x^{\gamma+\lambda+(n+1)\overline{\delta}-1}\right]\phi(x)\phi(y)\textup{d}y\textup{d}x$
$\displaystyle\leq c\iint_{\\{(x,y)\in[\delta,\infty)^{2}:x\leq
y\\}}x^{1-\lambda}y^{\gamma+\lambda}\left[(x+y)^{\gamma+\lambda+(n+1)\overline{\delta}-1}-x^{\gamma+\lambda+(n+1)\overline{\delta}-1}\right]\phi(x)\phi(y)\textup{d}y\textup{d}x$
$\displaystyle+c\iint_{\\{(x,y)\in[\delta,\infty)^{2}:x\geq
y\\}}x^{1+\gamma+\lambda}y^{-\lambda}\left[(x+y)^{\gamma+\lambda+(n+1)\overline{\delta}-1}-x^{\gamma+\lambda+(n+1)\overline{\delta}-1}\right]\phi(x)\phi(y)\textup{d}y\textup{d}x$
$\displaystyle\leq
c\int_{\delta}^{\infty}x^{1-\lambda}\phi(x)\textup{d}x\int_{\delta}^{\infty}y^{2(\gamma+\lambda)+(n+1)\overline{\delta}-1}\phi(y)\textup{d}y<\infty.$

Since we also have that

$\displaystyle\int_{0}^{\infty}\phi(z)\psi_{R}^{\prime}(z)z^{\gamma+\lambda+1}\phi(z)\textup{d}z\leq
c(n)\int_{0}^{\infty}\phi(z)z^{2(\gamma+\lambda)+(n+1)\delta-1}\phi(z)\textup{d}z<\infty$

the desired conclusion follows. ∎

######  Lemma 5.16 (Exponential bound).

Assume $K$ to be a homogeneous symmetric coagulation kernel, of homogeneity
$\gamma$, satisfying (1.5), (1.6), with $\gamma,\lambda$ such that (1.7) holds
and such that $\gamma+2\lambda\geq 1.$ Let
$\Phi\in\mathcal{M}_{+}(\mathbb{R})$ be as in Proposition 5.14 and let $\phi$
be its density. Then there exist two positive constants $L$ and $M$ such that

$\int_{M}^{\infty}\phi(z)e^{Lz}\textup{d}z<\infty.$ (5.49)

###### Proof.

Adapting the approach used in [11] we define the function $\Psi_{a}$ as

$\Psi_{a}(L):=\int_{M}^{\infty}\frac{e^{L\min\\{x,a\\}}}{\min\\{x,a\\}}x^{2}\phi(x)\textup{d}x$

where $M>\delta$. Hence, by its definition

$\Psi_{a}^{\prime}(L)=\int_{M}^{\infty}e^{L\min\\{x,a\\}}x^{2}\phi(x)\textup{d}x.$

We consider a function $\psi$ in (5.4) with
$\psi^{\prime}(x):=e^{L\min\\{x,a\\}}$ to deduce that

$\displaystyle\frac{2}{1-\gamma}\Psi_{a}^{\prime}(L)$
$\displaystyle\leq\frac{M^{\gamma+\lambda-1}}{M_{\gamma+\lambda}}\int_{M}^{\infty}e^{L\min\\{x,a\\}}x^{2}\phi(x)\textup{d}x+\frac{1}{M_{\gamma+\lambda}}\int_{\delta}^{M}e^{L\min\\{x,a\\}}x^{\gamma+\lambda+1}\phi(x)\textup{d}x$
$\displaystyle+\int_{\mathbb{R}_{*}}\int_{\mathbb{R}_{*}}xK(x,y)\int_{y}^{x+y}e^{L\min\\{w,a\\}}\textup{d}w\phi(x)\phi(y)\textup{d}x\textup{d}y$
$\displaystyle+\int_{\delta}^{\overline{z}}\left(\frac{1}{M_{\gamma+\lambda}}x^{\gamma+\lambda}-x\right)\psi(x)\phi(x)\textup{d}x,$

for
$\overline{z}:=\left(\frac{1}{M_{\gamma+\lambda}}\right)^{\frac{1}{1-\gamma-\lambda}}$,
where we are assuming without loss of generality that $\overline{z}>\delta$
and since for every $z\geq\overline{z}$ we have

$z-\frac{z^{\gamma+\lambda}}{M_{\gamma+\lambda}}>0.$ (5.50)

Thus, there exists a $c\geq 0$ such that

$\displaystyle\frac{2}{1-\gamma}\Psi_{a}^{\prime}(L)\leq\Psi_{a}^{\prime}(L)\frac{M^{\gamma+\lambda-1}}{M_{\gamma+\lambda}}+c+\int_{\mathbb{R}_{*}}\int_{\mathbb{R}_{*}}xK(x,y)\int_{y}^{x+y}e^{L\min\\{w,a\\}}\textup{d}w\phi(x)\phi(y)\textup{d}x\textup{d}y.$

As in [9] and [11] we can deduce, using Jensen’s inequality together with the
fact that $\gamma+\lambda<1$ and that $-\lambda<1$, that

$\int_{\mathbb{R}_{*}}\int_{\mathbb{R}_{*}}xK(x,y)\int_{y}^{x+y}e^{L\min\\{w,a\\}}\textup{d}w\phi(x)\phi(y)\textup{d}x\textup{d}y\leq\Psi_{a}(L)^{1-\gamma-\lambda}\Psi_{a}^{\prime}(L)^{\gamma+\lambda}.$

This, together with the fact that we can take $M$ arbitrary large, implies
that

$\displaystyle
c(\gamma,\lambda)\Psi_{a}^{\prime}(L)\leq\Psi_{a}(L)^{1-\gamma-\lambda}\Psi_{a}^{\prime}(L)^{\gamma+\lambda}+c,$
(5.51)

for a positive constant $c$ and for
$c(\gamma,\lambda)=\frac{2}{1-\gamma}-\frac{M^{\gamma+\lambda-1}}{M_{\gamma+\lambda}}>0.$

By the definition of $\Psi_{a}$ we have that

$\Psi_{a}(0)\leq\frac{2}{1-\gamma}\left(1+\frac{c_{1}}{a}\right)\rightarrow\frac{2}{1-\gamma}\text{
as }a\rightarrow\infty.$

If we prove that $\limsup_{a\rightarrow\infty}\Psi_{a}(L)<\infty$ for some
$L$, then we can conclude. Indeed we would have

$M\int_{M}^{\infty}e^{Lx}\phi(x)\textup{d}x\leq\int_{M}^{\infty}e^{Lx}x\phi(x)\textup{d}x<\infty.$

Let us prove that $\limsup_{a\rightarrow\infty}\Psi_{a}(L)<\infty$ for some
$L$. First of all notice that this is true if $\Psi_{a}^{\prime}(L)\leq 1$ as
in this case $\Psi_{a}(L)\leq L+\lim_{a\rightarrow\infty}\Psi_{a}(0)$ and the
desired conclusion follows. If instead $\Psi_{a}^{\prime}>1$, then

$c(\gamma,\lambda)\Psi_{a}^{\prime}(L)\leq\Psi_{a}(L)^{1-\gamma-\lambda}\Psi_{a}^{\prime}(L)^{\gamma+\lambda}+c\leq\alpha\Psi_{a}(L)^{1-\gamma-\lambda}\Psi_{a}^{\prime}(L)+c.$

Without loss of generality we can assume $\alpha>0$ such that
$\lim_{a\rightarrow\infty}\Psi_{a}(0)\neq\alpha$. This implies that a solution
of the ODE

$\Psi_{a}^{\prime}\left(1-\alpha\Psi_{a}^{1-\gamma-\lambda}\right)=c$

exists for small intervals around $0$ and is such that
$\lim_{a\rightarrow\infty}\Psi_{a}(L)<\infty$ for $L$ in that interval. ∎

######  Proposition 5.17 (Exponential decay).

Assume $K$ to be a homogeneous symmetric coagulation kernel, of homogeneity
$\gamma$, satisfying (1.5), (1.6), with $\gamma,\lambda$ such that (1.7) holds
and such that $\gamma+2\lambda\geq 1.$ Let
$\Phi\in\mathcal{M}_{+}(\mathbb{R})$ be as in Proposition 5.14 and let $\phi$
be its density. Then there exists a positive constant $\tilde{M}$ such that

$\limsup_{z\rightarrow\infty}\phi(z)e^{\tilde{M}z}<\infty.$ (5.52)

###### Proof.

From equation (5.43) we deduce that

$\displaystyle\left(\frac{2}{1-\gamma}z^{2}-\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}(\phi)}\right)\phi(z)\leq\int_{0}^{z}x\phi(x)\textup{d}x-\frac{1}{M_{\gamma+\lambda}(\phi)}\int_{0}^{z}x^{\gamma+\lambda}\phi(x)\textup{d}x+J_{\phi}(z).$

We now show that $J_{\phi}$ decays exponentially and that the term

$\int_{0}^{z}x\phi(x)\textup{d}x-\frac{1}{\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y}\int_{0}^{z}x^{\gamma+\lambda}\phi(x)\textup{d}x$

decays exponentially too.

First of all, we prove that there exists a constant $M_{1}>0$ such that
$J_{\phi}(z)\leq e^{-M_{1}z}$ for large $z$. To this end we notice that the
bound (5.49) implies that

$\displaystyle J_{\phi}(z)$
$\displaystyle=\int_{0}^{z}\int_{z-x}^{\infty}e^{-L(x+y)}e^{L(x+y)}xK(x,y)\phi(x)\phi(y)\textup{d}y\textup{d}x$
$\displaystyle\leq
e^{-Lz}\int_{0}^{z}\int_{z-x}^{\infty}e^{L(x+y)}xK(x,y)\phi(x)\phi(y)\textup{d}y\textup{d}x\leq
e^{-M_{1}z}.$

On the other side, for large values of $z$ we have that

$\displaystyle\int_{0}^{z}x\phi(x)\textup{d}x-\frac{1}{\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y}\int_{0}^{z}x^{\gamma+\lambda}\phi(x)\textup{d}x$
$\displaystyle=\frac{1}{\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y}\left(\int_{0}^{z}x\phi(x)\textup{d}x\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y-\int_{0}^{z}x^{\gamma+\lambda}\phi(x)\textup{d}x\right)$
$\displaystyle\leq\frac{1}{\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y}\left(\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y-\int_{0}^{z}x^{\gamma+\lambda}\phi(x)\textup{d}x\right)$
$\displaystyle\leq\frac{1}{\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y}\int_{z}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y\leq
C(\rho(M_{\gamma+\lambda}))e^{-M_{2}z},$

where $M_{2}>0$.

We deduce that for large values of $z$

$\displaystyle\left(\frac{2z^{2}}{1-\gamma}-\frac{z^{\gamma+\lambda+1}}{\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y}\right)\phi(z)\leq
c\max\\{e^{-zM_{1}},e^{-zM_{2}}\\}.$

The conclusion follows since, for every
$z>\left(\frac{1-\gamma}{2\int_{0}^{\infty}x^{\gamma+\lambda}\Phi(\textup{d}x)}\right)^{\frac{1}{1-\gamma-\lambda}}$,
we have

$\displaystyle\frac{2z^{2}}{1-\gamma}-\frac{z^{\gamma+\lambda+1}}{\int_{0}^{\infty}y^{\gamma+\lambda}\phi(y)\textup{d}y}>0.$

∎

##  6 Non-existence results and properties of the self-similar profiles

To study the non-existence of the self-similar solutions we proceed by
contradiction and start by assuming that a self-similar solution exists. To
find a contradiction we need to analyse the properties of each self-similar
profile. In the case $\gamma+2\lambda>1$ the fundamental properties that we
prove to be true are that $\Phi((0,\delta))=0$ for some positive $\delta>0$
and that $\Phi$ decays sufficiently fast for large sizes. When
$\gamma+2\lambda=1$, it is possible to prove that $\Phi$ decays fast for large
sizes, but we have not been able to prove that $\Phi((0,\delta))=0$ for some
$\delta>0$. This is the reason why we require the additional condition (3.5)
in Theorem 3.5.

###  6.1 Properties of the self-similar profiles

In this section we study the properties of each self-similar solution that do
not rely on the existence of the self-similar solution and hence do not rely
on the assumption $\gamma>-1$.

######  Theorem 6.1 (Properties of the self-similar profiles).

Let $K$ be a homogeneous symmetric coagulation kernel, of homogeneity
$\gamma$, satisfying (1.5), (1.6), with $\gamma,\lambda$ such that (1.7)
holds.

  1. 1.

If $\gamma+2\lambda>1$, then every self-similar profile $\Phi$ as in
Definition 3.1 is such that $\Phi((0,\rho(M_{\gamma+\lambda})))=0$ where
$\rho(M_{\gamma+\lambda})$ is given by (3.4). Additionally, $\Phi$ is such
that

$\int_{\mathbb{R}_{*}}e^{Lx}\Phi(\textup{d}x)<\infty$ (6.1)

for $L>0$ and it is absolutely continuous with respect to the Lebesgue
measure. Its density $\phi$ is such that

$\limsup_{x\in\mathbb{R}_{*}}\phi(x)e^{\overline{L}x}\leq c$ (6.2)

for constants $\overline{L},c>0$.

  2. 2.

If $\gamma+2\lambda=1$ and if, in addition, $\Phi((0,\delta))=0$ for some
$\delta>0$, then $\Phi([\delta,\rho(M_{\gamma+\lambda})))=0$ where
$\rho(M_{\gamma+\lambda})$ is given by (3.4). In addition, (6.1) and (6.2)
hold.

First of all we prove, see Proposition 6.2, that when $\gamma+2\lambda>1$,
each solution of equation (1.19) in the sense of Definition 3.1 is zero near
the origin. The statement of Theorem 6.1 then follows by Proposition 5.14 and
Proposition 5.17 and by the fact that when $\gamma+\lambda=1$ we are assuming
that there exists a $\delta>0$ such that $\Phi((0,\delta))=0$ .

######  Proposition 6.2 (Support of the self-similar solution).

Let $K$ be a homogeneous symmetric coagulation kernel, of homogeneity
$\gamma$, satisfying (1.5), (1.6), with $\gamma,\lambda$ such that (1.7) holds
and such that $\gamma+2\lambda>1.$ Let
$\Phi\in\mathcal{M}_{+}(\mathbb{R}_{*})$ be a self-similar profile as in
Definition 3.1. Then $\Phi((0,\overline{\xi}))=0$ for

$\overline{\xi}:=\min\left\\{\left(\frac{1}{M_{\gamma+\lambda}(\Phi)}\right)^{\frac{1}{1-\gamma-\lambda}},\left(\frac{1-\gamma}{2M_{\gamma+\lambda}(\Phi)}\right)^{\frac{1}{1-\gamma-\lambda}}\right\\}.$

###### Proof.

We start explaining the proof in a heuristic manner without entering in the
technical details that will be explained later. Equation (1.19) can be
rewritten in the following flux form

$\partial_{x}\left(J_{\Phi}(x)+\left(\frac{x^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2x^{2}}{1-\gamma}\right)\Phi(x)\right)=\left(\frac{x^{\gamma+\lambda}}{M_{\gamma+\lambda}}-x\right)\Phi(x),$
(6.3)

where $J_{\phi}$ is given by (5.44). For $x$ smaller than $\overline{\xi}$ we
have that $\frac{x^{\gamma+\lambda}}{M_{\gamma+\lambda}}-x\geq 0$ and hence
the function

$x\mapsto
J_{\Phi}(x)+\left(\frac{x^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2x^{2}}{1-\gamma}\right)\Phi(x)$

is increasing and right-continuous at $x=0$. This implies that

$\lim_{x\rightarrow
0^{+}}\left(J_{\Phi}(x)+\left(\frac{x^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2x^{2}}{1-\gamma}\right)\Phi(x)\right)=L.$

Since for every $x\leq\overline{\xi}$ we also have that
$\frac{x^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2x^{2}}{1-\gamma}\geq
0$ we deduce that $L\geq 0.$

We make a scaling argument to identify the value of $L.$ To this end we notice
that the units of the flux are given by
$[J_{\Phi}(x)]=[\Phi]^{2}[x]^{3+\gamma}$, where from now on we indicate with
$[\cdot]$ the dimensional properties of a quantity, hence

$[\Phi]=\left[x\right]^{-\frac{(3+\gamma)}{2}}.$

We deduce that

$[x^{\gamma+\lambda+1}\Phi]=[x^{\gamma+\lambda+1}][\Phi]=[x^{\frac{\gamma+2\lambda-1}{2}}].$

Since $\gamma+2\lambda>1$ this implies that as $x\rightarrow 0$ the dominant
term in equation (6.3) is $J_{\Phi}$, hence $J_{\Phi}(x)\sim L$ as
$x\rightarrow 0$. Finally, we prove that $L=0$, in agreement with the
statement proven in [10] that, when $\gamma+2\lambda\geq 1$, there are no
solutions to the constant flux equation. Integrating (6.3), we deduce that

$J_{\Phi}(x)+\left(\frac{x^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2x^{2}}{1-\gamma}\right)\Phi(x)=\int_{0}^{x}\left(\frac{z^{\gamma+\lambda}}{M_{\gamma+\lambda}}-z\right)\Phi(z)\textup{d}z.$

A detailed analysis of this ODE for $\Phi$ implies that $\Phi=0$ on the
interval
$\left(0,\left(\frac{1-\gamma}{2M_{\gamma+\lambda}}\right)^{\frac{1}{1-\gamma-\lambda}}\right)$,
see Step 4.

We now explain the proof in detail. Testing (3.3) with a function of the form
$\psi(\xi)=\xi\varphi(\xi),$ where
$\varphi\in\textup{C}^{1}_{c}(\mathbb{R}_{*}),$ we get to the following
equation:

$\displaystyle\int_{(0,\infty)}$
$\displaystyle\left(\xi-\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}\right)\varphi(\xi)\Phi(\textup{d}\xi)=\int_{(0,\infty)}\varphi^{\prime}(\xi)W[\Phi](\textup{d}\xi),$
(6.4)

where

$\displaystyle
W[\Phi](\textup{d}x)=J_{\Phi}(x)\textup{d}x+\left(\frac{x^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2x^{2}}{1-\gamma}\right)\Phi(\textup{d}x).$
(6.5)

Step 1: Regularity of $W[\Phi]$ and an integral representation formula

First of all we prove that, for every $\overline{\delta}>0$, the restriction
of the measure $W[\Phi]$ on the interval $[\overline{\delta},1]$ is absolutely
continuous with respect to the Lebesgue measure. This follows from the fact
that for every function $\varphi\in
C^{1}_{c}([\overline{\delta},\overline{\xi}])$ with $0<\overline{\delta}<1$ we
have that

$\int_{\mathbb{R}_{*}}\varphi^{\prime}(\xi)W[\Phi](\textup{d}\xi)\leq
c(\Phi,\delta,\gamma,\lambda)\|\varphi\|_{\infty};$

hence $W[\Phi]$, on the interval $[\overline{\delta},\overline{\xi}]$, has a
density $W_{\phi}\in BV([\overline{\delta},\overline{\xi}])$, see [1]. Since
$\overline{\delta}$ is arbitrary we deduce that $W[\Phi]$ has a density
$W_{\phi}\in L^{1}_{loc}((0,\overline{\xi}])$. This implies that the measure
$\Phi_{|(0,\overline{\xi}]}$ is absolutely continuous with respect to the
Lebesgue measure, its density is $\phi$.

We now prove that $W_{\phi}$ is increasing on $(0,\overline{\xi}]$. To this
end we consider a sequence of functions $\\{\varphi_{n}\\}_{n}\subset
C^{1}_{c}(\mathbb{R}_{*})$ such that
$\varphi_{n}\rightarrow\chi_{[\xi_{1},\xi_{2}]}$ pointwise, with
$0<\xi_{1},\xi_{2}<\overline{\xi}$, $\varphi_{n}^{\prime}(\xi)\geq 0$ for
$\xi\in\left(0,\frac{\xi_{1}+\xi_{2}}{2}\right)$ and such that
$\varphi_{n}^{\prime}(\xi)\leq 0$ for
$\xi\in\left(\frac{\xi_{1}+\xi_{2}}{2},\infty\right)$. Substituting
$\varphi_{n}$ in (6.4) we deduce that for every $n\geq 1$

$\displaystyle\int_{0}^{\infty}$
$\displaystyle\xi\varphi_{n}(\xi)\phi(\xi)\textup{d}\xi-\frac{1}{M_{\gamma+\lambda}}\int_{0}^{\infty}\varphi_{n}(\xi)\xi^{\gamma+\lambda}\phi(\xi)\textup{d}\xi=\int_{0}^{\infty}\varphi_{n}^{\prime}(\xi)W_{\phi}(\xi)\textup{d}\xi.$

Notice that

$\int_{0}^{\infty}\xi\varphi_{n}(\xi)\phi(\xi)\textup{d}\xi\rightarrow\int_{[\xi_{1},\xi_{2}]}\xi\phi(\xi)\textup{d}\xi\quad\text{
as }\quad n\rightarrow\infty$

and that

$\int_{0}^{\infty}\xi^{\gamma+\lambda}\varphi_{n}(\xi)\phi(\xi)\textup{d}\xi\rightarrow\int_{[\xi_{1},\xi_{2}]}\xi^{\gamma+\lambda}\phi(\xi)\textup{d}\xi\quad\text{
as }\quad n\rightarrow\infty.$

Additionally, since $W_{\phi}$ is in $L^{1}_{loc}((0,\overline{\xi}))$ and
satisfies (6.4) we deduce that $W_{\phi}\in W^{1}_{loc}(0,\overline{\xi})$ and
hence is continuous in $(0,\overline{\xi})$. Hence we have that

$\int_{(0,\infty)}\varphi_{n}^{\prime}(\xi)W_{\phi}(\xi)\textup{d}\xi\rightarrow
W_{\phi}(\xi_{1})-W_{\phi}(\xi_{2})\quad\text{ as }\quad n\rightarrow\infty.$

As a consequence, we deduce that

$\displaystyle
W_{\phi}(\xi_{2})-W_{\phi}(\xi_{1})=\int_{[\xi_{1},\xi_{2}]}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi.$
(6.6)

Since if $\xi\leq\overline{\xi}$, then
$\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\geq 0$, we deduce that
there exists a $\overline{\xi}$ such that $W_{\phi}$ is increasing in
$(0,\overline{\xi}]$. Hence, since $W_{\phi}$ is bounded, we have that the
following limit exists

$0\leq L:=\lim_{\xi\rightarrow 0^{+}}W_{\phi}(\xi)<\infty.$

Due to (3.2) we can take the limit as $\xi_{1}\rightarrow 0$ in (6.6) and
deduce that

$\displaystyle
W_{\phi}(\xi_{2})-L=\int_{(0,\xi_{2}]}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi.$
(6.7)

Step 2: we prove that $\int_{\mathbb{R}_{*}}J_{\phi}(\varepsilon
x)\varphi(x)\textup{d}x\rightarrow
L\int_{\mathbb{R}_{*}}\varphi(x)\textup{d}x$

Since $W_{\phi}$ is bounded in $(0,\overline{\xi}]$ and
$\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\frac{2}{1-\gamma}\xi\geq 0$
for every $\xi\leq\overline{\xi}$ we have that $J_{\phi}$ is also bounded in
$(0,\overline{\xi}]$. Since

$[2/3z,z]\times[2/3z,z]\subset\left\\{x\in\mathbb{R}_{*}^{2}:0<x\leq z,z-x\leq
y<\infty\right\\}=:\Omega_{z}$ (6.8)

we deduce that for every $z\in(0,\overline{\xi}]$ we have that

$z^{\gamma+1}\left(\int_{2z/3}^{z}\phi(x)\textup{d}x\right)^{2}\leq
J_{\phi}(z)\leq C$

hence

$\frac{1}{z}\int_{2z/3}^{z}\phi(x)\textup{d}x\leq\frac{C}{z^{\frac{\gamma+3}{2}}}\quad
z\in(0,\overline{\xi}].$ (6.9)

Equation (6.7) implies that

$\displaystyle\int_{\mathbb{R}_{*}}\varphi\left(x\right)J_{\phi}(\varepsilon
x)\textup{d}x$
$\displaystyle=\frac{1}{\varepsilon}\int_{\mathbb{R}_{*}}\varphi\left(\frac{x}{\varepsilon}\right)J_{\phi}(x)\textup{d}x=-\frac{1}{\varepsilon}\int_{\mathbb{R}_{*}}\varphi\left(\frac{x}{\varepsilon}\right)\left[\frac{x^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2}{1-\gamma}x^{2}\right]\phi(x)\textup{d}x$
(6.10)
$\displaystyle+\frac{1}{\varepsilon}\int_{\mathbb{R}_{*}}\varphi\left(\frac{x}{\varepsilon}\right)\int_{0}^{x}\left[\frac{z^{\gamma+\lambda}}{M_{\gamma+\lambda}}-z\right]\phi(z)\textup{d}z\textup{d}x+L\int_{\mathbb{R}_{*}}\varphi(x)\textup{d}x$

where $\varphi\in C_{c}((0,\overline{\xi}))$ and where $\varepsilon>0$.

Thanks to the bounds (3.2) for $\phi$ and to the fact that $\varphi$ is
compactly supported we can pass to the limit as $\varepsilon\rightarrow 0$ in
equation (6.10) and deduce that

$\displaystyle\lim_{\varepsilon\rightarrow
0}\int_{\mathbb{R}_{*}}J_{\phi}(\varepsilon\xi)\varphi(\xi)\textup{d}\xi=L\int_{\mathbb{R}_{*}}\varphi(x)\textup{d}x.$

for every $\varphi\in C^{1}_{c}((0,\overline{\xi}))$.

Step 3: we prove that $L=0$

We want to prove that for every $\varphi\in C^{1}_{c}((0,\overline{\xi}))$

$\lim_{\varepsilon\rightarrow
0}\int_{\mathbb{R}_{*}}J_{\phi}(\varepsilon\xi)\varphi(\xi)\textup{d}\xi=0.$

Using the formula for the fluxes we obtain that

$\displaystyle\frac{1}{\varepsilon}\int_{\mathbb{R}_{*}}J_{\phi}(\xi)\varphi\left(\frac{\xi}{\varepsilon}\right)\textup{d}\xi$
$\displaystyle=\int_{\mathbb{R}_{*}}\varphi(\xi)J_{\phi}(\varepsilon\xi)\textup{d}\xi$
$\displaystyle=\int_{\mathbb{R}_{*}}\int_{\mathbb{R}_{*}}xK(x,y)\left(\psi\left(\frac{x+y}{\varepsilon}\right)-\psi\left(\frac{x}{\varepsilon}\right)\right)\Phi(\textup{d}x)\Phi(\textup{d}y)$
$\displaystyle=I_{1}+I_{2}$

where $\psi(x)=\int_{0}^{x}\varphi(y)\textup{d}y$ and where

$I_{1}=\int_{\mathbb{R}_{*}}\int_{(y,\infty)}xK(x,y)\left(\psi\left(\frac{x+y}{\varepsilon}\right)-\psi\left(\frac{x}{\varepsilon}\right)\right)\Phi(\textup{d}x)\Phi(\textup{d}y)$

and

$I_{2}=\int_{\mathbb{R}_{*}}\int_{(0,y)}xK(x,y)\left(\psi\left(\frac{x+y}{\varepsilon}\right)-\psi\left(\frac{x}{\varepsilon}\right)\right)\Phi(\textup{d}x)\Phi(\textup{d}y).$

Since $\varphi$ is compactly supported we have, by definition, that
$\mbox{supp}\left[\psi\left(x+y\right)-\psi\left(x\right)\right]\subset\left\\{(x,y)\in\mathbb{R}_{*}^{2}:x\leq
R,\ x+y\geq\overline{\delta}>0\right\\}$ for suitable $\overline{\delta}>0$
and $R>0.$ As a consequence, using (3.2) we deduce that

$\displaystyle|I_{1}|$ $\displaystyle\leq c\int_{(0,\varepsilon
R)}\int_{(\overline{\delta}\varepsilon/2,\varepsilon
R)}y^{-\lambda}x^{1+\gamma+\lambda}\left|\psi\left(\frac{x+y}{\varepsilon}\right)-\psi\left(\frac{x}{\varepsilon}\right)\right|\Phi(\textup{d}x)\Phi(\textup{d}y)$
$\displaystyle\leq C(\psi)\frac{1}{\varepsilon}\int_{(0,\varepsilon
R)}y^{1-\lambda}\Phi(\textup{d}y)\int_{(\overline{\delta}\varepsilon/2,\varepsilon
R)}x^{1+\gamma+\lambda}\Phi(\textup{d}x)\rightarrow 0\text{ as
}\varepsilon\rightarrow 0.$

Using again (3.2) we have

$\displaystyle|I_{2}|$
$\displaystyle=\int_{\mathbb{R}_{*}}\int_{(0,y)}xK(x,y)\left|\psi\left(\frac{x+y}{\varepsilon}\right)-\psi\left(\frac{x}{\varepsilon}\right)\right|\Phi(\textup{d}x)\Phi(\textup{d}y)$
$\displaystyle\leq
C(\psi)\int_{(\varepsilon\delta/2,\infty)}y^{\gamma+\lambda}\Phi(\textup{d}y)\int_{(0,\varepsilon
R)}x^{1-\lambda}\Phi(\textup{d}x)\rightarrow 0\text{ as
}\varepsilon\rightarrow 0.$

Hence $\lim_{\varepsilon\rightarrow
0}\int_{\mathbb{R}_{*}}J_{\phi}(\varepsilon\xi)\varphi(\xi)\textup{d}\xi=0$
and therefore $L=0.$

Step 4: $\Phi$ is zero near zero

Since $L=0$ we deduce by (6.7) that

$\displaystyle
W_{\phi}(z)=\int_{0}^{z}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi\quad\text{
for }z\in(0,\overline{\xi}).$ (6.11)

Or, equivalently, since $J_{\phi}\in BV(\mathbb{R}_{*})$

$\displaystyle
J_{\phi}(z)+\left(\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2z^{2}}{1-\gamma}\right)\phi(z)=\int_{0}^{z}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi\quad\text{
for }a.e.z\in(0,\overline{\xi}).$ (6.12)

We now rewrite (6.12) in the following way

$\displaystyle\int_{0}^{z}\int_{z-x}^{z}xK(x,y)\Phi(\textup{d}y)\Phi(\textup{d}x)+\int_{0}^{z}\int_{z}^{\infty}xK(x,y)\Phi(\textup{d}y)\Phi(\textup{d}x)$
$\displaystyle+\left(\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2z^{2}}{1-\gamma}\right)\phi(z)=\int_{0}^{z}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi.$

If $\Phi$ is different to zero near the origin, then

$\displaystyle\int_{0}^{z}\int_{z}^{\infty}xK(x,y)\Phi(\textup{d}x)\Phi(\textup{d}y)$
$\displaystyle\geq
c_{1}\int_{0}^{z}x^{1-\lambda}\Phi(\textup{d}x)\int_{z}^{\infty}y^{\gamma+\lambda}\Phi(\textup{d}y)$
$\displaystyle+c_{1}\int_{0}^{z}x^{1+\gamma+\lambda}\Phi(\textup{d}x)\int_{z}^{\infty}y^{-\lambda}\Phi(\textup{d}y)\geq
C_{+}\int_{0}^{z}x^{1-\lambda}\Phi(\textup{d}x)$

where we have used that for sufficiently small $z$ we have that

$\int_{z}^{\infty}y^{\gamma+\lambda}\Phi(\textup{d}y)\geq C_{+}$

for a strictly positive constant $C_{+}$.

This implies that for $z$ small we have that

$\displaystyle 0$
$\displaystyle=J_{\phi}(z)+\left(\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2z^{2}}{1-\gamma}\right)\phi(z)-\int_{0}^{z}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi$
$\displaystyle\geq
C_{+}\int_{0}^{z}x^{1-\lambda}\phi(x)\textup{d}x+\left(\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2z^{2}}{1-\gamma}\right)\phi(z)-\int_{0}^{z}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi$
$\displaystyle\geq\left(\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2z^{2}}{1-\gamma}\right)\phi(z)+\int_{0}^{z}\left(C_{+}\xi^{1-\lambda}-\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}\right)\phi(\xi)\textup{d}\xi.$

Since we are assuming that $\Phi$ is different from zero near the origin and
since $\gamma+2\lambda>1$, we deduce that
$\int_{0}^{z}\left(C_{+}\xi^{1-\lambda}-\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}\right)\phi(\xi)\textup{d}\xi>0$.
This, together with the fact that for $z\leq\overline{\xi}$ we have
$\left(\frac{z^{\gamma+\lambda+1}}{M_{\gamma+\lambda}}-\frac{2z^{2}}{1-\gamma}\right)\geq
0$ leads to a contradiction. Hence there exists a $\delta>0$ such that
$\phi(x)=0$ for $x\in(0,\delta)$. We now prove that $\phi(x)=0$ for every
$x<\overline{\xi}$. To this end we define the function $H$ as

$H(z):=\int_{0}^{z}\left(\frac{1}{M_{\gamma+\lambda}}\xi^{\gamma+\lambda}-\xi\right)\phi(\xi)\textup{d}\xi.$

Due to (6.12) we have that $H$ satisfies

$\frac{d}{dz}H(z)=H(z)Q(z)-J_{\phi}(z)Q(z),\quad H(\delta)=0\quad\text{ for }\
a.e.\ z\in(\delta,\overline{\xi})$ (6.13)

where

$Q(z):=\frac{(1-\gamma)(z^{\gamma+\lambda}-M_{\gamma+\lambda}z)}{(1-\gamma)z^{\gamma+\lambda+1}-2M_{\gamma+\lambda}z^{2}}>0,\quad
z\in(\delta,\overline{\xi}).$ (6.14)

The solution of equation (6.13) is the function

$H(z)=-\int_{\delta}^{z}\frac{J_{\phi}(s)}{Q(s)}e^{\int_{s}^{z}Q(x)\textup{d}x}\textup{d}s.$
(6.15)

Expression (6.15) for $H$ implies that $H$ is negative for
$z\in(\delta,\overline{\xi})$. On the other hand, if we assume that $\phi$ is
different from zero in $(0,\overline{\xi})$ we deduce that $H$ is positive on
$z\in(\delta,\overline{\xi})$, by definition. This contradiction implies that
$H(z)=0$ for $z\in(\delta,\overline{\xi}).$ ∎

######  Remark 6.3.

We remark that equation (6.11) is valid both for $\gamma+2\lambda>1$ and for
$\gamma+2\lambda=1$, but the same argument we used to prove that
$\Phi((0,\delta))=0$ does not work in the case $\gamma+2\lambda=1$.

###### Proof of Theorem 6.1.

As in the proof of Theorem 5.1 we prove that if $\Phi((0,\delta))=0$ for some
$\delta>0$ then $\Phi((0,\rho(M_{\gamma+\lambda})))=0$ for
$\rho(M_{\gamma+\lambda})$ given by (3.4). To conclude the proof combine this
with Proposition 5.14 and Proposition 5.17 with Proposition 6.2. ∎

###  6.2 Non-existence results

###### Proof of Theorem 3.5.

We consider now the case $\gamma+2\lambda>1$ and $\gamma\leq-1$ in which
$\gamma+\lambda$ could be bigger than zero or smaller than zero, see Figure 1.
We will provide a proof that relies on the technical proof of Proposition 6.2.
We proceed by contradiction. Assume that a self-similar solution exists. Then
thanks to Proposition 6.2 we know that $\Phi((0,\delta))=0$ for some
$\delta>0$. Additionally, from Proposition 5.17 we know that
$\lim_{x\rightarrow\infty}\Phi(x)=0$ exponentially.

We can therefore consider the test function $\varphi\equiv 1$ in equation
(3.3). This implies

$\displaystyle
0\leq-\frac{1+\gamma}{1-\gamma}\int_{\mathbb{R}_{*}}\Phi(\textup{d}x)=-\frac{1}{2}\int_{\mathbb{R}_{*}}\int_{\mathbb{R}_{*}}K(x,y)\Phi(\textup{d}x)\Phi(\textup{d}y)<0,$

which is a contradiction.

Assume now that $\gamma+2\lambda=1$ and that $\gamma\leq-1$. Notice that in
this case we necessarily have that $\gamma+\lambda\leq 0$. Let us assume by
contradiction that a self-similar solution exists and is such that

$\int_{(0,1]}x^{-\lambda}\Phi(\textup{d}x)<\infty.$

We show that

$\int_{\mathbb{R}_{*}}\Phi(\textup{d}x)<\infty.$

Notice that this is true by assumption when $\gamma+\lambda=0$, therefore we
restrict the attention to the case $\gamma+\lambda<0.$ To this end we adapt
the argument used in the proof of Proposition 5.15, we refer there for the
technical aspects, and prove that

$\int_{\mathbb{R}_{*}}x^{\gamma+\lambda}\Phi(\textup{d}x)<\infty\quad\text{
implies
}\quad\int_{\mathbb{R}_{*}}x^{\gamma+\lambda+\overline{\delta}}\Phi(\textup{d}x)<\infty$

when $0<-\gamma-\lambda<\overline{\delta}<1-\gamma-\lambda$. To see this we
consider $\psi(x)=x^{\gamma+\lambda+\overline{\delta}-1}$ in equation (5.4) to
deduce that

$\displaystyle\frac{2\overline{\delta}}{1-\gamma}M_{\gamma+\lambda+\overline{\delta}}\leq\frac{M_{2(\gamma+\lambda)+\overline{\delta}-1}}{M_{\gamma+\lambda}}.$

Now we notice that our choice of $\overline{\delta}$ implies that

$-\lambda\leq 2(\gamma+\lambda)+\overline{\delta}-1\leq\gamma+\lambda.$

Hence $M_{2(\gamma+\lambda)+\overline{\delta}-1}<\infty$ and the desired
conclusion follows.

Similarly, we can prove that

$\int_{\mathbb{R}_{*}}x^{\gamma+\lambda+n\overline{\delta}}\Phi(\textup{d}x)<\infty\quad\text{
implies
}\quad\int_{\mathbb{R}_{*}}x^{\gamma+\lambda+(n+1)\overline{\delta}}\Phi(\textup{d}x)<\infty.$

As a consequence, taking $n$ sufficiently large, we deduce that
$M_{0}(\Phi)<\infty$. We can therefore consider the test function
$\varphi\equiv 1$ in equation (3.3). This implies

$\displaystyle
0\leq-\frac{1+\gamma}{1-\gamma}\int_{\mathbb{R}_{*}}\Phi(\textup{d}x)=-\frac{1}{2}\int_{\mathbb{R}_{*}}\int_{\mathbb{R}_{*}}K(x,y)\Phi(\textup{d}x)\Phi(\textup{d}y)<0$

which is a contradiction. ∎

#### Acknowledgements

The authors are grateful to A. Nota for many interesting mathematical
discussions that led to several of the main ideas contained in this paper. The
authors would also like to thank B. Kepka for his help clarifying some
technical points of this paper. The authors gratefully acknowledge the
financial support of the Hausdorff Research Institute for Mathematics (Bonn),
through the Junior Trimester Program on Kinetic Theory, of the CRC 1060,
Project-ID 211504053, The mathematics of emergent effects at the University of
Bonn funded through the German Science Foundation (DFG), of the Bonn
International Graduate School of Mathematics at the Hausdorff Center for
Mathematics (EXC 2047/1, Project-ID 390685813) funded through the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation), as well as of the
Atmospheric Mathematics (AtMath) collaboration of the Faculty of Science of
University of Helsinki, of the ERC Advanced Grant 741487, and of the Academy
of Finland, via the Finnish centre of excellence in Randomness and STructures
(project No. 346306). The funders had no role in study design, analysis,
decision to publish, or preparation of the manuscript.

## References

  * [1] L. Ambrosio, N. Fusco, and D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, Clarendon Press, 2000.
  * [2] R. C. Ball, C. Connaughton, P. P. Jones, R. Rajesh, and O. Zaboronski, Collective oscillations in irreversible coagulation driven by monomer inputs and large-cluster outputs, Physical Review Letters, 109 (2012), p. 168304.
  * [3] J. Banasiak, W. Lamb, and P. Laurençot, Analytic Methods for Coagulation-Fragmentation Models, Volume II, CRC Press, 2019.
  * [4] F. P. da Costa, A finite-dimensional dynamical model for gelation in coagulation processes, Journal of NonLinear Science, 8 (1998), pp. 619–653.
  * [5] S. C. Davies, J. R. King, and J. A. D. Wattis, The Smoluchowski coagulation equations with continuous injection, Journal of Physics A: Mathematical and General, 32 (1999), pp. 7745–7763.
  * [6] P. B. Dubovskiĭ, Mathematical Theory of Coagulation, Lecture Notes Series No. 23, Seoul National University, Seoul, 1994.
  * [7] M. Escobedo and S. Mischler, Dust and self-similarity for the Smoluchowski’s coagulation equation, Annales de l’Institut Henri Poincaré C, Analyse Non Linéaire, 23 (2006), pp. 331–362.
  * [8] M. Escobedo, S. Mischler, and M. R. Ricard, On self-similarity and stationary problem for fragmentation and coagulation models, Annales de l’Institut Henri Poincaré C, Analyse Non Linéaire, 22 (2005), pp. 99–125.
  * [9] M. A. Ferreira, E. Franco, and J. J. L. Velázquez, On the self-similar behaviour of coagulation systems with injection, To appear in Annales de l’Institut Henri Poincaré C, Analyse Non Linéaire, (2022).
  * [10] M. A. Ferreira, J. Lukkarinen, A. Nota, and J. J. L. Velázquez, Stationary non-equilibrium solutions for coagulation systems, Archive for Rational Mechanics and Analysis, 240 (2021), pp. 809–875.
  * [11] N. Fournier and P. Laurençot, Local properties of self-similar solutions to Smoluchowski’s coagulation equation with sum kernels, Proceedings of the Royal Society of Edinburgh-A-Mathematics, 136 (2006), pp. 485–508.
  * [12] S. K. Friedlander, Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics, Oxford University Press, New York, 2000.
  * [13] H. Hayakawa, Irreversible kinetic coagulations in the presence of a source, Journal of Physics A: Mathematical and General, 20 (1987), pp. L801–L805.
  * [14] P. L. Krapivsky and C. Connaughton, Driven Brownian coagulation of polymers, Journal Chemical Physics, 136 (2012), p. 204901.
  * [15] P. L. Krapivsky, J. F. F. Mendes, and S. Redner, Logarithmic islanding in submonolayer epitaxial growth, The European Physical Journal B-Condensed Matter and Complex Systems, 4 (1998), pp. 401–404.
  * [16] , Influence of island diffusion on submonolayer epitaxial growth, Physical Review B, 59 (1999), p. 15950.
  * [17] G. Menon and R. L. Pego, Approach to self-similarity in Smoluchowski’s coagulation equations, Communications on Pure and Applied Mathematics, 57 (2004), pp. 1197–1232.
  * [18] T. Olenius, O. Kupiainen-Määttä, I. K. Ortega, T. Kurtén, and H. Vehkamäki, Free energy barrier in the growth of sulfuric acid-ammonia and sulfuric acid-dimethylamine clusters, The Journal of Chemical Physics, 139 (2013), p. 084312.
  * [19] M. Smoluchowski, Drei Vorträge über Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteilchen, Zeitschrift fur Physik, 17 (1916), pp. 557–585.
  * [20] H. V. Tran and T.-S. Van, Coagulation-fragmentation equations with multiplicative coagulation kernel and constant fragmentation kernel, Communications on Pure and Applied Mathematics, 75 (2022), pp. 1292–1331.
  * [21] H. Vehkamäki and I. Riipinen, Thermodynamics and kinetics of atmospheric aerosol particle formation and growth, Chemical Society Reviews, 41 (2012), pp. 5160–73.
  * [22] V. E. Zakharov, V. S. L’vov, and G. Falkovich, Kolmogorov Spectra of Turbulence I, Springer Series in Nonlinear Dynamics, 1992.


# Excitation of electrostatic solitary waves and surface waves in ion beam
neutralization process

Nakul Nuwal#<EMAIL_ADDRESS>Igor D. Kaganovich† Deborah A. Levin#
#University of Illinois Urbana-Champaign, IL, USA, †Princeton Plasma Physics
Lab, NJ, USA

###### Abstract

Unusually long electrostatic solitary waves (ESWs) are discovered in 2D and 3D
Particle-in-Cell studies of the process of ion beam neutralization by electron
emission from filaments. These ESWs are long because trapped and untrapped
electron density perturbations nearly compensate each other. Surface waves
were discovered in the process of neutralization but were only observed in 3D
simulations. This is because the phase velocity of surface waves in a 2D
geometry is higher than in 3D giving the high-energy electrons generated
upstream near the electron source enough energy to excite such waves only in
3D cylindrical beams.

Ion beams are used in various applications such as particle accelerators, ion-
thrustersWang2012 ; jambunathan2018chaos ; nuwal2020kinetic ,
nanopantographychen2019nearly , and ion-implantationconrad1989plasma . The
space-charge in the beam has to be neutralized to generate a focused ion beam,
which is best done by a background plasmadudin2009simultaneous ;
kaganovich2001nonlinear ; berdanier2015intense . However, the complexity
incurred by the generation of background plasma inside of a high-vacuum
accelerator can be avoided by using an external filament electron source,
which is easy to install in experiments. The neutralization by a filament
electron source, however, leads to the excitation of space-charge waves, such
as electrostatic solitary waves (ESWs), which limit the beam’s space-charge
neutralizationlan2020neutralization2 . While the ESWs have been studied
theoretically for many decadesng2005bernstein ; chen2003width ;
chen2004bernstein , their behavior has gained recent interest since the
discovery of their ubiquitous formation in the process of charge
neutralization of ion beams lan2019electrostatic ; lan2020neutralization2 .
The ESWs, which we also observed in our numerical simulations, have mostly
been studied in 1D in the pastbernstein1957exact ; hutchinson2017electron ;
chen2003width . In this letter, however, using our heterogeneous high-
performance-computing capable PIC code CHAOSjambunathan2018chaos , we studied
the formation of these ESWs in 2D planar and 3D cylindrical beams. We show
that the previously developed 1D theoretical models can describe the long ESWs
that form in 2D planar beams. However, they cannot describe the oblong
axisymmetric ESWs that form in 3D cylindrical beams (see Fig. 13 of Ref.
nuwal2021PRE ). Further, by conducting a theoretical analysis of these modes
similar to the original BGK treatmentbernstein1957exact , we show that because
of a non-Maxwellian electron distribution function (EVDF) formed by
neutralizing electrons in the beam, _the ESWs are much longer_ than those
observed in the previous theoretical chen2002bgk ; goldman2007theory  and
experimental workslefebvre2010laboratory .

Figure 1: Electric potential profiles along the beam axis at $t=4.72\ \mu s$,
when all cases have achieved at least $80\%$ of neutralization. Cases: 2D
(blue), $a=2.5$ mm 3D (green), $a=5.0$ mm 3D (red), and $a=7.5$ mm 3D (black);
$a$ is the beam radius. Wave features: ESWs (square boxes), electron-holes
(triangles), and surface waves (circles).

Additionally, in 3D simulations, we observed the excitation of the
Trivelpiece-Gould (TG) surface wavestrivelpiece1959space  in the process of
space charge neutralization of an ion beam. The TG waves have been
experimentally studied for plasma columns in waveguides lynov1979observations
; moisan1982experimental  where they were excited by an externally applied
field. In this work, however, an ion beam together with neutralizing electrons
form a plasma column that can serve as a plasma waveguide for surface waves
moisan1982experimental ; i.e., emitted electrons from a filament form an
electron beam at early stages of neutralization process that is uniquely
suited to spontaneously excite such surface waves with a high phase speed. At
later time in the neutralization process, these high energy electrons leave
the simulation’s domain due to wall losses. This prevents the decay of TG
waves due to the Landau damping by fast electrons, which allows TG waves to
exist for very long times. Finally, we derived a dispersion relation for the
surface waves in a 2D planar ion beam in a metal chamber and identified that
the phase speed of the surface waves in 2D planar ion beams is sufficiently
different than 3D, explaining why surface waves are only observed in 3D.

Particle-in-Cell (PIC) simulations were performed in a 2D domain of $6\times
40$ $\rm cm^{3}$ and a 3D domain of $3\times 3\times 40$ $\rm cm^{3}$ (see
Figs. 2(a) and 2(b) of Ref. nuwal2021PRE ), where a monoenergetic argon ion
beam of 38 keV was injected at the center of the $z=0$ mm plane from a
rectangular(circular) region for 2D(3D) case(s). A uniform grid with cell
sizes ranging from $0.3\lambda_{\rm D}$ to $0.91\lambda_{\rm D}$ was used for
both 2D and 3D cases, here, $\lambda_{\rm D}$ is the local Debye length at the
end of neutralization process when the electron density is highest. The 40 cm
length of our domain and about $2.5$ mm beam inlet radius (half-width for 2D)
is used to be comparable to the previously conducted plasma neutralization
experimentsstepanov2016dynamics . Note that this study is also relevant to
nanopantography applicationschen2019nearly . The electrons with a $T_{e}=2$ eV
Maxwellian EVDF and current, $I_{e}=I_{\rm Ar^{+}}/3$, $I_{\rm Ar^{+}}$ being
the ion beam current, were introduced at the beam axis at $z=10$ cm. The
neutralization process of the ion beam happens in two phases: (1) rapid
neutralization, which occurs when the beam potential, $\phi>>k_{B}T_{e}/e$,
and (2) slow neutralization when $\phi\approx k_{B}T_{e}/e$. After the initial
rapid neutralization phase, the rate of neutralization slows down considerably
because the remaining low unneutralized potential can no longer trap all
emitted electrons and high-energy electrons escape to the
wallslan2020neutralization .

More than $80\%$ of neutralization is achieved by $t=4.0\ \mu s$. Figure 1
shows the electric potential profiles at the end of the first phase in the
neutralization process, in which ESWs form both upstream and downstream of the
electron source for the 2D and 3D beam cases with $a=2.5$ mm, whereas ESWs
only form in the upstream and surface waves in the downstream regions of the
electron source for 3D cases with $a=5.0$ and 7.5 mm, as shown in the figure
inset. Space-charge sheaths form at $z_{\rm min}$ and $z_{\rm max}$ boundaries
and reflect most electrons and produce two electron streams formed inside the
ion beam. The two-stream instability between two electron streams results in
the formation of ESWs, as was reported in Refs.lan2020neutralization2 ;
omura1996electron ; roberts1967nonlinear . For the neutralization process, the
ESWs were observed to be very robust: they are ubiquitously excited during the
neutralization process, survive many collisions between themselves and
reflection from the boundaries.lan2020neutralization2

These ESWs can be easily recognized as vortex-like structures in the electron
phase space, as shown in Fig. 2. These ESWs (their positions shown by the red-
dots in Fig. 2 for discrete number of times for 2D planar beam case) move,
collide, and merge with each other to eventually form new ESWs that are of
larger lengths than the original ones. Among the multiple ESWs in Fig. 2, the
ESW shown by the yellow ‘X’ is of particular interest because it moves with
nearly a constant speed, and does not collide or merge with other ESWs for a
long time. Because its electric potential amplitude was found to be similar to
the other ESWs in the system, the yellow ‘X’ ESW is representative of other
ESWs in the 2D planar beam case and is further analyzed below.

Figure 2:  (a) Electron distribution in phase space, $f_{e-}$ normalized by
$n_{0}=1.75\times 10^{14}\ {\rm m^{-3}}$, at $t=3.12\ {\rm\mu s}$. The black
solid line in (a) shows electric potential on the right vertical-axis. (b)
Time evolution of ESWs positions, shown by red circles, along the beam axis,
$(x,y)=(0.18,3.0)$ cm, for the 2D planar beam case. In (b), ESWs shown by blue
and yellow ‘X’s eventually merge with the large electron-hole near the
electron source at $z=10$ cm, and ESWs shown by white and green ‘X’s merge
with each other.

Because of their collisions and merging, ESW lengths in both 2D and 3D beams
increase to become between 30-50 local Debye lengths long, which is unusually
long when compared to the other studieslefebvre2010laboratory ; chen2002bgk ;
hutchinson2017electron . Since the beam radius is small compared with the ESW
length, we analyze the yellow ‘X’ ESW at $t=3.12\ \mu s$ in Fig. 3 using 1D
BGK theorybernstein1957exact , by fitting an analytical electric potential
profile to the PIC resultschamel1979theory ,

$\displaystyle\Delta\phi(\Delta z)={\phi_{0}}{\rm
sech^{4}}\left(\frac{b_{s}\Delta z}{\lambda_{\rm D_{0}}}\right),$ (1)

where we obtain $\phi_{0}=0.876$ V and $b_{s}=0.05$ as the fit parameters,
$\lambda_{\rm D_{0}}=0.78$ mm is the Debye length at $n_{e}=1.75\times
10^{14}\ {\rm m^{-3}}$ and $T_{e}=2$ eV. The comparison of the fit with PIC is
shown in Fig. 3, where it can be seen that it compares well for $\Delta z<10$
mm. In Eq. 1, $\Delta\phi=\phi(\Delta z)-\phi_{\rm base}$, where $\phi_{\rm
base}$ is the potential outside the ESW at $\Delta z=15$ mm. In Fig. 3(b), the
contour lines of the constant total energy, $E_{\rm tot}$, are shown by dashed
black lines and are superimposed on color-plot of the EVDF of the trapped-in-
the wave electrons, $f_{t}$, in the phase space for the ESW marked by yellow
‘X’ in Fig. 2. They indicate a small variation in electron populations at
different $E_{\rm tot}$ values, which is qualitatively similar to the
‘shallow’ ESWs discussed in Ref. hutchinson2017electron . Given that the EVDF
of electrons far from the ESW, $f_{\rm FF}$, is known, the untrapped EVDF
reads,

$\displaystyle f_{ut}\left(v_{z},\Delta z\right)=f_{\rm
FF}\left(\text{sgn}{(v_{z})}\sqrt{v_{z}^{2}-2\frac{e}{m_{e}}\Delta\phi(\Delta
z)}\right)$ (2)

where $f_{ut}(v_{z},\Delta z)$ is the EVDF of untrapped electrons at location
$\Delta z$. We can determine the electron density profiles of both the trapped
and untrapped electrons as a function of potential using Eq. 2 for untrapped
and a similar equation for trapped electrons where the EVDF is taken from the
PIC simulations results at the location of the potential maximum. We now
compare the effect of different $f_{\rm FF}$ on ESW: 1. a Maxwellian EVDF with
$T_{e}=1$ eV chosen to fit to the PIC result, and 2. non-Maxwellian EVDF taken
directly from the PIC simulation results in Fig. 3(b) at $\Delta z=15$ mm
(also see Fig. 10(c) of Ref. nuwal2021PRE  for further details).

Figure 3:  In (a), the potential profile and the trapped electron distribution
in phase space, $f_{t}$ (normalized by $n_{0}=1.75\times 10^{14}\ {\rm
m^{-3}}$), are shown for ESW marked with yellow ‘X’ in Fig. 2. (b): Dashed
black lines show constant $E_{\rm tot}$ (eV) contour lines with increments of
0.2 eV. (c) shows the electron density profile normalized by $1.5\times
10^{14}\ {\rm m^{-3}}$. (d) Comparison of predicted total (untrapped) electron
density based on Eq. 2 using a 1. Maxwellian (blue) and 2. non-Maxwellian
(red) $f_{\rm FF}$ and 3. the PIC result (green).

The ESWs generated in our PIC simulations seem to have formed having a very
small electron density perturbation of the order of 3-5%, as shown by
depletion in the normalized electron density in Fig. 3(c). This small density
perturbation in ESWs is allowed due to the non-Maxwellian background EVDF in
the beam. This is demonstrated in Fig. 3(d) where the total electron density
perturbation from Eq. 2 applied to a non-Maxwellian $f_{\rm FF}$ (‘2’) is
closer to the PIC result (‘3’) than a Maxwellian $f_{\rm FF}$ (‘1’). This
difference in the magnitude of these perturbations is entirely due to a lower
electron density decline of untrapped electrons towards the center of the ESW
for a non-Maxwellian $f_{\rm FF}$ than a Maxwellian one, as shown on the
secondary axis of Fig. 3(d). Since the EVDF at the center and the electric
potential profile are taken from the PIC result, trapped electron profiles are
identical for both $f_{\rm FF}$ cases. This small untrapped density decline
nearly compensates the trapped density in the non-Maxwellian case, resulting
in a very small total density perturbation. This shows that a Maxwellian
background EVDF would likely result in a larger electron density and electric
potential perturbation that a non-Maxwellian one found in ion beams and this
is the reason why such long ESWs with a small density and potential
perturbations appear in our ion beam. Present 1D BGK theory needs to be
revised to be able to describe ESWs in 3D due to their highly complex electron
trajectories (they are shown in Fig. 16 of Ref. nuwal2021PRE ).

Turning to understanding the importance of 3D effects, ESWs were not found
downstream of the electron source for the 3D simulations with beam radii
$a=5.0$ and $a=7.5$ mm. In these cases, the potential of the un-neutralized
ion beam is higher than in the case with $a=2.5$ mm, because there are more
ions in beams of thicker radii. Correspondingly, in these cases, electrons
gain velocities much higher than thinner beams. Our PIC simulations indicate
that, as a result, the energetic electrons excite high-phase-speed surface
waves, as shown in Fig. 4. These waves have a long wavelength
$\lambda_{w}=6.1$ cm and high phase velocity, $v_{\phi}=2.42\times 10^{6}$
m/s, for the $a=5.0$ mm case. Figure 4 shows the temporal evolution of the
EVDF for the $a=5.0$ mm case where a small number of electrons gain velocity
close to $2.42\times 10^{6}$ m/s at about $t=2.40\ \mu s$ and excite the
surface wave of that phase velocity.

We compared the surface waves parameters obtained from the simulations and the
theory using the dispersion relation for axisymmetric surface waves derived by
Trivelpiece and Gould (TG)trivelpiece1959space , as shown in Fig. 4. This
shows that axisymmetric TG surface waves become excited for 3D cases only with
$a=5.0$ and $7.5$ mm. Although high-velocity electrons excite these waves,
with the decrease in electric potential, the ion beam space charge’s ability
to generate high-energy electrons declineslan2020neutralization  with time.
This leads to a decline in the electron temperature inside the beam, shown by
the narrowing of EVDF with time in Fig. 4, which leaves the previously excited
TG waves in the domain with no high-velocity electrons in resonance with the
waves that could quench them through Landau damping. This also explains the
absence of TG waves in the $a=2.5$ mm 3D case where the ion beam space-charge
generates electrons with high-enough velocities to dampen the low phase speed
surface waves that may exist for such small radius beam(see Fig. 18(b) of Ref.
nuwal2021PRE ). For more details see Sec. 4 of Ref. nuwal2021PRE .

Figure 4: (a) Time evolution of electric potential profile $\Delta t=t-4.72\
\mu s$. (b) Evolution of EVDF, $f_{e-}$ normalized by $n_{0}=1.75\times
10^{14}\ {\rm m^{-3}}$, for 3D beam case of $a=5.0$ mm with a few electrons
reaching the surface wave phase velocity, $v_{\phi}=2.42\times 10^{6}$ m/s at
$t=2.40\ \mu s$. (c) Dispersion curves for the surface waves in 3D
cylindricaltrivelpiece1959space  and 2D planar beams.

Such surface waves, however, were not observed in the 2D planar
beamlan2020neutralization . Following a procedure similar to
Vedenovvedenov1965solid , Krallkrall1973principles , and Trivelpiece-
Gouldtrivelpiece1959space , we derived the dispersion relation of a 2D planar
beam,

$\displaystyle\frac{\omega}{\omega_{pe}}=\frac{1}{\left[1+\frac{\tanh{[(b-a)\kappa]}}{\tanh{a\kappa}}\right]^{1/2}},$
(3)

where $a$ and $b$ are the half-widths of plasma beam and channel, respectively
(See Appendix C of Ref. nuwal2021PRE  for details). Equation 3 converges to
$1/\sqrt{2}$ for $a\to\infty$, $b\to\infty$, and $a<<b$, as expected for a
planar surface wave in a vacuum with no conducting wallsvedenov1965solid . As
seen from Fig. 4, the surface waves in the 2D planar beam have a higher phase
speed, $v_{\phi}=\omega/\kappa$, than the TG surface waves for cylindrical
beams. This means that, when comparing 2D and 3D surface waves dispersion for
$a=2.5$ mm, electrons with about 2.5 times higher velocities than 3D beam
would be required to excite these waves in a 2D beam. Since such high velocity
electrons are not present when the electrons neutralize the ion beam in the 2D
simulation, the surface waves were not observed in that case.

Finally, in our PIC simulations, we were able to artificially excite these
surface waves when we introduced a small amount of emitted electrons with
velocities higher than the electron thermal velocity of $T_{e}=2$ eV in 2D and
3D beam cases of $a=2.5$ mm, where the surface waves were not spontaneously
excited originally (see Fig. 21 of Ref. nuwal2021PRE ). This shows that a
small amount of high velocity electrons may excite a surface wave and affirms
that the TG waves in our other 3D cases were excited by high energy electrons.

In summary, in our kinetic PIC modeling we observed a curious fact that ESWs
with axial lengths of 30-50 Debye lengths can be generated during
neutralization of an ion beam by electron emission. In a 2D geometry, such
long ESWs can be approximately described by a modified 1D BGK theory because
the ESW length is much larger than the beam width. The large length of ESW is
determined to be a consequence of a surprising near compensation of density
variation along the ESW by trapped and untrapped electrons because of the non-
Maxwellian background EVDF. We also observed the generation of TG surface
waves in cylindrical ion beams during neutralization process simulations in 3D
geometry; whereas these waves were not observed in 2D geometry. Analysis shows
that this can be explained by the fact that surface waves have a much higher
phase velocity in 2D than in 3D.

This research was supported by the US Department of Energy Award DE-
SC00021348. I.D.K. was supported by US Department of Energy contract DE-
AC02-09CH11466.

## References

  * (1) Wang, J., Chang, O., and Cao, Y., “Electron–ion coupling in mesothermal plasma beam emission: Full particle PIC simulations,” IEEE Transactions on Plasma Science, Vol. 40, No. 2, 2012, pp. 230–236.
  * (2) Jambunathan, R. and Levin, D. A., “CHAOS: An octree-based PIC-DSMC code for modeling of electron kinetic properties in a plasma plume using MPI-CUDA parallelization,” Journal of Computational Physics, Vol. 373, 2018, pp. 571–604.
  * (3) Nuwal, N., Jambunathan, R., and Levin, D. A., “Kinetic Modeling of Spacecraft Surfaces in a Plume Backflow Region,” IEEE Transactions on Plasma Science, Vol. 48, No. 12, 2020, pp. 4305–4325.
  * (4) Chen, Y.-M., Sawadichai, R., Tian, S., Donnelly, V. M., Ruchhoeft, P., and Economou, D. J., “Nearly monoenergetic positive ion beam with self-neutralized space charge extracted from a pulsed plasma,” Journal of Physics D: Applied Physics, Vol. 52, No. 35, 2019, pp. 355205.
  * (5) Conrad, J., “Plasma source ion implantation: A new approach to ion beam modification of materials,” Materials Science and Engineering: A, Vol. 116, 1989, pp. 197–203.
  * (6) Dudin, S. and Rafalskyi, D., “On the simultaneous extraction of positive ions and electrons from single-grid ICP source,” EPL (Europhysics Letters), Vol. 88, No. 5, 2009, pp. 55002.
  * (7) Kaganovich, I. D., Shvets, G., Startsev, E., and Davidson, R. C., “Nonlinear charge and current neutralization of an ion beam pulse in a pre-formed plasma,” Physics of Plasmas, Vol. 8, No. 9, 2001, pp. 4180–4192.
  * (8) Berdanier, W., Roy, P. K., and Kaganovich, I., “Intense ion beam neutralization using underdense background plasma,” Physics of Plasmas, Vol. 22, No. 1, 2015, pp. 013104.
  * (9) Lan, C. and Kaganovich, I. D., “Neutralization of ion beam by electron injection: Excitation and propagation of electrostatic solitary waves,” Physics of Plasmas, Vol. 27, No. 4, 2020, pp. 043104.
  * (10) Ng, C. and Bhattacharjee, A., “Bernstein-Greene-Kruskal modes in a three-dimensional plasma,” Physical review letters, Vol. 95, No. 24, 2005, pp. 245004.
  * (11) Chen, L.-J., Thouless, D. J., and Tang, J.-M., “Width-amplitude relation of Bernstein-Greene-Kruskal solitary waves,” arXiv preprint physics/0303021, 2003.
  * (12) Chen, L.-J., Thouless, D. J., and Tang, J.-M., “Bernstein–Greene–Kruskal solitary waves in three-dimensional magnetized plasma,” Physical Review E, Vol. 69, No. 5, 2004, pp. 055401.
  * (13) Lan, C. and Kaganovich, I., “Electrostatic solitary waves in ion beam neutralization,” Physics of Plasmas, Vol. 26, No. 5, 2019, pp. 050704.
  * (14) Bernstein, I. B., Greene, J. M., and Kruskal, M. D., “Exact nonlinear plasma oscillations,” Physical Review, Vol. 108, No. 3, 1957, pp. 546.
  * (15) Hutchinson, I. H., “Electron holes in phase space: What they are and why they matter,” Physics of Plasmas, Vol. 24, No. 5, 2017, pp. 055601.
  * (16) Nuwal, N., Levin, D. A., and Kaganovich, I. D., “Kinetic modeling of three-dimensional electrostatic-solitary and surface waves in beam neutralization,” Companion paper on arXiv, 2022.
  * (17) Chen, L.-J. and Parks, G., “BGK electron solitary waves: 1D and 3D,” Nonlinear Processes in Geophysics, Vol. 9, No. 2, 2002, pp. 111–119.
  * (18) Goldman, M. V., Newman, D. L., and Mangeney, A., “Theory of weak bipolar fields and electron holes with applications to space plasmas,” Physical review letters, Vol. 99, No. 14, 2007, pp. 145002.
  * (19) Lefebvre, B., Chen, L.-J., Gekelman, W., Kintner, P., Pickett, J., Pribyl, P., Vincena, S., Chiang, F., and Judy, J., “Laboratory measurements of electrostatic solitary structures generated by beam injection,” Physical review letters, Vol. 105, No. 11, 2010, pp. 115001.
  * (20) Trivelpiece, A. and Gould, R., “Space charge waves in cylindrical plasma columns,” Journal of Applied Physics, Vol. 30, No. 11, 1959, pp. 1784–1793.
  * (21) Lynov, J.-P., Michelsen, P., Pécseli, H., Rasmussen, J. J., Saeki, K., and Turikov, V., “Observations of solitary structures in a magnetized, plasma loaded waveguide,” Physica Scripta, Vol. 20, No. 3-4, 1979, pp. 328.
  * (22) Moisan, M., Shivarova, A., and Trivelpiece, A., “Experimental investigations of the propagation of surface waves along a plasma column,” Plasma physics, Vol. 24, No. 11, 1982, pp. 1331.
  * (23) Stepanov, A. D., Gilson, E. P., Grisham, L. R., Kaganovich, I. D., and Davidson, R. C., “Dynamics of ion beam charge neutralization by ferroelectric plasma sources,” Physics of Plasmas, Vol. 23, No. 4, 2016, pp. 043113.
  * (24) Lan, C. and Kaganovich, I. D., “Neutralization of ion beam by electron injection: Accumulation of cold electrons,” Physics of Plasmas, Vol. 27, No. 4, 2020, pp. 043108.
  * (25) Omura, Y., Matsumoto, H., Miyake, T., and Kojima, H., “Electron beam instabilities as generation mechanism of electrostatic solitary waves in the magnetotail,” Journal of Geophysical Research: Space Physics, Vol. 101, No. A2, 1996, pp. 2685–2697.
  * (26) Roberts, K. and Berk, H. L., “Nonlinear evolution of a two-stream instability,” Physical Review Letters, Vol. 19, No. 6, 1967, pp. 297.
  * (27) Schamel, H., “Theory of electron holes,” Physica Scripta, Vol. 20, No. 3-4, 1979, pp. 336.
  * (28) Vedenov, A., “Solid state plasma,” Soviet Physics Uspekhi, Vol. 7, No. 6, 1965, pp. 809.
  * (29) Krall, N. A. and Trivelpiece, A. W., Principles of plasma physics, McGraw-Hill, New York, 1973.


††institutetext: Department of Physics, Shinshu University,
3-1-1 Asahi, Matsumoto 390-8621, Japan

# Capacity of entanglement in random pure state

Kazumi Okuyama<EMAIL_ADDRESS>
###### Abstract

We compute the capacity of entanglement in the bipartite random pure state
model using the replica method. We find the exact expression of the capacity
of entanglement which is valid for a finite dimension of the Hilbert space. We
argue that in the gravitational path integral, the capacity of entanglement
receives contributions only from the sub-leading saddle points corresponding
to the partially connected geometries.

##  1 Introduction

In recent papers Penington:2019kki ; Almheiri:2019qdq , the Page curve of the
Hawking radiation Page:1993wv  is reproduced from the replica computation of
the entanglement entropy (see also Almheiri:2020cfm  for a review). As argued
by Page in Page:1993df , the entropy computation of the Hawking radiation is
nicely modeled by the random pure state $|\Psi\rangle$ in a bipartite Hilbert
space

$\displaystyle\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}.$ (1)

Here the subsystems $A$ and $B$ correspond to the Hawing radiation and the
black hole, respectively. From the reduced density matrix $\rho_{A}$ of the
subsystem $A$

$\displaystyle\rho_{A}=\operatorname{Tr}_{B}|\Psi\rangle\langle\Psi|,$ (2)

we can compute the entanglement entropy $S_{A}$

$\displaystyle S_{A}=-\langle\operatorname{Tr}\rho_{A}\log\rho_{A}\rangle,$
(3)

where the bracket $\langle\cdots\rangle$ denotes the ensemble average over the
random pure state $|\Psi\rangle$. The exact form of $S_{A}$ (42) is obtained
in Page:1993df  as a function of the dimensions $d_{A},d_{B}$ of the Hilbert
spaces (1)

$\displaystyle d_{A}=\dim\mathcal{H}_{A},\quad d_{B}=\dim\mathcal{H}_{B}.$ (4)

As discussed in Page:1993wv , $S_{A}$ in the random pure state model exhibits
a similar behavior as the Page curve for the Hawking radiation from an
evaporating black hole.

In a recent paper Kawabata:2021hac , it is argued that the capacity of
entanglement $C_{A}$ introduced in qi  is a useful quantity to diagnose the
phase transition around the Page time. $C_{A}$ is defined by

$\displaystyle
C_{A}=\operatorname{Tr}(\rho_{A}K^{2})-(\operatorname{Tr}\rho_{A}K)^{2},$ (5)

where $K=-\log\rho_{A}$ is the modular Hamiltonian. In other words, $C_{A}$
measures the fluctuation of the modular Hamiltonian. See e.g.
Nakaguchi:2016zqi ; Nakagawa:2017wis ; deBoer:2018mzv  for the study of
$C_{A}$ in various models.

In this paper, we compute the capacity of entanglement in the random pure
state model using the replica method. From the ensemble average of
$\operatorname{Tr}\rho_{A}^{n}$ over the random pure state, we can compute the
entanglement entropy $S_{A}$ and the capacity of entanglement $C_{A}$ as the
derivative with respect to the replica number $n$ at $n=1$

$\displaystyle S_{A}$
$\displaystyle=-\partial_{n}\log\langle\operatorname{Tr}\rho_{A}^{n}\rangle\Big{|}_{n=1}=-\partial_{n}\langle\operatorname{Tr}\rho_{A}^{n}\rangle\Big{|}_{n=1},$
(6) $\displaystyle C_{A}$
$\displaystyle=\partial_{n}^{2}\log\langle\operatorname{Tr}\rho_{A}^{n}\rangle\Big{|}_{n=1}=\partial_{n}^{2}\langle\operatorname{Tr}\rho_{A}^{n}\rangle\Big{|}_{n=1}-(S_{A})^{2}.$

We find the exact form of $C_{A}$ as a function of the dimensions
$d_{A},d_{B}$ of the Hilbert spaces $\mathcal{H}_{A},\mathcal{H}_{B}$. The
exact expression of $C_{A}$ in (46) is the main result of this paper.

As discussed in Penington:2019kki ; Almheiri:2019qdq , fully disconnected or
fully connected geometries dominate in the replica computation of the entropy
$S_{A}$, and their contributions exchange dominance around the Page time. In
the case of the capacity $C_{A}$, it turns out that the leading saddle point
from the fully disconnected or fully connected geometries does not contribute
to $C_{A}$, and it receives non-zero contributions only from the sub-leading
saddle points corresponding to the partially connected geometries. This is
consistent with the result of deBoer:2018mzv  that $C_{A}$ is a measure of the
partial entanglement.

This paper is organized as follows. In section 2, we review the random pure
state model and the known exact result of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$. We find a new formula of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ in terms of the Narayana number
(19), which is useful for the replica computation of the entanglement entropy
$S_{A}$ and the capacity of entanglement $C_{A}$. In section 3, we compute
$S_{A}$ and $C_{A}$ in the planar limit using the replica method. Our
computation shows that the leading saddle point does not contributes to
$C_{A}$ and it receives contributions only from sub-leading saddle points. In
section 4, we compute the exact $S_{A}$ and $C_{A}$ using the replica method.
Finally we conclude in section 5.

##  2 Random pure state model

In this section, let us briefly review the random pure state model. We
consider a pure state $|\Psi\rangle$ in the bipartite Hilbert space
$\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. This models the black
hole evaporation where $A$ corresponds to the Hawking radiation while $B$
corresponds to the black hole. In the model discussed in Penington:2019kki ,
$A$ and $B$ correspond to the end of the world brane and the bulk JT gravity,
respectively. We can expand the state $|\Psi\rangle$ in terms of the
orthonormal basis of $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$

$\displaystyle|\Psi\rangle=\mathcal{N}\sum_{i=1}^{d_{A}}\sum_{\alpha=1}^{d_{B}}X_{i\alpha}|i\rangle_{A}\otimes|\alpha\rangle_{B},$
(7)

where $\mathcal{N}$ is the normalization factor to ensure the unit norm of
$|\Psi\rangle$

$\displaystyle\langle\Psi|\Psi\rangle=1.$ (8)

It is useful to regard the coefficient $X_{i\alpha}$ in (7) as a component of
the $d_{A}\times d_{B}$ complex matrix $X$

$\displaystyle X=(X_{i\alpha}).$ (9)

Then the normalization factor $\mathcal{N}$ is written as

$\displaystyle\mathcal{N}=\frac{1}{\sqrt{\operatorname{Tr}(XX^{\dagger})}}.$
(10)

We are interested in the reduced density matrix $\rho_{A}$ defined in (2)
obtained by tracing out $B$. In terms of the matrix $X$ in (9), $\rho_{A}$ is
written as a $d_{A}\times d_{A}$ matrix

$\displaystyle\rho_{A}=\frac{XX^{\dagger}}{\operatorname{Tr}(XX^{\dagger})}=\frac{W}{\operatorname{Tr}W},$
(11)

where $W=XX^{\dagger}$. The ensemble average over the random pure state
$|\Psi\rangle$ can be defined by the Gaussian integral over the matrix $X$

$\displaystyle\langle\mathcal{O}(W)\rangle=\frac{\int
dXdX^{\dagger}\mathcal{O}(W)e^{-\operatorname{Tr}(XX^{\dagger})}}{\int
dXdX^{\dagger}e^{-\operatorname{Tr}(XX^{\dagger})}}.$ (12)

As a distribution of the matrix $W=XX^{\dagger}$, this is known as the
Wishart-Laguerre ensemble. See Collins2015  for a nice review on this subject.

The matrix integral (12) can be written as the eigenvalue integral by
diagonalizing the matrix $W$. In the original paper by Page Page:1993df , the
entropy $S_{A}$ was computed by evaluating the eigenvalue integral of
$\log\rho_{A}$ directly. In this paper, we will compute the entropy $S_{A}$
and the capacity $C_{A}$ using the replica method (6). To do this, we need the
expectation value of the moment $\operatorname{Tr}\rho_{A}^{n}$. Fortunately,
the exact result of $\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ is already
obtained in Nechita2007 111 As discussed in Nechita2007 , the moments of
$\rho_{A}=\frac{W}{\operatorname{Tr}W}$ and those of $W$ for the Wishart-
Laguerre ensemble are related by
$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle=\frac{\Gamma(d_{A}d_{B})}{\Gamma(d_{A}d_{B}+n)}\langle\operatorname{Tr}W^{n}\rangle.$
(13) The exact result of $\langle\operatorname{Tr}W^{n}\rangle$ in the
Wishart-Laguerre ensemble was obtained in hanlon1992some ; Haagerup2003 ,
which is equivalent to (14) via the relation (13).

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle$
$\displaystyle=\frac{\Gamma(d_{A}d_{B})}{\Gamma(d_{A}d_{B}+n)}\cdot\frac{1}{n}\sum_{j=1}^{n}(-1)^{j-1}\frac{[d_{A}+n-j]_{n}[d_{B}+n-j]_{n}}{(j-1)!(n-j)!},$
(14)

where $[a]_{n}=a(a-1)\cdots(a-n+1)$ denotes the falling factorial. For
instance, the first few terms of $\langle\operatorname{Tr}\rho_{A}^{n}\rangle$
read

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{2}\rangle$
$\displaystyle=\frac{d_{A}+d_{B}}{d_{A}d_{B}+1},$ (15)
$\displaystyle\langle\operatorname{Tr}\rho_{A}^{3}\rangle$
$\displaystyle=\frac{d_{A}^{2}+3d_{A}d_{B}+d_{B}^{2}+1}{(d_{A}d_{B}+1)(d_{A}d_{B}+2)},$
$\displaystyle\langle\operatorname{Tr}\rho_{A}^{4}\rangle$
$\displaystyle=\frac{d_{A}^{3}+6d_{A}^{2}d_{B}+6d_{A}d_{B}^{2}+d_{B}^{3}+5d_{A}+5d_{B}}{(d_{A}d_{B}+1)(d_{A}d_{B}+2)(d_{A}d_{B}+3)},$

which agree with the known exact results of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ Lubkin1978 ; Sommers2004 . Note
that $\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ in (14) is symmetric under
the exchange of $d_{A}$ and $d_{B}$, which implies that $S_{A}$ and $C_{A}$
are also symmetric functions of $d_{A}$ and $d_{B}$

$\displaystyle S_{A}(d_{A},d_{B})=S_{A}(d_{B},d_{A}),\quad
C_{A}(d_{A},d_{B})=C_{A}(d_{B},d_{A}).$ (16)

In what follows, we will assume $d_{A}\leq d_{B}$ without loss of generality.
$S_{A}$ and $C_{A}$ in the opposite regime $d_{A}>d_{B}$ can be obtained by
exchanging $d_{A}$ and $d_{B}$ using the symmetry (16).

When $n$ is a positive integer, the summation of $j$ in (14) can be extended
to $j=\infty$ since the summand vanishes for $j>n$. Then (14) becomes

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle$
$\displaystyle=\frac{\Gamma(d_{A}d_{B})}{\Gamma(d_{A}d_{B}+n)}\cdot\frac{1}{n}\sum_{j=1}^{\infty}(-1)^{j-1}\frac{[d_{A}+n-j]_{n}[d_{B}+n-j]_{n}}{(j-1)!(n-j)!}$
(17)
$\displaystyle=\frac{\Gamma(d_{A}+n)\Gamma(d_{B}+n)\Gamma(d_{A}d_{B}+1)}{\Gamma(n+1)\Gamma(d_{A}+1)\Gamma(d_{B}+1)\Gamma(d_{A}d_{B}+n)}$
$\displaystyle~{}~{}~{}\times{}_{3}F_{2}\bigl{(}\\{1-d_{A},1-d_{B},1-n\\},\\{1-d_{A}-n,1-d_{B}-n\\};1\bigr{)}.$

The last expression makes sense for non-integer $n$ and it defines an analytic
continuation of $\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ away from the
integer $n$. One can in principle compute the derivative of the last
expression in (17) with respect to $n$ to find $S_{A}$ and $C_{A}$. However,
it is not straightforward to simplify the derivative of the hypergeometric
function ${}_{3}F_{2}$ in (17).

It turns out that it is useful to rewrite (17) as the following form using the
identity of ${}_{3}F_{2}$ 222 See the identity in the Wolfram Functions Site
http://functions.wolfram.com/07.27.17.0046.01.

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle$
$\displaystyle=\frac{\Gamma(d_{B}+n)\Gamma(d_{A}d_{B}+1)}{\Gamma(d_{B}+1)\Gamma(d_{A}d_{B}+n)}$
(18)
$\displaystyle\times{}_{3}F_{2}\bigl{(}\\{1-d_{A},1-n,-n\\},\\{2,1-d_{B}-n\\};1\bigr{)}.$

We find that this is expanded as

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle=\sum_{k=1}^{\infty}N_{n,k}\frac{\Gamma(d_{A})\Gamma(d_{B}+1+n-k)\Gamma(d_{A}d_{B}+1)}{\Gamma(d_{A}+1-k)\Gamma(d_{B}+1)\Gamma(d_{A}d_{B}+n)},$
(19)

where $N_{n,k}$ is the Narayana number

$\displaystyle
N_{n,k}=\frac{1}{n}\binom{n}{k}\binom{n}{k-1}=\frac{\Gamma(n)\Gamma(n+1)}{k!(k-1)!\Gamma(1+n-k)\Gamma(2+n-k)}.$
(20)

Indeed, one can show that the summation in (19) reproduces the hypergeometric
function in (18). This expression (19) makes contact with the planar limit of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ where the Narayana number
naturally appears from the number of non-crossing permutations Kudler-
Flam:2021rpr . When $n$ is a positive integer, the summation of $k$ in (19) is
truncated at $k=n$ and one can easily check that (19) reproduces the result
(15) for small $n$.

Using the analytic continuation of the Narayana number by the last expression
in (20), we can define a natural analytic continuation of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ in (19) for non-integer $n$.
When $d_{A}\leq d_{B}$ and $d_{A}$ is a positive integer, the summand in (19)
vanishes for $k>d_{A}$ and hence (19) becomes

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle=\sum_{k=1}^{d_{A}}N_{n,k}\frac{\Gamma(d_{A})\Gamma(d_{B}+1+n-k)\Gamma(d_{A}d_{B}+1)}{\Gamma(d_{A}+1-k)\Gamma(d_{B}+1)\Gamma(d_{A}d_{B}+n)}.$
(21)

In section 4, we will use this expression of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ for the replica computation of
the exact $S_{A}$ and $C_{A}$.

##  3 Planar limit

Before discussing the exact result of $C_{A}$, in this section we will compute
$C_{A}$ in the planar limit

$\displaystyle
d_{A},d_{B}\to\infty~{}~{}\text{with}~{}~{}\alpha=\frac{d_{A}}{d_{B}}:~{}\text{fixed}.$
(22)

We will assume $\alpha\leq 1$ without loss of generality. $C_{A}$ in the
opposite regime $\alpha>1$ can be obtained by sending $\alpha\to\alpha^{-1}$
using the symmetry (16). The computation of $C_{A}$ in this limit (22) has
been already done in Kawabata:2021hac  using the planar eigenvalue density of
the Wishart-Laguerre ensemble, known as the Marchenko–Pastur distribution.
Here we will use the replica method to compute $C_{A}$, which clarifies the
role of replica wormholes in $C_{A}$.

In the planar limit (22), the exact result of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ in (21) reduces to Kudler-
Flam:2021rpr

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}=d_{A}^{1-n}\sum_{k=1}^{\infty}N_{n,k}\alpha^{k-1}.$
(23)

When $n$ is a positive integer, the sum over $k$ is truncated to $1\leq k\leq
n$ since the Narayana number $N_{n,k}$ in (20) vanishes for $k\geq n+1$. The
first few terms of the planar expectation values of
$\operatorname{Tr}\rho_{A}^{n}$ are given by

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{2}\rangle_{\text{planar}}$
$\displaystyle=\frac{d_{A}+d_{B}}{d_{A}d_{B}},$ (24)
$\displaystyle\langle\operatorname{Tr}\rho_{A}^{3}\rangle_{\text{planar}}$
$\displaystyle=\frac{d_{A}^{2}+3d_{A}d_{B}+d_{B}^{2}}{(d_{A}d_{B})^{2}},$
$\displaystyle\langle\operatorname{Tr}\rho_{A}^{4}\rangle_{\text{planar}}$
$\displaystyle=\frac{d_{A}^{3}+6d_{A}^{2}d_{B}+6d_{A}d_{B}^{2}+d_{B}^{3}}{(d_{A}d_{B})^{3}}.$

One can see that (24) is obtained from the planar limit of the exact result
(15), as expected.

When $n$ is a positive integer, (23) is expanded as

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}=d_{A}^{1-n}+\frac{1}{2}n(n-1)d_{A}^{2-n}d_{B}^{-1}+\cdots+d_{B}^{1-n}.$
(25)

In the gravitational path integral, the first term $d_{A}^{1-n}$ and the last
term $d_{B}^{1-n}$ come from fully disconnected and fully connected
geometries, respectively. If we assume $d_{A}<d_{B}$, the dominant
contribution is the first term $d_{A}^{1-n}$ and
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}$ is written as

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}=d_{A}^{1-n}\Bigl{[}1+f(n,\alpha)\Bigr{]},$
(26)

where $f(n,\alpha)$ is given by

$\displaystyle
f(n,\alpha)=\sum_{k=2}^{\infty}N_{n,k}\alpha^{k-1}={}_{2}F_{1}(1-n,-n;2;\alpha)-1.$
(27)

In other words, in the gravitational picture $f(n,\alpha)$ summarizes all
contributions from the sub-dominant, partially connected geometries. Note that
$f(n,\alpha)$ vanishes at $n=1$ by our definition in (27)

$\displaystyle f(1,\alpha)=0.$ (28)

Plugging (26) into (6) we find

$\displaystyle S_{A,\text{planar}}$
$\displaystyle=\log{d_{A}}-f^{\prime}(1,\alpha),$ (29) $\displaystyle
C_{A,\text{planar}}$
$\displaystyle=f^{\prime\prime}(1,\alpha)-f^{\prime}(1,\alpha)^{2},$

where the prime in $f^{\prime}$ and $f^{\prime\prime}$ denotes the derivative
with respect to $n$. Note that the capacity is completely determined by the
sub-leading contributions $f(n,\alpha)$ in (26). In fact, if we use the
leading approximation of the trace

$\displaystyle\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}\approx
d_{A}^{1-n},$ (30)

the capacity vanishes

$\displaystyle C_{A,\text{planar}}$
$\displaystyle=\partial_{n}^{2}\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}\Big{|}_{n=1}-\left(\partial_{n}\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}\Big{|}_{n=1}\right)^{2}$
(31) $\displaystyle\approx(\log d_{A})^{2}-(\log d_{A})^{2}=0.$

The same conclusion holds in the opposite regime $d_{A}>d_{B}$ as well if we
use the leading approximation
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle_{\text{planar}}\approx
d_{B}^{1-n}$. This implies that the dominant saddle point of gravitational
path integral does not contribute to the capacity. In other words, the
capacity of entanglement is sensitive to the sub-dominant saddle points
corresponding to the partially connected geometries. Thus, $C_{A}$ is a useful
probe of the contributions of replica wormholes which are not fully connected
nor fully disconnected geometries, but some “intermediate” geometries. This is
consistent with the result in deBoer:2018mzv  that $C_{A}$ takes a non-zero
value for partially entangled states and $C_{A}$ vanishes for the pure state
or a maximally entangled state. Namely, $C_{A}$ is a measure of partial
entanglement deBoer:2018mzv .

Let us evaluate $f^{\prime}(1,\alpha)$ and $f^{\prime\prime}(1,\alpha)$. From
(27) they are written as

$\displaystyle f^{\prime}(1,\alpha)$
$\displaystyle=\sum_{k=2}^{\infty}\partial_{n}N_{n,k}\Big{|}_{n=1}\alpha^{k-1},$
(32) $\displaystyle f^{\prime\prime}(1,\alpha)$
$\displaystyle=\sum_{k=2}^{\infty}\partial_{n}^{2}N_{n,k}\Big{|}_{n=1}\alpha^{k-1}.$

Thus, we need to compute the derivative of Narayana number $N_{n,k}$ at $n=1$.
From (20), one can easily show that $N_{n,k}$ is expanded around $n=1$ as

$\displaystyle N_{n,k}=\left\\{\begin{aligned} &1,&\quad(k=1),\\\
&\frac{1}{2}(n-1)+\frac{1}{2}(n-1)^{2},&\quad(k=2),\\\
&-\frac{(n-1)^{2}}{k(k-2)(k-1)^{2}}+\mathcal{O}\big{(}(n-1)^{3}\big{)},&\quad(k\geq
3).\end{aligned}\right.$ (33)

This implies that the first derivative $\partial_{n}N_{n,k}$ at $n=1$ vanishes
unless $k=2$, and $f^{\prime}(1,\alpha)$ in (32) becomes

$\displaystyle
f^{\prime}(1,\alpha)=\alpha\partial_{n}N_{n,2}\Big{|}_{n=1}=\frac{\alpha}{2}.$
(34)

In a similar manner, from (33) we find that $f^{\prime\prime}(1,\alpha)$ in
(32) becomes

$\displaystyle f^{\prime\prime}(1,\alpha)$
$\displaystyle=\alpha-\sum_{k=3}^{\infty}\frac{2}{k(k-2)(k-1)^{2}}\alpha^{k-1}$
(35)
$\displaystyle=-1-\frac{3\alpha}{2}+(\alpha-\alpha^{-1})\log(1-\alpha)+2\text{Li}_{2}(\alpha).$

Finally, plugging the result of $f^{\prime}(1,\alpha)$ in (34) and
$f^{\prime\prime}(1,\alpha)$ in (35) into the definition of
$S_{A,\text{planar}}$ and $C_{A,\text{planar}}$ in (29) we find

$\displaystyle S_{A,\text{planar}}$ $\displaystyle=\log
d_{A}-\frac{\alpha}{2},$ (36) $\displaystyle C_{A,\text{planar}}$
$\displaystyle=-1-\frac{3\alpha}{2}-\frac{\alpha^{2}}{4}+(\alpha-\alpha^{-1})\log(1-\alpha)+2\text{Li}_{2}(\alpha).$

This agrees with the result in Kawabata:2021hac  obtained from the
Marchenko–Pastur distribution. Our replica computation reveals the importance
of the sub-leading contribution $f(n,\alpha)$ to the capacity of entanglement.

We note in passing that $C_{A,\text{planar}}$ in (36) takes the maximal value
at $\alpha=1$, or $d_{A}=d_{B}$ Kawabata:2021hac ; deBoer:2018mzv

$\displaystyle
C_{A,\text{planar}}^{(\text{max})}=C_{A,\text{planar}}\Big{|}_{\alpha=1}=\frac{\pi^{2}}{3}-\frac{11}{4}.$
(37)

##  4 Exact capacity of entanglement at finite $d_{A},d_{B}$

In this section we will compute the exact $S_{A}$ and $C_{A}$ using the exact
result of $\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ in (21).333 The exact
computation of $S_{A}$ by the replica method is also considered in
dyer2014divergence  using the expression of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ in (14).  Here we assume that
$d_{A}$ and $d_{B}$ are both integers and $d_{A}\leq d_{B}$.

Let us first compute the entanglement entropy $S_{A}$. Plugging (21) into the
definition of $S_{A}$ in (6) we find

$\displaystyle S_{A}$ $\displaystyle=-\lim_{n\to
1}\sum_{k=1}^{d_{A}}\frac{\Gamma(d_{A})\Gamma(d_{B}+1+n-k)\Gamma(d_{A}d_{B}+1)}{\Gamma(d_{A}+1-k)\Gamma(d_{B}+1)\Gamma(d_{A}d_{B}+n)}$
(38)
$\displaystyle\times\Biggl{[}N_{n,k}\Bigl{(}\psi(d_{B}+1+n-k)-\psi(d_{A}d_{B}+n)\Bigr{)}+\partial_{n}N_{n,k}\Biggr{]},$

where $\psi(z)$ denotes the digamma function

$\displaystyle\psi(z)=\frac{d}{dz}\log\Gamma(z).$ (39)

From the behavior (33) of the Narayana number $N_{n,k}$ near $n=1$, (38)
becomes

$\displaystyle
S_{A}=\psi(d_{A}d_{B}+1)-\psi(d_{B}+1)-\frac{1}{2}\frac{\Gamma(d_{A})\Gamma(d_{B})}{\Gamma(d_{A}-1)\Gamma(d_{B}+1)}.$
(40)

Using the property of the digamma function

$\displaystyle\psi(m+1)=\sum_{k=1}^{m}\frac{1}{k}-\gamma,\quad(m\in\mathbb{N}),$
(41)

with $\gamma$ being the Euler’s constant, we arrive at the exact entanglement
entropy $S_{A}$ for $d_{A}\leq d_{B}$

$\displaystyle
S_{A}=\sum_{k=d_{B}+1}^{d_{A}d_{B}}\frac{1}{k}-\frac{d_{A}-1}{2d_{B}}.$ (42)

This agrees with the famous Page’s result Page:1993df .444 The exact $S_{A}$
in (42) is conjectured in Page:1993df  and it is later proved in Foong1994 ;
Sanchez1995 ; Sen1996 . $S_{A}$ in the opposite regime $d_{A}>d_{B}$ is
obtained from (42) by exchanging the role of $d_{A}$ and $d_{B}$. Note that
the first term of (42) is written as

$\displaystyle\sum_{k=d_{B}+1}^{d_{A}d_{B}}\frac{1}{k}=H_{d_{A}d_{B}}-H_{d_{B}},$
(43)

where $H_{m}=\sum_{k=1}^{m}1/k$ denotes the harmonic number.

Figure 1:  Plot of the entanglement entropy $S_{A}$ as a function of $\log
d_{A}$. We have set $d_{B}=20$ in this figure. The blue solid curve represents
the exact result of $S_{A}$ in (42). The orange dashed curve is the leading
approximation $S_{A}=\log d_{A}~{}(d_{A}\leq d_{B})$ and $S_{A}=\log
d_{B}~{}(d_{A}>d_{B})$.

In Fig. 1, we show the plot of $S_{A}$ as a function of $\log d_{A}$ with a
fixed $d_{B}$. As we can see from Fig. 1, $S_{A}$ grows like $\log d_{A}$ for
small $d_{A}$ and approaches $\log d_{B}$ for large $d_{A}$. At least
qualitatively, this reproduces the Page curve of the Hawking radiation
Page:1993wv  if we regard the subsystem $A$ as the radiation and the subsystem
$B$ as the black hole and $t=\log d_{A}$ as time. Around the Page time $t=\log
d_{B}$, the contributions from the fully disconnected and the fully connected
geometries exchange dominance in the replica computation of $S_{A}$
Penington:2019kki ; Almheiri:2019qdq .

Next, let us compute the capacity of entanglement $C_{A}$ by the replica
method (6). To do this, we need to compute the second derivative of
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ at $n=1$

$\displaystyle\partial_{n}^{2}\langle\operatorname{Tr}\rho_{A}^{n}\rangle\Big{|}_{n=1}=$
$\displaystyle\lim_{n\to
1}\sum_{k=1}^{d_{A}}\frac{\Gamma(d_{A})\Gamma(d_{B}+1+n-k)\Gamma(d_{A}d_{B}+1)}{\Gamma(d_{A}+1-k)\Gamma(d_{B}+1)\Gamma(d_{A}d_{B}+n)}$
(44) $\displaystyle\times$
$\displaystyle\Biggl{[}N_{n,k}\Bigl{(}\psi(d_{B}+1+n-k)-\psi(d_{A}d_{B}+n)\Bigr{)}^{2}$
$\displaystyle+2\partial_{n}N_{n,k}\Bigl{(}\psi(d_{B}+1+n-k)-\psi(d_{A}d_{B}+n)\Bigr{)}$
$\displaystyle+N_{n,k}\Bigl{(}\psi_{1}(d_{B}+1+n-k)-\psi_{1}(d_{A}d_{B}+n)\Bigr{)}+\partial_{n}^{2}N_{n,k}\Biggr{]},$

where $\psi_{1}(z)=\frac{d}{dz}\psi(z)$ denotes the trigamma function. Using
the relation

$\displaystyle\psi_{1}(m+1)=\frac{\pi^{2}}{6}-\sum_{k=1}^{m}\frac{1}{k^{2}},\quad(m\in\mathbb{N}),$
(45)

and the behavior (33) of $N_{n,k}$ near $n=1$, after some algebra we find the
exact result of capacity $C_{A}$ for $d_{A}\leq d_{B}$

$\displaystyle
C_{A}=\sum_{k=d_{B}+1}^{d_{A}d_{B}}\frac{1}{k^{2}}-\frac{(d_{A}-1)(d_{A}+3)}{4d_{B}^{2}}+\frac{d_{A}-1}{d_{B}}-\sum_{k=3}^{d_{A}}\frac{2}{k(k-2)(k-1)^{2}}\prod_{i=0}^{k-2}\frac{d_{A}-1-i}{d_{B}-i}.$
(46)

This is our main result. $C_{A}$ in the opposite regime $d_{A}>d_{B}$ is
obtained from (46) by exchanging $d_{A}$ and $d_{B}$ using the symmetry (16).
Note that the first term of (46) is written as

$\displaystyle\sum_{k=d_{B}+1}^{d_{A}d_{B}}\frac{1}{k^{2}}=H_{d_{A}d_{B}}^{(2)}-H_{d_{B}}^{(2)},$
(47)

where $H_{m}^{(2)}=\sum_{k=1}^{m}1/k^{2}$ denotes the generalized harmonic
number of the 2nd order. This is similar to the first term of $S_{A}$ in (42),
but the other terms in $C_{A}$ are more complicated than $S_{A}$.

One can easily check that (46) reduces to $C_{A,\text{planar}}$ in (36) in the
planar limit (22). Also, one can check that $C_{A}(d_{A},d_{B})$ in (46) for
$d_{A},d_{B}=2,3$ agree with the result in deBoer:2018mzv

$\displaystyle C_{A}(2,2)$ $\displaystyle=\frac{13}{36},$ (48) $\displaystyle
C_{A}(2,3)$ $\displaystyle=\frac{1169}{3600},$ $\displaystyle C_{A}(3,3)$
$\displaystyle=\frac{2898541}{6350400}.$

From the exact result (46), we find the small $d_{A}$ and the large $d_{A}$
behavior of $C_{A}$

$\displaystyle C_{A}\approx\left\\{\begin{aligned}
&\left(d_{A}-\frac{1}{d_{A}}\right)\frac{1}{d_{B}},&\quad&(1\leq d_{A}\ll
d_{B}),\\\ &\left(d_{B}-\frac{1}{d_{B}}\right)\frac{1}{d_{A}},&\quad&(d_{A}\gg
d_{B}).\end{aligned}\right.$ (49)

Our exact $C_{A}$ in (46) takes the maximal value at $d_{A}=d_{B}$

$\displaystyle
C_{A}^{(\text{max})}=\sum_{k=1}^{d_{B}^{2}}\frac{1}{k^{2}}+\sum_{k=1}^{d_{B}}\frac{1}{k^{2}}+\frac{1}{d_{B}}-\frac{1}{4d_{B}^{2}}-\frac{11}{4}.$
(50)

In the large $d_{B}$ limit this is expanded as

$\displaystyle
C_{A}^{(\text{max})}=\frac{\pi^{2}}{3}-\frac{11}{4}-\frac{3}{4d_{B}^{2}}+\mathcal{O}(d_{B}^{-3}),$
(51)

where the first two terms agree with the maximal value of capacity in the
planar limit (37). For finite $d_{A},d_{B}$, one can show that the exact
$C_{A}$ is bounded from above

$\displaystyle C_{A}\leq C_{A}^{(\text{max})}<\frac{\pi^{2}}{3}-\frac{11}{4}.$
(52)

Figure 2:  Plot of the capacity of entanglement $C_{A}$ as a function of $\log
d_{A}$. We have set $d_{B}=20$ in this figure. The blue solid curve is the
exact result of $C_{A}$ in (46). The orange dashed curves represent the
asymptotic behavior of $C_{A}$ in (49). The dashed vertical line is at
$d_{A}=d_{B}$ where $C_{A}$ becomes maximal.

In Fig. 2, we show the plot of the exact capacity $C_{A}$ in (46). We can see
that $C_{A}$ vanishes for $d_{A}=1$ and approaches zero at large $d_{A}\gg
d_{B}$. This is qualitatively similar to the result of the planar limit found
in Kawabata:2021hac . We emphasize that our result (46) is exact at finite
$d_{A},d_{B}$ and (46) includes all the non-planar corrections. As we argued
in the previous section, $C_{A}$ is sensitive to the sub-leading terms in
$\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ corresponding to the partially
connected geometries. Indeed, $C_{A}$ vanishes at the early and late “time”
$t=\log d_{A}$ where the fully connected or fully disconnected geometry is
dominant. $C_{A}$ takes a non-zero value near the Page time $t=\log d_{B}$ (or
$d_{A}=d_{B}$) which is interpreted that the partially connected geometries
give substantial contributions to $C_{A}$ near the Page time.

##  5 Conclusions and outlook

In this paper, we have computed the exact capacity of entanglement $C_{A}$
(46) at finite $d_{A},d_{B}$ using the replica method (6). At the technical
level, the important ingredient in our computation is the new exact formula
(21) of $\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ written in terms of the
Narayana number $N_{n,k}$. This formula (21) makes the relation to the planar
limit manifest. We argued that $C_{A}$ vanishes for the fully connected or
fully disconnected geometries, and $C_{A}$ is sensitive to the sub-leading
contributions to $\langle\operatorname{Tr}\rho_{A}^{n}\rangle$ coming from the
partial connected geometries in the gravitational path integral. This suggests
that $C_{A}$ is a good probe of the partial entanglement, as discussed in
deBoer:2018mzv .

There are several open questions. The capacity of entanglement is introduced
in qi  as an analogue of the heat capacity. Indeed, if we introduce the
modular Hamiltonian $K=-\log\rho_{A}$, the moment
$\operatorname{Tr}\rho_{A}^{n}$ looks like the partition function

$\displaystyle Z_{n}=\operatorname{Tr}\rho_{A}^{n}=\operatorname{Tr}e^{-nK},$
(53)

and $n$ plays the role of the inverse temperature $\beta$. In this picture,
our definition (6) of $S_{A}$ and $C_{A}$ is based on the “annealed” free
energy $\log\langle Z_{n}\rangle$

$\displaystyle S_{A}=-\partial_{n}\log\langle Z_{n}\rangle\Big{|}_{n=1},\qquad
C_{A}=\partial_{n}^{2}\log\langle Z_{n}\rangle\Big{|}_{n=1}.$ (54)

One could consider the quenched version of $S_{A}$ and $C_{A}$ as well

$\displaystyle S_{A}^{\text{qu}}=-\partial_{n}\langle\log
Z_{n}\rangle\Big{|}_{n=1},\qquad C_{A}^{\text{qu}}=\partial_{n}^{2}\langle\log
Z_{n}\rangle\Big{|}_{n=1}.$ (55)

We leave the computation of the quenched version of $S_{A}$ and $C_{A}$ as an
interesting future problem.

It would be interesting to study the gravitational picture of the capacity of
entanglement. $C_{A}$ is related to the quantum fluctuation of the modular
Hamiltonian and the prescription of the gravitational computation of $C_{A}$
is proposed in Nakaguchi:2016zqi . We have argued that $C_{A}$ receives
contributions only from the sub-leading partially connected geometries in the
replica computation. It would be interesting to related this picture to the
prescription in Nakaguchi:2016zqi .

###### Acknowledgements.

This work was supported in part by JSPS KAKENHI Grant No. 19K03845.

## References

  * (1) G. Penington, S. H. Shenker, D. Stanford, and Z. Yang, “Replica wormholes and the black hole interior,” arXiv:1911.11977 [hep-th].
  * (2) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian, and A. Tajdini, “Replica Wormholes and the Entropy of Hawking Radiation,” JHEP 05 (2020) 013, arXiv:1911.12333 [hep-th].
  * (3) D. N. Page, “Information in black hole radiation,” Phys. Rev. Lett. 71 (1993) 3743–3746, arXiv:hep-th/9306083.
  * (4) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian, and A. Tajdini, “The entropy of Hawking radiation,” arXiv:2006.06872 [hep-th].
  * (5) D. N. Page, “Average entropy of a subsystem,” Phys. Rev. Lett. 71 (1993) 1291–1294, arXiv:gr-qc/9305007.
  * (6) K. Kawabata, T. Nishioka, Y. Okuyama, and K. Watanabe, “Probing Hawking radiation through capacity of entanglement,” arXiv:2102.02425 [hep-th].
  * (7) H. Yao and X.-L. Qi, “Entanglement entropy and entanglement spectrum of the Kitaev model,” Phys. Rev. Lett. 105 (2010) 080501, arXiv:1001.1165 [cond-mat].
  * (8) Y. Nakaguchi and T. Nishioka, “A holographic proof of Rényi entropic inequalities,” JHEP 12 (2016) 129, arXiv:1606.08443 [hep-th].
  * (9) Y. O. Nakagawa and S. Furukawa, “Capacity of entanglement and the distribution of density matrix eigenvalues in gapless systems,” Phys. Rev. B 96 no. 20, (2017) 205108, arXiv:1708.08924 [cond-mat.str-el].
  * (10) J. De Boer, J. Järvelä, and E. Keski-Vakkuri, “Aspects of capacity of entanglement,” Phys. Rev. D 99 no. 6, (2019) 066012, arXiv:1807.07357 [hep-th].
  * (11) B. Collins and I. Nechita, “Random matrix techniques in quantum information theory,” J. Math. Phys. 57 (2016) 015215, arXiv:1509.04689 [quant-ph].
  * (12) I. Nechita, “Asymptotics of random density matrices,” Ann. Henri Poincaré 8 (2007) 1521, arXiv:quant-ph/0702154.
  * (13) P. J. Hanlon, R. P. Stanley, and J. R. Stembridge, “Some combinatorial aspects of the spectra of normally distributed random matrices,” Contemp. Math 138 (1992) 151–174.
  * (14) U. Haagerup and S. Thorbjørnsen, “Random matrices with complex Gaussian entries,” Expo. Math. 21 (2003) 293.
  * (15) E. Lubkin, “Entropy of an $n$‐system from its correlation with a $k$‐reservoir,” J. Math. Phys. 19 (1978) 1028.
  * (16) H.-J. Sommers and K. Zyczkowski, “Statistical properties of random density matrices,” J. Phys. A: Math. Gen. 34 (2004) 8457, arXiv:quant-ph/0405031.
  * (17) J. Kudler-Flam, “Relative Entropy of Random States and Black Holes,” arXiv:2102.05053 [hep-th].
  * (18) J. P. Dyer, “Divergence of Lubkin’s series for a quantum subsystem’s mean entropy,” arXiv:1406.5776 [cond-mat].
  * (19) S. K. Foong and S. Kanno, “Proof of Page’s conjecture on the average entropy of a subsystem,” Phys. Rev. Lett. 72 (1993) 1148.
  * (20) J. Sánchez-Ruiz, “Simple proof of Page’s conjecture on the average entropy of a subsystem,” Phys. Rev. E 52 (1995) 5653.
  * (21) S. Sen, “Average Entropy of a Quantum Subsystem,” Phys. Rev. Lett. 77 (1996) 1.