# A lower bound for the spectral gap of the conjugate Kac process with 3
interacting particles

Luís Simão Ferreira Supported by CMAF-CIO grant UIDP/04561/2020.

###### Abstract

In this paper, we proceed as suggested in the final section of [2] and prove a
lower bound for the spectral gap of the conjugate Kac process with 3
interacting particles. This bound turns out to be around $0.02$, which is
already physically meaningful, and we perform Monte Carlo simulations to
provide a better empirical estimate for this value via entropy production
inequalities. This finishes a complete quantitative estimate of the spectral
gap of the Kac process.

###### Acknowledgements

The author would like to thank CMAF-CIO for the opportunity to develop this
work and professor José Francisco Rodrigues for bringing this topic to my
attention.

I also would like to extend my deepest gratitude to professor Eric Carlen, my
supervisor, for all the patience, clear and intuitive explanations and
everything that I have learned regarding functional analysis and mathematical
physics for the duration of this project.

###### Contents

  1. 1 Introduction
    1. 1.1 The Kac process
    2. 1.2 Structure of the paper
  2. 2 The conjugate Kac process
    1. 2.1 Spectral gap bounds via entropy production inequalities
    2. 2.2 Monte Carlo simulations
  3. 3 Proof of the main result
    1. 3.1 The correlation operator on the sphere
    2. 3.2 Anti-symmetric sector
    3. 3.3 Symmetric sector
      1. 3.3.1 Large $\ell$ case
      2. 3.3.2 Small $\ell$
  4. 4 Conclusion
  5. A Small $\ell$ computations
  6. B Python code for the Monte Carlo simulations

##  1 Introduction

$N$-particle systems have been a target of flourishing research for the better
part of the last 50 years, in particular due to their connection to classical
statistical physics PDE’s via hydrodynamic limits. Thoroughly inspecting these
collision processes allows us to better understand limiting cases such as the
Boltzmann equation, while avoiding dealing with the equation explicitly, since
it involves strong non-linearity.

In 1956, Marc Kac [1] introduced what we refer to as the Kac process, a
reversible Markov process for $N$ spatially homogeneous particles with one
dimensional velocities, which describes the evolution of the system as
particles randomly collide in a momentum and energy conserving fashion.

A behaviour exhibited in many naturally occurring particle systems, and a
desirable property of these models, is that of exponential decay to
equilibrium. The easiest way one can prove such property is by showing the
existence of a spectral gap for the generator of the process, and a uniform
bound lets us draw conclusions as the number of particles becomes infinite.
Recall that the spectral gap of an operator is defined as the infimum of its
positive spectrum, which may coincide with the lowest non-zero eigenvalue, if
these exist.

In essence, under the Kac process, if $L$ is the infinitesimal generator and
$f_{0}$ is an initial distribution in our state space, then it evolves as
$f_{t}=e^{tL}f_{0}$, and the existence of a spectral gap $\Delta$ implies an
exponential decay with rate at least equal to $\Delta$.

It is hard to relate relaxation in the $N$-particle system to relaxation in
the limit equation, unless one has uniform bounds on the spectral gap,
independent of $N$. Kac conjectured this to be true in his toy model, and this
was successfuly proven later by Janvresse [5], but with one dimensional
velocities and uniform jump rates, which nevertheless allows us to efficiently
generate uniform points on the $n$-sphere.

In a series of papers, Carlen, Carvalho and Loss [2] managed to prove the
existence of a spectral gap for a one parameter family of Kac processes in
3-dimensions, undergoing more general collisions with rates that are not
bounded away from zero. We will explore next the background of these papers.

###  1.1 The Kac process

More precisely, consider an additional physical parameter $0\leq\alpha\leq 2$
and an $N$-tuple of velocities $v_{j}$ in ${{\mathbb{R}}}^{3}$ such that

$\frac{1}{N}\sum_{j=1}^{N}v_{j}=0,\frac{1}{N}\sum_{j=1}^{N}|v_{j}|^{2}=1,$ (1)

meaning that the total momentum is zero and the average kinetic energy per
particle (with mass 2) is 1. These constraints describe our state space
$\mathcal{S}_{N}$, isometric to a sphere with $3N-4$ dimensions and radius
$\sqrt{N}$.

Before we move on, it is important to fix some notation. We will refer to
elements of this state space $\mathcal{S}_{N}$ as a single letter $v$, $y$,
and sometimes as a vector $\vec{v}$, $\vec{y}$, if it’s not clear from
context. Letters with subscripts such as $v_{j}$ shall always denote
velocities in ${{\mathbb{R}}}^{3}$.

When the process begins, there is an associated exponential clock $T_{i,j}$
for each pair of velocities $(v_{i},v_{j})$, with parameter

$\lambda_{i,j}=N\binom{N}{2}^{-1}|v_{i}-v_{j}|^{\alpha},$ (2)

and the minimum

$T:=\min_{i,j}T_{i,j}$ (3)

yields our first collision time (and colliding pair). Given the restrictions
above, these collisions can be parametrized in terms of a vector
$\sigma\in\mathbb{S}^{2}$:

$\displaystyle v_{i}^{*}(\sigma)$
$\displaystyle=\frac{v_{i}+v_{j}}{2}+\frac{|v_{i}-v_{j}|}{2}\sigma,$ (4)
$\displaystyle v_{j}^{*}(\sigma)$
$\displaystyle=\frac{v_{i}+v_{j}}{2}-\frac{|v_{i}-v_{j}|}{2}\sigma.$ (5)

Another source of randomness on the system is then the determination of this
direction vector $\sigma$. The particular choice we follow is uniform
sampling, as it turns out to be more physically relevant, even though other
distributions or rules could be used.

The cases $0\leq\alpha\leq 1$, with 0 being the Maxwellian molecule case and 1
the hard sphere case, are the ones of physical interest, but it turns out, as
shown by Villani (regarding entropy production inequalities) [6] and as seen
in the work of Carlen, Carvalho and Loss, that the ”super hard sphere”
$\alpha=2$ case provides powerful insights for lower values of the parameter,
and will later be our main focus.

Their proof is based on chaos estimates that quantify the ”near independence”
of a finite set of coordinates. If $\sigma_{N}$ denotes the uniform measure on
the sphere $\mathbb{S}^{N}(\sqrt{N})$ with radius $\sqrt{N}$, it is known that

$\lim_{N\to+\infty}\int_{\mathbb{S}^{N}(\sqrt{N})}\phi(v_{1},...,v_{k})d\sigma_{N}=\int_{{{\mathbb{R}}}^{k}}\phi(v_{1},...,v_{k})d\gamma^{\otimes
k},$ (6)

with $\gamma$ being the Gaussian measure in ${{\mathbb{R}}}$, and we say that
the measure $\sigma_{N}$ is $\gamma$-chaotic. In general, if the above
equality holds for $\mu$ and a sequence of measures $\mu_{N}$, $\mu_{N}$ is
said to be $\mu$-chaotic.

Another important step of the proof is the introduction of a conjugate
process, which allows for a clever induction argument. With the right choice
of coordinates, integrals over the state space can be factored as

$\int_{\mathcal{S}_{N}}\phi(\vec{v})d\sigma_{N}=\int_{B}\bigg{[}\int_{S_{N-1}}\phi(T_{j}(\vec{y},v))d\sigma_{N-1}\bigg{]}d\nu_{N}(v),$
(7)

with $B$ the unit ball in ${{\mathbb{R}}}^{3}$. $T_{1}:\mathcal{S}_{N-1}\times
B\to\mathcal{S}_{N}$ is given by

$T_{1}(y,v)=\bigg{(}\sqrt{N-1}v,\beta(v)y_{1}-\frac{1}{\sqrt{N-1}}v,...,\beta(v)y_{N-1}-\frac{1}{\sqrt{N-1}}v\bigg{)},$
(8)

while $T_{j}$ is defined analogously but with $v$ going to the $j$-th
coordinate, and

$\beta^{2}(v):=\frac{N}{N-1}(1-|v|^{2}).$ (9)

Finally, the measure $\nu_{N}$ on $B$ is defined as

$d\nu_{N}(v)=\frac{|\mathbb{S}^{3N-7}|}{|\mathbb{S}^{3N-4}|}(1-|v|^{2})^{(3N-8)/2}dv,$
(10)

which is what we refer to as the equilibrium radial measure, as it is the
radial marginal of the velocities for the uniform measure in
$\mathcal{S}_{N}$.

As we will verify later, this formula allows for an inductive inequality of
the type

$\Delta_{N,\alpha}\geq\frac{N}{N-1}\Delta_{N-1,\alpha}\hat{\Delta}_{N,\alpha},$
(11)

where $\hat{\Delta}_{N,\alpha}$ is the spectral gap of a conjugate Kac
process, with vastly simplified dynamics. Furthermore, the authors make use of
this fact to prove that there are $K_{0}>0$ and $N_{0}\in{{\mathbb{N}}}$ such
that

$\Delta_{N,\alpha}\geq K_{0}\Delta_{N_{0},\alpha}>0,$ (12)

leaving only the base case to show. For $N=2$, the process is trivial, as the
constraints completely determine one velocity from the other. Finally, making
use of a compactness argument to show that $\hat{\Delta}_{3,\alpha}>0$, it
follows that

$\Delta_{3,\alpha}\geq\frac{3}{2}\Delta_{2,\alpha}\hat{\Delta}_{3,\alpha}>0.$
(13)

However, unlike previous bounds that could be explicitly computed, this does
not give any quantitative information about $\hat{\Delta}_{3,\alpha}$, besides
positivity, and it remains to be shown that this gap is physically meaningful
when interpreted as an expected relaxation time. In this paper, we compute a
lower bound for the spectral gap of the conjugate Kac process with 3
particles.

Even though the proof provided in [2] makes use of asymptotic methods in the
form of quantitative chaos estimates, the results are strong enough to hold
already for all $N\geq 4$, and provide a positive gap. This fails for $N=3$.

Due to the fact that $N=2$ is already trivial, and most methods allow for good
estimates for $N>3$, one would expect this case to be accessible.
Unfortunately, this is not true, as most methods fail to provide a positive
gap in this number of dimensions.

###  1.2 Structure of the paper

This paper is organized as follows.

In section 2, we introduce our main object of study, the conjugate Kac
process, building up the necessary tools and finally stating our main result.
We also explore the possibility of obtaining spectral gap bounds via entropy
production inequalities and, while not being able to prove it analytically,
present some numerical evidence of such results.

The proof is developed over section 3. We begin by inspecting the spectrum of
our main tool, the correlation operator $K$, and proving suitable upper bounds
on the eigenvalue sequence $\kappa_{n,\ell}$.

The proof itself is then broken down over the anti-symmetric and symmetric
sectors, as they are orthogonal to each other. The first one is handled fairly
quickly, while the symmetric sector, which is of bigger physical interest due
to its relation to the limit equation, requires a more hands-on approach. We
treat different angular momentum sectors, where different methods handle
’large’ and ’small’ values of $\ell$, leaving only a finite amount of
eigenvalues to be computed.

##  2 The conjugate Kac process

We will now define the family of reversible Markov jump processes on
$\mathcal{S}_{N}$, conjugate to the Kac process. Fixing a number of
interacting particles $N$, $0<\alpha\leq 2$, and given $\vec{v}\in S_{N}$, let
$\\{\hat{T}_{1},...,\hat{T}_{N}\\}$ be a set of $N$ independent exponential
variables with parameter $\lambda_{k}(\vec{v})$ of $\hat{T}_{k}$ given by

$\lambda_{k}(\vec{v})=\frac{1}{N}\bigg{[}\frac{N^{2}-(1+|v_{k}|^{2})N}{(N-1)^{2}}\bigg{]}^{\alpha/2}.$
(14)

Due to the restraints on energy and momentum, the maximum of $|v_{k}|^{2}$ is
$N-1$, and thus the parameter is always positive.

Let $\hat{T}=\min\\{\hat{T}_{1},...,\hat{T}_{N}\\}$ be the first jump time,
and let $k$ be the corresponding index. At this time, the process makes a
conditional jump to uniform to a new point $\vec{v}\in S_{N}$, fixing $v_{k}$,
and all other velocities are re-sampled uniformly, conditional on $v_{k}$.
Even though this process is trivial for $N=2$, since the momentum equation
completely determines one velocity from the other, our $N=3$ case is subtle.

The generator of this process can be obtained via projection (or conditional
expectation) operators. Recall that, for any $\phi\in L^{2}(\sigma_{N})$,
there is a unique element $f(v_{k})$ of $L^{2}(\sigma_{N})$ such that

$\int_{S_{N}}\phi(\vec{v})g(v_{k})d\sigma_{N}=\int_{S_{N}}f(v_{k})g(v_{k})d\sigma_{N},$
(15)

and we define this way the conditional expectation of $\phi$, given $v_{k}$:

$P_{k}\phi:={{\mathbb{E}}}[\phi|v_{k}].$ (16)

This can also be understood as an orthogonal projection onto the subspace of
square integrable functions that depend only on the $k$-th coordinate.
Furthermore, the factorization formula (7) gives us an expression for this
operator:

$P_{k}\phi(v)=\int_{\mathcal{S}_{N-1}}\phi(T_{k}(y,v_{k}/\sqrt{N-1})d\sigma_{N-1}.$
(17)

The generator of this process can then be written as

$\hat{L}_{N,\alpha}f=-\frac{1}{N}\sum_{k=1}^{N}\frac{1}{N}\bigg{[}\frac{N^{2}-(1+|v_{k}|^{2})N}{(N-1)^{2}}\bigg{]}^{\alpha/2}[f-P_{k}f],$
(18)

and it is also useful to define

$w_{N}(v_{k}):=\frac{N^{2}-(1+|v_{k}|^{2})N}{(N-1)^{2}}.$ (19)

It is easy to prove self-adjointness, as we would expect from the
reversibility of the process, and its null space is spanned by the constants.
Similarly as before, the spectral gap of the conjugate Kac process is defined
as

$\hat{\Delta}_{N,\alpha}:=\inf\\{\mathcal{D}_{N,\alpha}(f,f)|\langle
f,1\rangle_{L^{2}(\sigma_{N})},\|f\|^{2}_{L^{2}(\sigma_{N})}=1\\},$ (20)

with

$\displaystyle\mathcal{D}_{N,\alpha}(f,f)$ $\displaystyle:=-\langle
f,\hat{L}_{N,\alpha}f\rangle_{L^{2}(\sigma_{N})}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\frac{1}{N}\int_{S_{N}}\bigg{[}\frac{N^{2}-(1+|v_{k}|^{2})N}{(N-1)^{2}}\bigg{]}^{\alpha/2}[f^{2}-fP_{k}f]d\sigma_{N}.$
(21)

We are at this moment in conditions of explaining how this process arises. Let
$L_{N,\alpha}$ be the infinitesimal generator of the Kac process,

$\mathcal{E}_{N,\alpha}(f,f)=-\langle f,L_{N,\alpha}f\rangle$ (22)

the Dirichlet form associated with $L_{N,\alpha}$ and
$\mathcal{E}_{N,\alpha}(f,f|v_{k})$ the conditional Dirichlet form obtained by
factorizing, so that

$\mathcal{E}_{N,\alpha}(f,f)=\frac{N}{N-1}\bigg{(}\frac{1}{N}\sum_{k=1}^{N}\int_{B}w_{N}^{\alpha/2}(v_{k})\mathcal{E}_{N,\alpha}(f,f|v_{k})d\nu_{N}(v_{k}/\sqrt{N-1})\bigg{)}.$
(23)

Then, using the spectral gap for $N-1$ particles, one has

$\mathcal{E}_{N,\alpha}(f,f)\geq\frac{N}{N-1}\Delta_{N-1,\alpha}\bigg{(}\frac{1}{N}\sum_{k=1}^{N}\int_{\mathcal{S}_{N}}w_{N}^{\alpha/2}(v_{k})[f-P_{k}f]^{2}d\sigma_{N}\bigg{)},$
(24)

and the quantity in brackets corresponds exactly to (2).

Now, define the self-adjoint operators

$W^{(\alpha)}=W^{(\alpha)}(v):=\frac{1}{N}\sum_{k=1}^{N}w_{N}^{\alpha/2}(v_{k})$
(25)

and

$P^{(\alpha)}=\frac{1}{N}\sum_{k=1}^{N}w_{N}^{\alpha/2}(v_{k})P_{k},$ (26)

an average of weighted projections, so that we have
$\hat{L}_{N,\alpha}=W^{(\alpha)}-P^{(\alpha)}$. For $\alpha=2$ and $N=3$, we
see that $W^{(2)}=\frac{3}{4}$, and (20) becomes

$\hat{\Delta}_{3,2}=\frac{3}{4}-\mu_{3},$ (27)

with

$\mu_{3}:=\sup\\{\langle
f,P^{(2)}f\rangle_{L^{2}(d\sigma_{3})}|\|f\|_{L^{2}(d\sigma_{3})}=1,\langle
f,1\rangle_{L^{2}(d\sigma_{3})}=0\\}.$ (28)

The spectrum of $P^{(\alpha)}$ is also studied on [3]. More importantly, for
$\alpha=2$, the following is true

###### Lemma 2.1

The essential spectrum of $P^{(2)}$ is the interval $[0,\frac{1}{2}]$, and the
remaining spectrum can only consist of isolated eigenvalues in
$\big{(}\frac{1}{2},\frac{3}{4}\big{)}$.

Note that, in particular, $\frac{3}{4}$ cannot be an accumulation point, and
thus this proves that the spectral gap exists. Moreover, if
$\mu_{3}\leq\frac{1}{2}$, we immediately have that $\hat{\Delta}_{3,2}\geq
0.25$, and if $\frac{1}{2}<\mu_{3}<\frac{3}{4}$, it is an eigenvalue of
$P^{(2)}$.

We then turn to the eigenvalue equation $P^{(2)}f=\lambda f$ to obtain the
only missing estimate for a fully quantitative result on the three particle
spectral gap $\hat{\Delta}_{3,2}$. In [2], the authors break this down to the
symmetric and anti-symmetric eigen-sectors, as they are mutually orthogonal. A
computation for the anti-symmetric case is provided, and a suggestion is made
on how to proceed for the remaining case. In section 3.2 we explain and
slightly improve the former, while the latter is treated in section 3.3 \-
this last case proved to be much harder, due to the analytical difficulties in
dealing with the spectrum of the conditional expectation operator (see section
3.1).

Using this method, we were able to prove the following main result:

###### Theorem 2.1

The spectral gap of the conjugate Kac process in 3 dimensions, with
$\alpha=2$, is at least $0.01984$.

The actual computations we provide show a result slightly less then $0.02$ in
order to simplify and reduce the amount of cases, but could easily be adapted
to achieve this quantity at the expense of explicitly determining a larger
number of eigenvalues. Note that it is expected that the true value of the gap
is larger, and we will illustrate this next.

###  2.1 Spectral gap bounds via entropy production inequalities

Cedric Villani [6] was the first to prove entropy production inequalities for
variants of the Kac process, namely with one dimensional velocities. Even
though obtaining such bounds for the 3-dimensional Kac process is still an
open problem, one can adopt a similar reasoning and attempt to derive entropy
production inequalities for the simpler conjugate process. Recall that the
entropy functional of a suitable density $f$ on $\mathcal{S}_{N}$ is defined
as

$H(f):=\int_{\mathcal{S}_{N}}f\log fd\sigma_{N},$ (29)

and the entropy production of a process with generator $\hat{L}_{N,\alpha}$ as

$D_{\hat{L}_{N,\alpha}}(f):=-\frac{d}{dt}\bigg{(}H(e^{t\hat{L}_{N,\alpha}}f)\bigg{)}\geq
0,$ (30)

so that we can look for estimates of the type

$D_{\hat{L}_{N,\alpha}}(f)\geq C_{N,\alpha}H(f),$ (31)

with $C_{N,\alpha}$ a constant independent of $f$, and depending on
$\hat{L}_{N,\alpha}$ only via the parameters describing the process.

Results of the above type are stronger than spectral gap estimates, as a
classical linearization procedure implies the latter. More precisely, taking
$f:=1+\varepsilon h$, with $\varepsilon>0$ and

$\int_{\mathcal{S}_{N}}hd\sigma_{N}=0,$

and expanding in first order leads us to

$\hat{\Delta}_{N,\alpha}\geq\frac{1}{2}C_{N,\alpha}.$ (32)

One advantage of the entropy production approach is that we can actually
measure the empirical entropy in Monte Carlo simulations, and obtain an
estimate of the true spectral gap by analyzing the log-entropy decay, which we
will do next.

First, by symmetrizing one can compute the entropy production functional of
the conjugate process:

$D_{N,\alpha}(f):=\frac{1}{N}\sum_{k=1}^{N}\int_{\mathcal{S}_{N}}w_{N}^{\alpha/2}(v_{k})\log
f[f-P_{k}f]d\sigma_{N}.$ (33)

Specializing to $\alpha=2$, and using the convexity of $-\log t$ and the fact
that

$\frac{1}{N}\sum_{k=1}^{N}w_{N}(v)=1,$ (34)

we obtain

$\displaystyle D_{N,2}(f)\geq\bar{D}_{N,2}(f):=$
$\displaystyle\frac{1}{N}\sum_{k=1}^{N}\int_{\mathcal{S}_{N}}w_{N}(v_{k})[f\log
f-P_{k}f\log P_{k}f]d\sigma_{N}$ (35) $\displaystyle=$ $\displaystyle
H(f)-\frac{1}{N}\sum_{k=1}^{N}\int_{\mathcal{S}_{N}}w_{N}(v_{k})[P_{k}f\log
P_{k}f-P_{k}f+1]d\sigma_{N}.$ (36)

Now, as $t\log t-t+1\geq 0$ for all $t>0$ and $w_{N}(v)\leq\frac{N}{N-1}$ for
all $v$, we see that the $k$-th summand in the second term is less than
$\frac{N}{N-1}H(P_{k}f)$, and thus we have

$\bar{D}_{N,2}(f)\geq H(f)-\frac{1}{N-1}\sum_{k=1}^{N}H(P_{k}f)$ (37)

To continue, we need to somehow relate the entropy of the marginals $P_{k}f$
with the entropy of the joint distribution $f$. Fortunately, this is precisely
developed in [7] and [8], and we make use of the CCL inequality (in 3
dimensions),

$\frac{N-1}{N}\sum_{k=1}^{N}H(P_{k}f)\leq 2H(f),$ (38)

thus giving us

$\bar{D}_{N,2}(f)\geq\bigg{(}1-\frac{2N}{(N-1)^{2}}\bigg{)}H(f).$ (39)

As for the much simpler case of $\alpha=0$ with constant jump rates, one can
easily verify by slightly adapting the computations above that this yields

$C_{N,0}=\frac{N-3}{N-1}.$ (40)

With this quick argument, we can then state the following:

###### Proposition 2.2

The entropy production functional of the Kac process in 3 dimensions,
$D_{N,\alpha}$, satisfies

$D_{N,2}(f)\geq\bigg{(}1-\frac{2N}{(N-1)^{2}}\bigg{)}H(f),$ (41)

and

$D_{N,0}(f)\geq\frac{N-3}{N-1}H(f).$ (42)

Thus, for $N\geq 4$, the spectral gap of the conjugate Kac process in 3
dimensions is at least $\frac{1}{18}$ for the super hard sphere case, and
$\frac{1}{6}$ for the Maxwellian molecule case.

Sadly, this is not quite enough, and the bounds degenerate at $N=3$, and
adapting Villani’s approach seems to give degenerate bounds at $N=3$ as well,
even though the technique is much stronger than the one we use above. The main
culprit seems to be the CCL inequality (38) as the extra $N/(N-1)$ factor,
which isn’t present for the 1-dimensional case, loses us enough for the method
to fail.

This makes it clear that we’re dealing with an edge case, just outside the
reach of most methods, even if we can simultaneously handle both the entropy
production and the spectral gap problems for all $N$ above 3.

Thus, a careful analysis of the dynamics and geometry of the problem is needed
to obtain such a bound, and we make use of the information we have about the
spectrum of the conditional expectation operator $K$ on the sphere and
properties of orthogonal polynomials for this purpose.

###  2.2 Monte Carlo simulations

In order to observe this relaxation to equilibrium, we perform Monte Carlo
simulations of the Kac process with 3 interacting particles. Even if we can’t
prove entropy production inequalities for the conjugate Kac process when
$N=3$, this method allows us to obtain an estimate of 0.3 for the exponential
entropy decay rate, and thus an estimate of about 0.15 for the spectral gap,
which is only about an order of magnitude larger than the theoretical result
we achieve.

Based on the factorization formula (7), we notice that the uniform measure in
$S_{3}$, $\sigma_{3}$ factorizes as

$d\sigma_{3}=d\sigma_{2}\times d\nu_{3},$

with

$d\nu_{3}:=\frac{|\mathbb{S}^{2}|}{|\mathbb{S}^{5}|}(1-|v|^{2})^{\frac{1}{2}}dv,$

induced by the transformation $T_{1}:S_{2}\times B\to S_{3}$,

$T_{1}(y,v):=\bigg{(}\sqrt{2}v,\beta(v)y_{1}-\frac{v}{\sqrt{2}},\beta(v)y_{2}-\frac{v}{\sqrt{2}}\bigg{)},$
$\beta(v):=\sqrt{\frac{3}{2}(1-|v|^{2})}.$

As the system decays to equilibrium and the distribution of the velocities
approaches uniform, the marginal distributions become closer do $d\nu_{3}$. In
particular, by symmetry, we will only look at the radial component of this
measure, which is proportional to $r^{2}(1-r^{2})^{\frac{1}{2}}dr$.
Additionally, the factorization formula and its symmetries reduce the
conditional sampling problem to sampling from $S_{2}$, which is trivial. We
refer to the radial distribution of the first particle as the sampled
distribution, and the others as the implied distributions.

Thus, we can proceed as follows. Fixing the first coordinate, we sample
$r=|v|$ from a given distribution using an acceptance/rejection mechanism, and
sample $y_{1}$, $y_{2}$ uniformly, so that the point is in $S_{3}$. We then
let the system evolve sampling first the jump time and the particle that stays
fixed, identified with $\sqrt{2}v$, and sample $y_{1}$, $y_{2}$ uniformly on
$S_{2}$.

We perform 1 million simulations of the process in Python, for $\alpha=0,2$,
and use these simulations to build a histogram approximating the probability
distribution over 100 bins and several timestamps. We then measure the
relative entropy between the empirical measure and the equilibrium radial
measure over time, and use a log-plot to showcase and estimate the exponential
decay. This was done over a selection of initial distributions, and a ”worst
case scenario” is showcased here with the sampling distribution being
$2(1-x)$, as most of the probability mass is concentrated close to zero.

In figures 1, 2 and 3, the rather quick relaxation is evident. By computing
the relative entropy using a first order Euler method, we see that there is a
log-linear relationship up to a certain point, where the distributions are
close enough so that the numerical approximations and random noise prevent the
entropy from decaying below $e^{-6}$. Meanwhile, figures 4, 5 and 6 show the
evolution of the distributions over several timestamps.

Figure 1: $\alpha=2$, sampled entropy

Figure 2: $\alpha=2$, implied entropy

Figure 3: $\alpha=0$, sampled entropy

t=0s

t=0.5s

t=2s

t=3.5s

t=5s

t=10s

Figure 4: $\alpha=0$, sampled velocity

t=0s

t=2s

t=6s

t=10s

t=16s

t=24s

Figure 5: $\alpha=2$, sampled velocity

t=0s

t=2s

t=6s

t=10s

t=16s

t=24s

Figure 6: $\alpha=2$, implied velocity

##  3 Proof of the main result

###  3.1 The correlation operator on the sphere

The main result is based on careful analysis of the eigenfunctions and the
spectrum of $K$. We turn our attention to the study of this operator.

Let $B$ be the unit ball in ${{\mathbb{R}}}^{3}$, $N\geq 3$, $\nu_{N}$ as
defined in (10) and consider the projections $\pi_{k}(\vec{v}):S_{N}\to B$,

$\pi_{k}(\vec{v}):=\frac{1}{\sqrt{N-1}}v_{k}.$ (43)

By formula (7), we then have that

$\int_{B}\psi(v)d\nu_{N}=\int_{S_{N}}\psi(\pi_{k}(\vec{v}))d\sigma_{N},$

for any integrable function $\psi$ on $B$. We then define the operator $K$ on
$L^{2}(B,\nu_{N})$ by

$\langle\psi_{1},K\psi_{2}\rangle_{L^{2}(B,\nu_{N})}:=\int_{S_{N}}\psi_{1}^{*}(\pi_{1}(\vec{v}))\psi_{2}(\pi_{2}(\vec{v}))d\sigma_{N}.$
(44)

Note that the symmetries of the variables make $K$ self-adjoint. Furthermore,
a simple rescaling, to extend $K$ to $L^{2}(\sigma_{N})$, allows us to
interpret this operator simply as a conditional expectation on the uniform
sphere

$K\xi(v)={{\mathbb{E}}}[\xi(v_{1})|v_{2}=v],$

but we will not consider this scaling as the former simplifies the treatment
of eigenvalues.

Carlen, Carvalho and Loss make use of this operator to prove quantitative
chaos estimates; that is, to quantify how far from independent the coordinates
are, asymptotically in $N$, and most of their work is based on extensive study
of the spectrum of $K$, built up over [3], [4].

Before we continue, define the functions

$\eta_{0}=1,\eta_{j}(v)=v_{j},\text{ for }1\leq j\leq 3,\text{ and
}\eta_{4}(v)=(|v|^{2}-1)/(N-1).$ (45)

We will now present the most important facts obtained in these papers for
general $N$, and later turn to the our $N=3$ case. Following the results in
[3], we have this next lemma:

###### Lemma 3.1

Let $N\geq 3$. The operator $K$ is compact. The function $\eta_{0}$ is an
eigenfunction of $K$ with eigenvalue 1, and spans the corresponding
eigenspace. The functions $\eta_{j}$, for $1\leq j\leq 4$, are eigenfunctions
with eigenvalue $-\frac{1}{N-1}$, and they form an orthogonal basis for this
subspace. No other eigenvalues of $K$ are larger in absolute value then
$\frac{5N-3}{3(N-1)^{3}}$. Thus, for all $\psi\in L^{2}(B,\nu_{N})$ orthogonal
to these $\eta_{j}$,

$\|K\psi\|_{2}^{2}\leq\frac{5N-3}{3(N-1)^{3}}\|\psi\|^{2}_{2}.$ (46)

Finally, every eigenvalue $\kappa$ of $K$, other than 1, $-\frac{1}{N-1}$ and
$\frac{5N-3}{3(N-1)^{3}}$, satisfies

$-\frac{7N-3}{3(N-1)^{4}}\leq\kappa\leq\frac{5N-3}{3(N-1)^{3}}.$ (47)

Unfortunately, these bounds are not enough when $N=3$ and are just slightly
off to allow us to obtain the main result more easily. On the other hand, the
authors also manage to establish that $K$ has a complete basis of
eigenfunctions of the form

$g_{n,\ell,m}=h_{n,\ell}(|v|)|v|^{\ell}\mathcal{Y}_{\ell,m}(v/|v|),$ (48)

where $\mathcal{Y}_{\ell,m}$ denotes the $\ell$-th degree spherical harmonic
and, with $t=|v|^{2}-1$, $\alpha=\frac{3N-8}{2}\bigg{(}=\frac{1}{2}\bigg{)}$
and $\beta=\ell+\frac{1}{2}$,

$h_{n,\ell}(|v|)=P_{n}^{(\alpha,\beta)}(t),$

the $n$-th degree Jacobi polynomial (recall that the Jacobi polynomials are
the orthogonal polynomials on the interval $[-1,1]$ for the weight
$(1-x)^{\alpha}(1+x)^{\beta}$). For convenience, we let capital $P$ denote
these polynomials under Rodrigues’ normalization,

$P_{n}^{(\alpha,\beta)}(t)=\frac{(-1)^{n}}{2^{n}n}(1-x)^{-\alpha}(1+x)^{-\beta}\frac{d^{n}}{dx^{n}}\bigg{(}(1-x)^{\alpha+n}(1+x)^{\beta+n}\bigg{)},$
(49)

while $p_{n}^{(\alpha,\beta)}$ will correspond to the $L^{2}$ orthonormal
basis. Equation (48) also allows us to conclude that the spectrum of $K$
consists of eigenvalues $\kappa_{n,\ell}$, independent of $m$ and given by the
explicit formula

$\kappa_{n,\ell}=\frac{P_{n}^{(\alpha,\beta)}\Big{(}-1+\frac{2}{(N-1)^{2}}\Big{)}}{P_{n}^{(\alpha,\beta)}(1)}\bigg{(}-\frac{1}{N-1}\bigg{)}^{\ell}.$
(50)

While similar formulas derived from the above were useful in the asymptotic
$N$ case, a large amount of cancellation happens and makes the eigenvalues
take values much smaller than implied by naive bounds. Furthermore, it appears
that the eigenvalues display an explicit $\mod 3$ behavior. By this, we mean
that for each $\ell$, plotting the eigenvalues for $n\mod 3$ separately yields
3 eventually monotonous sequences, as we can observe bellow:

Figure 7: First 300 values of $\kappa_{n,5}$

Figure 8: First 300 values of $\kappa_{n,15}$

Figure 9: Overlap of the plots from $\ell=6$ to $10$

As these 3 different behaviours are included simultaneously in the sequences,
it further makes difficult the tasks of obtaining good uniform bounds.

One can also make use of the recurrence relations for Jacobi polynomials [10]
to obtain the same type of relations for the eigenvalues $\kappa_{n,\ell}$.
More precisely, for each $\ell$, we have that

$-\frac{1}{2}\kappa_{n,\ell}=\frac{n+\ell+\frac{1}{2}}{2n+\ell+1}\kappa_{n,\ell-1}+\frac{n+\frac{1}{2}}{2n+\ell+1}\kappa_{n+1,\ell-1}\
.$ (51)

Since we also know that

$\kappa_{n,0}=\frac{1}{n+1}(-1)^{n}\sin\bigg{(}\frac{n\pi}{3}\bigg{)},$ (52)

this provides a fast and efficient way to at least numerically compute all the
needed eigenvalues. We can, however, solve this recurrence relation in terms
of the $\ell=0$ values, even if the obtained representation isn’t very useful:

$\kappa_{n,\ell}=(-2)^{\ell}\sum_{j=0}^{\ell}\binom{\ell}{j}C_{n,\ell}^{(j)}\kappa_{n+j,0}\
,$ (53)

with

$C_{n,\ell}^{(j)}=\frac{\prod_{i=1}^{\ell}\\{n+i+1/2\\}}{\prod_{i=0,i\neq
j}^{\ell}\\{2n+i+j+2\\}}\ .$ (54)

One can see that the missing term in the denominator, when $i=j$, is exactly 2
times the numerator of $\kappa_{n+j,0}$. Unfortunately, this representation
shows no indication of decay, meaning that there is a large amount of
cancellation occurring (the terms are large, but very close in absolute
value). Collecting the sine terms may, however, shed some light on the $\mod
3$ question.

Thus, while explicit values of Jacobi polynomials are not trivial to compute,
the authors make use of uniform pointwise bounds obtained by Nevai, Erdelyi
and Magnus [9], a surprisingly deep result, to obtain monotonous bounds on
these eigenvalues:

###### Lemma 3.2

We have that

$\kappa^{2}_{n,\ell}\leq\frac{2e}{\pi}g_{1}(n,\ell,N)g_{2}(N)g_{3}(n,N)g_{4}(n,\ell,N)\
,$ (55)

with

$\displaystyle g_{1}(n,\ell,N)$
$\displaystyle=\bigg{(}\frac{4+\sqrt{9N^{2}-48N+65+4\ell^{2}+4\ell}}{3N+4n+2\ell-5}\bigg{)}\
,$ $\displaystyle g_{2}(N)$
$\displaystyle=\bigg{(}\frac{(N-1)^{2}}{N(N-2)}\bigg{)}^{(3N-7)/2}\ ,$
$\displaystyle g_{3}(n,N)$
$\displaystyle=\frac{\Gamma(n+1)\Gamma(\frac{3}{2}N-3)}{\Gamma(n+\ell+\frac{3}{2}N-\frac{5}{2})}\
,$ $\displaystyle g_{4}(n,\ell,N)$
$\displaystyle=\frac{(N-1)^{2}\Gamma(n+\ell+\frac{3}{2})\Gamma(\frac{3}{2}N-3)}{\Gamma(n+\ell+\frac{3}{2}N-\frac{5}{2})}\
.$ (56)

This bound can be enlarged to allow for the monotonous behaviour. This time,
however, when $N=3$, we are actually in a better position to further simplify
the achieved bound, yielding

$\kappa_{n,\ell}^{2}\leq\frac{8e}{3}\frac{\Gamma(n+1)}{\Gamma(n+3/2)}\frac{\Gamma(n+\ell+3/2)}{\Gamma(n+\ell+2)}\leq\frac{8e}{3}\frac{1}{(n+1)^{1/2}}\frac{1}{(n+\ell+3/2)^{1/2}}=:\hat{\kappa}_{n,\ell}^{2}\
.$ (57)

Now notice that, for all $\ell\geq 4$,

$(n+1)(n+\ell+3/2)\geq(n+\sqrt{\ell})^{2}\ ,$ (58)

so that if we define

$\widetilde{\kappa}_{n,\ell}=\sqrt{\frac{8e}{3}}\frac{1}{\sqrt{n+\sqrt{\ell}}}\
,$ (59)

we have $\hat{\kappa}_{n,\ell}\leq\widetilde{\kappa}_{n,\ell}$. Both upper
bounds will be useful later on. Inequality (57) is more precise, especially
when we are looking for bounds independent of $\ell$. Nevertheless, inequality
(59) has a far simpler representation and will be considered when we need to
take derivatives.

###  3.2 Anti-symmetric sector

Now that we have built the necessary tools, we first treat the simpler anti-
symmetric case, which can be done mostly through clever manipulation and
simple bounds. We reproduce here the approach done in [2], with some minor
numerical improvements.

###### Lemma 3.3

For $N=3$, the largest eigenvalue of $P^{(2)}$ on the orthogonal complement of
the symmetric sector is no larger than $0.729$. Thus, either we have
$\hat{\Delta}_{3,2}\geq 0.021$ or else the gap eigenfunction is symmetric.

Suppose $f$ is anti-symmetric and we assume, without loss of generality, that
$f(v)=\varphi(v_{1})-\varphi(v_{2})$, with $\varphi$ orthogonal to the
constants. Then, $P^{(2)}f=\lambda f$ becomes

$\frac{1}{N}w_{N,2}(v_{1})(1-K)\varphi(v_{1})-\frac{1}{N}w_{N,2}(v_{2})(1-K)\varphi(v_{2})=\lambda(\varphi(v_{1})-\varphi(v_{2}))$
(60)

Multiplying by $\varphi(v_{1})$ and integrating over $\mathcal{S}_{N}$, we
obtain

$\frac{1}{N}\int_{\mathcal{S}_{N}}w_{N,2}(v_{1})|(1-K)\varphi(v_{1})|^{2}d\sigma_{N}=\lambda\langle\varphi,(1-K)\varphi\rangle.$
(61)

From

$w_{N,2}(v_{k})=\frac{N}{N-1}-\frac{N}{(N-1)^{2}}|v_{k}|^{2},$ (62)

we can write (61) as

$\frac{1}{N-1}\langle\varphi,(1-K)^{2}\varphi\rangle-\frac{1}{(N-1)^{2}}\langle(1-K)\varphi,|v|^{2}(1-K)\varphi\rangle=\lambda\langle\varphi,(1-K)\varphi\rangle$
(63)

After decomposing $\sqrt{1-K}\varphi=\psi+\zeta$, where $\psi$ is orthogonal
to the constants, the components $v_{j}$ and $|v|^{2}$, we have that:

$\displaystyle\langle\varphi,(1-K)^{2}\varphi\rangle$
$\displaystyle=\langle\psi+\zeta,(1-K)(\psi+\zeta)\rangle$
$\displaystyle=\langle\psi,(1-K)\psi\rangle+(1+1/(N-1))\|\zeta\|^{2}_{2}$
$\displaystyle=\|\sqrt{1-K}\psi\|^{2}_{2}+\frac{N}{N-1}\|\zeta\|^{2}_{2}$ (64)

and

$\displaystyle\langle(1-K)\varphi,|v|^{2}(1-K)\varphi\rangle=$
$\displaystyle\langle\sqrt{1-K}(\psi+\zeta),|v|^{2}\sqrt{1-K}(\psi+\zeta)\rangle$
$\displaystyle=$
$\displaystyle\langle\sqrt{1-K}\psi,|v|^{2}\sqrt{1-K}\psi\rangle+\langle\sqrt{1-K}\zeta,|v|^{2}\sqrt{1-K}\zeta\rangle$
$\displaystyle+$ $\displaystyle
2\langle|v|\sqrt{1-K}\psi,|v|\sqrt{1-K}\zeta\rangle,$ (65)

since $\zeta$ is an eigenfunction of $K$ with eigenvalue $-1/(N-1)$, it is
also an eigenfunction of $\sqrt{1-K}$, with eigenvalue
$\sqrt{1+\frac{1}{N-1}}=\sqrt{\frac{N}{N-1}}$.

It then follows that

$\displaystyle\langle\sqrt{1-K}\varphi,|v|^{2}\sqrt{1-K}\varphi\rangle$
$\displaystyle\geq\langle\sqrt{1-K}\psi,|v|^{2}\sqrt{1-K}\psi\rangle+\langle\sqrt{1-K}\zeta,|v|^{2}\sqrt{1-K}\zeta\rangle$
$\displaystyle-2\||v|\sqrt{1-K}\psi\||_{2}\sqrt{\frac{N}{N-1}}\||v|\zeta\|_{2}$
$\displaystyle\geq\bigg{(}1-\frac{1}{t}\bigg{)}\langle\sqrt{1-K}\psi,|v|^{2}\sqrt{1-K}\psi\rangle+(1-t)\frac{N}{N-1}\||v|\zeta\|_{2}^{2},$
(66)

using the Cauchy-Schwarz and arithmetic-geometric mean inequalities.

Before we move on, it is useful to note that
$\||v|\zeta\|_{2}^{2}\geq\|\zeta\|^{2}_{2}$. To see this for $N=3$, write

$\zeta=\sum_{j=1}^{4}a_{j}\eta_{j},$ (67)

with $\eta_{j}(v)=\boldsymbol{e}_{j}\cdot v$ for $j=1,2,3$ and
$\eta_{4}=|v|^{2}-1$ as discussed in section 3.1.

From (67), a simple computation yields

$\||v|\zeta\|^{2}_{2}=\sum_{j=1}^{4}|a_{j}|^{2}\||v|\eta_{j}\|^{2}_{2}\ ,$
(68)

and while a bit more cumbersome, one can also show that

$\int_{\mathcal{S}_{3}}|v_{1}|^{4}d\sigma_{3}=\frac{5}{4}$ (69)

and

$\int_{\mathcal{S}_{3}}|v_{1}|^{6}d\sigma_{3}=\frac{7}{4}.$ (70)

Together, the equations above give us

$\||v|\eta_{j}\|^{2}_{2}=\frac{5}{4}\text{ for }j=1,2,3\text{ and
}\||v|\eta_{4}\|=1,$ (71)

which prove what we needed.

Finally, from the trivial bound
$\langle\sqrt{1-K}\psi,|v|^{2}\sqrt{1-K}\psi\rangle\leq(N-1)\|\sqrt{1-K}\psi\|^{2}_{2}$,
we have that, for all $0<t<1$,

$\langle\sqrt{1-K}\varphi,|v|^{2}\sqrt{1-K}\varphi\rangle\geq\bigg{(}1-\frac{1}{t}\bigg{)}(N-1)\|\sqrt{1-K}\psi\|^{2}_{2}+(1-t)\frac{N}{N-1}\|\zeta\|_{2}^{2}\
.$ (72)

Specializing to $N=3$, the second most negative eigenvalue of $K$ is
$-\frac{3}{8}$, which yields
$\|\sqrt{1-K}\psi\|^{2}_{2}\leq\frac{11}{8}\|\psi\|^{2}_{2}$. Combining all of
the above,

$\displaystyle\lambda(\|\psi\|^{2}_{2}+\|\zeta\|^{2}_{2})$
$\displaystyle\leq\frac{1}{2}\cdot\frac{11}{8}\|\psi\|^{2}_{2}+\frac{3}{4}\|\zeta\|^{2}_{2}-\frac{1}{2}\cdot\frac{11}{8}\bigg{(}1-\frac{1}{t}\bigg{)}\|\psi\|^{2}_{2}-\frac{1}{4}\frac{3}{2}(1-t)\|\zeta\|^{2}_{2}$
$\displaystyle=\frac{11}{16t}\|\psi\|_{2}^{2}+\frac{3}{8}\bigg{(}1+t\bigg{)}\|\zeta\|^{2}_{2}\
.$ (73)

With $t=0.943$, this yields $\lambda\leq 0.729$.

###  3.3 Symmetric sector

The remaining task is to bound the second largest eigenvalue of $P^{(2)}$ in
the symmetric sector. Assume $f$ is of the form

$f(\vec{v})=\sum_{j=1}^{N}\varphi(v_{j}).$ (74)

Then, we can write $P^{(2)}f=\lambda f$ as

$\frac{1}{N}\sum_{j=1}^{N}w_{N,2}(v_{j})(\varphi(v_{j})+(N-1)K\varphi(v_{j}))=\lambda\sum_{j=1}^{N}\varphi(v_{j}).$
(75)

Theorem 2.5 in [2] is here extremely important, as it allows us to treat the
above sum component wise:

$\frac{1}{N}w_{N,2}(v)(\varphi(v)+(N-1)K\varphi(v))=\lambda\varphi(v).$ (76)

By integrating against $\varphi(v)$, it follows that

$\frac{1}{N}\int_{\mathcal{S}_{N}}\varphi(v_{1})w_{N,2}(v_{1})(\varphi(v_{1})+(N-1)K\varphi(v_{1}))d\sigma_{N}=\lambda\|\varphi\|^{2}_{2}.$
(77)

As, once again,

$w_{N,2}(v)=\frac{N}{N-1}-\frac{N}{(N-1)^{2}}|v|^{2},$ (78)

and focusing on $N=3$, (77) becomes

$\displaystyle\lambda\|\varphi\|^{2}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{\mathcal{S}_{3}}\varphi(v_{1})(1-\tfrac{1}{2}|v|^{2})(1+2K)\varphi(v_{1})d\sigma_{N}$
(79) $\displaystyle=$
$\displaystyle\frac{1}{2}\langle(1-\tfrac{1}{2}|v|^{2})\varphi,(1+2K)\varphi\rangle.$

Since

$-1\leq\bigg{(}1-\frac{1}{2}|v|^{2}\bigg{)}\leq 1,\text{
}\bigg{\|}\bigg{(}1-\frac{1}{2}|v|^{2}\bigg{)}\varphi\bigg{\|}_{2}\leq\|\varphi\|_{2},$
(80)

we obtain the first bound

$\lambda\leq\frac{1}{2}\|1+2K\|_{2}.$ (81)

Unfortunately, this bound is useless as is, since the largest eigenvalue of
$K$, other than $1$, is $\frac{1}{2}$, giving us $\lambda\leq 1$. Thus, we
need to be more incisive on the spectrum of $K$ and its interaction with the
multiplication operator $(1-\frac{1}{2}|v|^{2})$.

Consider the operator

$M:=|v|^{2}(1+(N-1)K).$ (82)

Recalling the eigenbasis of $K$,

$g_{n,l,m}=h_{n,l}(|v|^{2})|v|^{l}\mathcal{Y}_{l,m}(v/|v|),$ (83)

we see that $M$ commutes with rotations and thus different angular momentum
sectors will be mutually orthogonal, so we can consider separately the
dependence in $\ell$ and $m$ (in fact, we will ignore $m$ completely, as
$g_{n,\ell,m}$, $g_{n,\ell,m^{\prime}}$ belong to the same $\ell$-eigenspace).
Thus, the quantities we are interested in estimating are

$\lambda_{\ell}:=\sup_{\|\varphi\|=1}\frac{1}{4}\langle(1-t)\varphi,(1+2K)\varphi\rangle,$
(84)

with the supremum being taken over the $\ell$ sector. Moreover, as
$h_{n,l}(|v|^{2})=P_{n}^{(\alpha,\beta)}(t)$, with $t=|v|^{2}-1$,
$\beta=l+\frac{1}{2}$, we see that $1-\frac{1}{2}|v|^{2}=\frac{1}{2}(1-t)$ in
the Jacobi variables, and we can exploit the properties of Jacobi polynomials.

The slow simultaneous decay in $n$ and $\ell$ we obtain from expressions like
(59) pose an obstacle for computing meaningful bounds. Thus, we break down the
problem over several ranges.

First, we look for a uniform bound given large enough $\ell>\ell^{*}$, but
such that $\ell^{*}$ is as small as possible.

For $5<\ell<\ell^{*}$, the eigenvalues are already small enough for (81) to
work as well, and we make use of a uniform bound in $n$ to reduce this problem
to explicitly computing a finite amount of these. Note that the recurrence
formula (recurrence) makes this fairly trivial and efficient.

Finally, we provide an argument based on linear algebra for $\ell\leq 5$,
where closed form expressions are simple enough for this.

####  3.3.1 Large $\ell$ case

Let $f_{n,\ell,m}$ be a re-scaling of $g_{n,\ell,m}$ such that the basis is
orthonormal, and fix $\ell$, $m$.

Writing $P_{n}:=P_{n}^{(\alpha,\beta)}$, the three-term relation for Jacobi
polynomials [10] is

$P_{n+1}(t)=(A_{n}t+B_{n})P_{n}(t)+C_{n}P_{n-1}(t)\ ,$ (85)

with

$\displaystyle A_{n}$
$\displaystyle=\frac{(2n+\alpha+\beta+1)(2n+\alpha+\beta+2)}{2(n+1)(n+\alpha+\beta+1)}\
,$ $\displaystyle B_{n}$
$\displaystyle=\frac{(\alpha^{2}-\beta^{2})(2n+\alpha+\beta+1)}{2(n+1)(n+\alpha+\beta+1)(2n+\alpha+\beta)}\
,$ $\displaystyle C_{n}$
$\displaystyle=\frac{(n+\alpha)(2n+\alpha+\beta+2)(n+\beta)}{(n+1)(n+\alpha+\beta+1)(2n+\alpha+\beta)}\
,$ (86)

which gives

$tP_{n}=(1/A_{n})P_{n+1}-(B_{n}/A_{n})P_{n}+(C_{n}/A_{n})P_{n-1}\ .$ (87)

As we are working with orthonormal functions, we have to re-scale the
equality. By the orthogonality relations for Jacobi polynomials,

$\langle
P_{n}^{(\alpha,\beta)},P_{n}^{(\alpha,\beta)}\rangle_{L^{2}(\mu^{(\alpha,\beta)})}=\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}\frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)n!}\delta_{nm}\
,$ (88)

we first compute

$\displaystyle\|P_{n}^{(\alpha,\beta)}\|_{2}^{2}=\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}\frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)n!}\
,$
$\displaystyle\frac{\|P_{n-1}^{(\alpha,\beta)}\|_{2}^{2}}{\|P_{n}^{(\alpha,\beta)}\|_{2}^{2}}=\frac{2n+\ell+2}{2n+\ell}\frac{n(n+\ell+1)}{(n+1/2)(n+\ell+1/2)}\
,$
$\displaystyle\frac{\|P_{n+1}^{(\alpha,\beta)}\|_{2}^{2}}{\|P_{n}^{(\alpha,\beta)}\|_{2}^{2}}=\frac{2n+\ell+2}{2n+\ell+4}\frac{(n+3/2)(n+\ell+3/2)}{(n+1)(n+\ell+2)}\
.$ (89)

Defining the scaling factors

$\displaystyle F_{n}^{-}$
$\displaystyle:=\frac{\|P_{n-1}\|}{\|P_{n}\|}=\sqrt{\frac{n(2n+\alpha+\beta+1)(n+\alpha+\beta)}{(n+\alpha)(n+\beta)(2n+\alpha+\beta-1)}}\
,$ $\displaystyle F^{+}_{n}$
$\displaystyle:=\frac{\|P_{n+1}\|}{\|P_{n}\|}=\sqrt{\frac{(n+\alpha+1)(n+\beta+1)(2n+\alpha+\beta+1)}{(n+1)(n+\alpha+\beta+1)(2n+\alpha+\beta+3)}}\
,$ (90)

the three-term relation for the multiplication operator in the orthonormal
Jacobi polynomials is

$tp_{n}=c_{n}p_{n-1}+a_{n}p_{n}+b_{n}p_{n+1},$ (91)

with

$a_{n}:=-\frac{B_{n}}{A_{n}}\text{, }b_{n}:=\frac{F^{+}_{n}}{A_{n}}\text{,
}c_{n}:=\frac{F^{-}_{n}C_{n}}{A_{n}}.$ (92)

Notice that, by the self-adjointness of multiplication by $t$, we must have
$c_{n}=b_{n-1}$. After some algebraic manipulation, we verify this and also
obtain

$b_{n,\ell}=2\sqrt{\frac{(n+1)(n+3/2)(n+\ell+3/2)(n+\ell+2)}{(2n+\ell+2)(2n+\ell+3)^{2}(2n+\ell+4)}},$
(93) $a_{n,\ell}=\frac{\ell(\ell+1)}{(2n+\ell+1)(2n+\ell+3)}.$ (94)

Evidently $a_{n,\ell}\geq\widetilde{a}_{n,\ell}$ where

$\widetilde{a}_{n,\ell}:=\left(\frac{\ell}{2n+\ell+3}\right)^{2}=\left(1-\frac{2n+3}{2n+\ell+3}\right)^{2}\
.$ (95)

It is clear that $\widetilde{a}_{n,\ell}$ is decreasing in $n$ and increasing
in $\ell$. It will also be useful to have a bound on $b_{n,\ell}$ that is
simpler to estimate. For this, notice that we can pair the ’close’ terms in
the numerator, and compare with those in the denominator. We then have

$\displaystyle b_{n,\ell}\leq$ $\displaystyle
2\sqrt{\frac{(n+3/2)^{2}(n+\ell+2)^{2}}{(2n+\ell+2)^{2}(2n+\ell+3)^{2}}}=\frac{2n+3}{2n+\ell+3}\frac{n+\ell+2}{2n+\ell+2}$
$\displaystyle=$
$\displaystyle\bigg{(}1-\frac{\ell}{2n+\ell+3}\bigg{)}\bigg{(}1-\frac{n}{2n+\ell+2}\bigg{)}\leq\bigg{(}1-\frac{\ell}{2n+\ell+3}\bigg{)}\bigg{(}1-\frac{n}{2n+\ell+3}\bigg{)}\
.$ (96)

We therefore define

$\widetilde{b}_{n,\ell}:=\bigg{(}1-\frac{\ell}{2n+\ell+3}\bigg{)}\bigg{(}1-\frac{n}{2n+\ell+3}\bigg{)}\
.$ (97)

In order to better understand the behaviour of these sequences, we treat them
as functions of continuous variables. Differentiating in $\ell$ we find

$\frac{{\partial}}{{\partial}\ell}\widetilde{b}_{n,\ell}=-\frac{(2\ell+6)n+3\ell+9}{(2n+\ell+3)^{3}}\
,$ (98)

so that $\widetilde{b}_{n,\ell}$ is decreasing in $\ell$. Differentiating in
$n$, we see that

$\frac{{\partial}}{{\partial}n}\widetilde{b}_{n,\ell}=\frac{2\ell^{2}+3\ell-(6n+9)}{(2n+\ell+3)^{3}}\
,$ (99)

so that $\widetilde{b}_{n,\ell}$ is increasing in $n$ for

$n\leq\frac{2\ell^{2}+3\ell-9}{6}\ .$

Suppose now $\varphi$ belongs to the $(\ell,m)$ angular momentum sector, so
that we can expand it as

$\varphi=\sum_{n=0}^{+\infty}\varphi_{n}f_{n,\ell,m}\ .$ (100)

Observing that

$\displaystyle(1+2K)\varphi=\sum_{n=1}^{+\infty}(1+2k_{n,\ell})\varphi_{n}f_{n,\ell,m}\
,$ (101)
$\displaystyle(1-t)\varphi=\sum_{n=1}^{+\infty}(b_{n-1}\varphi_{n-1}+(1-a_{n})\varphi_{n}+c_{n+1}\varphi_{n+1})f_{n,\ell,m}\
,$ (102)

we finally derive

$\displaystyle\langle(1-t)\varphi,(1+2K)\varphi\rangle=$
$\displaystyle-\sum_{n=1}^{+\infty}b_{n-1,\ell}(1+2k_{n,\ell})\varphi_{n-1}\varphi_{n}$
$\displaystyle+\sum_{n=1}^{+\infty}(1-a_{n,\ell})(1+2k_{n,\ell})\varphi_{n}^{2}$
$\displaystyle-\sum_{n=1}^{+\infty}b_{n,\ell}(1+2k_{n,\ell})\varphi_{n+1}\varphi_{n}\
.$ (103)

Taking the supremum on each term and using Cauchy-Schwarz on the
$\varphi_{n}\varphi_{n\pm 1}$ terms yields, and recalling the upper bound (59)
$\kappa_{n,\ell}\leq\widetilde{\kappa}_{n,\ell}$,

$\displaystyle\lambda_{\ell}$
$\displaystyle\leq\frac{1}{4}\bigg{[}2\sup_{n}\\{b_{n,\ell}(1+2\kappa_{n,\ell})\\}+\sup_{n}\\{(1-a_{n,\ell})(1+2\kappa_{n,\ell})\\}\bigg{]}$
$\displaystyle\leq\frac{1}{4}\bigg{[}2\sup_{n}\\{\widetilde{b}_{n,\ell}(1+2\widetilde{\kappa}_{n,\ell})\\}+\sup_{n}\\{(1-\widetilde{a}_{n,\ell})(1+2\widetilde{\kappa}_{n,\ell})\\}\bigg{]}.$
(104)

We can then turn our attention to the sequences
$\widetilde{b}_{n,\ell}(1+2\widetilde{\kappa}_{n,\ell})$ and
$(1-\widetilde{a}_{n,\ell})(1+2\widetilde{\kappa}_{n,\ell})$.

###### Proposition 3.1

For all $\ell\geq 70$, $\lambda_{\ell}\leq 0.73016$.

Proof –  We first compute the derivative in $n$ of
$(1-\widetilde{a}_{n,\ell})(1+2\widetilde{\kappa}_{n,\ell})$, which we find

$\frac{4\ell^{2}}{3}\frac{4\sqrt{6e}+3\sqrt{n+\sqrt{\ell}}}{(2n+\ell+3)^{3}\sqrt{n+\sqrt{\ell}}}-\frac{2}{3}\frac{\left(1-\frac{k^{2}}{(2n+\ell+3)^{2}}\right)\sqrt{6e}}{(n+\sqrt{\ell})^{3/2}}\
.$

which has the same sign as

${2\ell^{2}}\frac{4\sqrt{6e}+3\sqrt{n+\sqrt{\ell}}}{(2n+\ell+3)^{3}}-\frac{\left(1-\frac{k^{2}}{(2n+\ell+3)^{2}}\right)\sqrt{6e}}{n+\sqrt{\ell}}\
.$ (105)

Now assume $n\geq\ell^{3/2}$. Then, we have that
${2\ell^{2}}\leq(2n+\ell+3)^{4/3}$,
$\sqrt{n+\sqrt{\ell}}\leq(2n+\ell+3)^{1/2}$, and

$\frac{k^{2}}{(2n+\ell+3)^{2}}\leq\frac{1}{4\ell}\ .$

Thus, when $n\geq\ell^{3/2}$, the quantity in (105) is bounded above by

$\frac{4\sqrt{6e/\sqrt{\ell}}+3}{(n+\sqrt{\ell})^{7/6}}-\frac{\left(1-\frac{1}{4\ell}\right)\sqrt{6e}}{n+\sqrt{\ell}}\
.$

Hence our derivative is negative when $n\geq\ell^{3/2}$ and

$(n+\sqrt{\ell})^{1/6}\geq\frac{4\sqrt{6e/\sqrt{\ell}}+3}{\left(1-\frac{1}{4\ell}\right)\sqrt{6e}}\
.$

For $\ell=70$, this is satisfied for $n\geq 86$, and $70^{3/2}\leq 586$.
Hence, for $\ell=70$,
$(1-\widetilde{a}_{n,\ell})(1+2\widetilde{\kappa}_{n,\ell})$ is decreasing in
$n$ for al $n\geq 586$, and these values of $n$ and $\ell$ is it less than
$1.22$. Checking the finitely many remaining cases, one finds that the maximum
occurs at $n=66$, and is less than $1.4351$.

Proceeding similarly for $2\widetilde{b}_{n,\ell}(1+2\tilde{\kappa}_{n,\ell})$
we find the for $\ell=70$, the maximum occurs at $n=53$, and is no more than
$1.4855$. Altogether we find,

$\lambda_{70}\leq\frac{1.4855+1.4351}{4}\leq 0.73016\ .$

Since it is clear that $b_{n,\ell}(1+2\tilde{\kappa}_{,\ell})$ and
$(1-a_{n,\ell})(1+2\tilde{\kappa}_{n,\ell})$ are decreasing in $\ell$, this
proves that $\lambda_{\ell}\leq 0.73016$ for all $\ell\geq 70$.

####  3.3.2 Small $\ell$

By inequality (59), and ignoring the dependence in $\ell$, we see that for
$n\geq 151$ we have $\kappa_{n,\ell}\leq 0.23$. Numerically checking, this
holds as well for all $6\leq\ell\leq 50$ and $0\leq n\leq 151$. Thus, for
$5<\ell<70$, the norm of $K_{\ell}$, the restriction of $K$ to the $\ell$
angular momentum sector satisfies

$\|K_{\ell}\|\leq 0.23\ .$

It follows from (81) that for $5<\ell<70$,

$\lambda_{\ell}\leq\frac{1}{2}(1+2\cdot 0.23)=0.73\ .$

Taking into consideration the $\ell$ dependency can also help reduce the
amount of cases needed to check.

This leaves only the $\ell=0,1,2,3,4,5$ cases. Below, we provide a computation
for $\ell=0$, making use of the explicit expressions obtained via (51) and
(52). The same computation can be carried out analogously for the remaining
cases, at the expense of explicitly computing the three term relation factors
and the eigenvalues, and taking an appropriate choice of index for a break
point.

For $\ell=0$, the operator
$\varphi(v)\mapsto(1-\tfrac{1}{2}|v|^{2})\varphi(v)$ is given in the Jacobi
basis by the matrix

$X:=\frac{1}{4}\left[\begin{array}[]{cccccc}\phantom{-}2&-1&\phantom{-}0&\phantom{-}0&\phantom{-}0&\dots\\\
-1&\phantom{-}2&-1&\phantom{-}0&\phantom{-}0&\dots\\\
\phantom{-}0&-1&\phantom{-}2&-1&\phantom{-}0&\dots\\\
\phantom{-}0&\phantom{-}0&-1&\phantom{-}2&-1&\dots\\\
\phantom{-}0&\phantom{-}0&\phantom{-}0&-1&\phantom{-}2&\dots\\\
\phantom{-}\vdots&\phantom{-}\vdots&\phantom{-}\vdots&\phantom{-}\vdots&\vdots&\ddots\end{array}\right]\
.$

Let $Y$ be the diagonal matrix whose $n$th diagonal entry is
$1+2\kappa_{n+1,0}$. The upper $5\times 5$ block of $Z:=\frac{1}{2}(YX+XY)$ is

$\left[\begin{array}[]{ccccc}\phantom{-}\frac{1}{2}&-\frac{5}{16}&\phantom{-}0&\phantom{-}0&\phantom{-}0\\\
-\frac{5}{16}&\phantom{-}\frac{3}{4}&-\frac{21}{80}&\phantom{-}0&\phantom{-}0\\\
\phantom{-}0&-\frac{21}{80}&\phantom{-}\frac{3}{10}&-\frac{1}{5}&\phantom{-}0\\\
\phantom{-}0&\phantom{-}0&-\frac{1}{5}&\phantom{-}\frac{1}{2}&-\frac{2}{7}\\\
\phantom{-}0&\phantom{-}0&\phantom{-}0&-\frac{2}{7}&\phantom{-}\frac{9}{14}\end{array}\right]\
.$

According to Maple, the largest eigenvalue of this matrix is less than
$1.0412$. We also have that

$Z_{6,5}=Z_{5,6}=-\frac{57}{224}\ .$

Let $\varphi$ be any unit vector in $\ell^{2}$, and let $\eta$ be the vector
with $\eta_{j}=\varphi_{j}$ for $j=1,2,3,4$, and $\eta_{j}=0$ for $j>4$. Let
$\zeta:=\varphi-\eta$. Then

$\langle\varphi,Z\varphi\rangle=\langle\eta,\widehat{Z}\eta\rangle+\langle\zeta,Z\zeta\rangle-\frac{57}{224}2\varphi_{4}\varphi_{5}\
.$

We have

$\langle\eta,\widehat{Z}\eta\rangle\leq 1.0412\|\eta\|^{2}\ ,$

and by the Schwarz argument, since the largest relevant eigenvalue of $K$ is
$\tfrac{1}{10}$,

$\langle\zeta,Z\zeta\rangle\leq\frac{6}{5}\|\zeta\|^{2}\ .$

Thus, $\langle\varphi,Z\varphi\rangle$ is less than the largest eigenvalue of

$\left[\begin{array}[]{cc}1.0412&-\frac{57}{224}\\\
-\frac{57}{224}&\frac{6}{5}\end{array}\right]\ ,$

which is no more than $1.388$. Dividing by $2$, we get our bound on the top
eigenvalue in the $\ell=0$ sector, $0.694$.

##  4 Conclusion

Comparing the bounds obtained in each sector, theorem 2.1 follows. Obviously,
the result is only as good as the worst bound we obtain. While a bound of
$0.02$ is certainly achievable as we have seen, it is not clear if this method
could yield a much better result.

Nevertheless, we appear to be off by one order of magnitude, which is decent
enough, and we saw how the erratic behaviour of the eigenvalues
$\kappa_{n,\ell}$ poses an obstacle towards this.

Thus, a future goal would be proving entropy production inequalities with non-
degenerate ’constants’, first for the conjugate process and possibly for the
Kac process itself following Villani’s work. This is out of reach at the
moment, and careful geometrical analysis is needed due to the state space
constraints.

It is possible, however, that these quantities can be related analogously to
their respective spectral gaps, and new insights can be obtained this way.

##  Appendix A Small $\ell$ computations

For $\ell=1,...,5$, we may proceed exactly the same way as for the $\ell=0$
case. Recall that, through the recurrence relation (51), we are able to obtain
explicit formulas for eigenvalues in each $\ell$ sector. For instance, we can
always write $\kappa_{n,\ell}$ as a linear combination of $\kappa_{n,0}$,
$\kappa_{n+1,0}$, …, $\kappa_{n+\ell,0}$, and we see that

$\displaystyle\kappa_{n,1}$ $\displaystyle=-(\kappa_{n,0}+\kappa_{n+1,0})\ ,$
(106) $\displaystyle\kappa_{n,2}$
$\displaystyle=\frac{n+5/2}{n+2}\kappa_{n,0}+\kappa_{n+1,0}+\frac{n+3/2}{n+2}\kappa_{n+2,0}\
.$ (107)

Once again, these simplify if we consider separate cases given $n\mod 3$, and
analogous formulas hold for $\ell=3,4,5$, but with 4, 5 and 6 terms
respectively.

We begin by computing the upper $5\times 5$ block of the $Z$ matrix, as well
as its largest eigenvalue, and note the remainder entries that separate it
from the lower block (the coefficient of the $\varphi_{5}\varphi_{6}$ term).
Finally, we compute the largest relevant eigenvalue $\kappa_{n,\ell}$ for the
Schwarz argument.

This procedure is represented in the table below, with fractions denoting the
exact values and decimals as close upper bounds.

$\ell$ | Upper block bound | Remainder | Largest $\kappa_{n,\ell}$ | Eigenvalue bound
---|---|---|---|---
0 | 1.0412 | -57/224 | 1/10 | 0.694
1 | 0.946 | -0.2729 | 1/8 | 0.7052
2 | 0.895 | -0.2084 | 0.12 | 0.669
3 | 0.81 | -0.254 | 0.105 | 0.667
4 | 0.784 | -0.2275 | 0.12 | 0.6671
5 | 0.754 | -0.206 | 0.1 | 0.6403

##  Appendix B Python code for the Monte Carlo simulations

We leave in this appendix the Python code that was used for the Monte Carlo
simulations of the $\alpha=2$ case, as the only difference in the Maxwellian
case are the constant jump rates $1/N$. Please note that the code is not
optimized and vast memory and run time improvements could be made. The helper
functions are defined as close as possible to our theoretical conventions and
notations.

⬇

N = 3

def norm_2(v):

return sum([x**2 for x in v])

def exp_sample(lamb=1):

u = np.random.rand()

X = -np.log(1-u)/lamb

return X

def lamb(v):

l = [0]*N

for i in range(N):

l[i] = (N**2 - (1 + norm_2(v[i]))*N)/( N*(N-1)**2 )

return l

def sample_jump_time(v):

params_ = lamb(v)

jump_time = exp_sample( params_[0] )

jump_index = 0

for i in range(1,N):

candidate = exp_sample( params_[i] )

if candidate < jump_time:

jump_time = candidate

jump_index = i

return jump_time, jump_index

def beta(v):

return np.sqrt(3/2 * (1-norm_2(v)))

def conditional_jump_to_uniform3(vel,index):

new_vel = vel

v = vel[index]/np.sqrt(2)

y1 = np.random.normal(size=3)

y2 = -y1

scale = (np.sqrt(norm_2(y1)+norm_2(y2)))/np.sqrt(2)

y1 = y1/scale

y2 = y2/scale

if index == 0:

new_vel[1] = beta(v)*y1 - v/np.sqrt(2)

new_vel[2] = beta(v)*y2 - v/np.sqrt(2)

elif index == 1:

new_vel[0] = beta(v)*y1 - v/np.sqrt(2)

new_vel[2] = beta(v)*y2 - v/np.sqrt(2)

elif index == 2:

new_vel[0] = beta(v)*y1 - v/np.sqrt(2)

new_vel[1] = beta(v)*y2 - v/np.sqrt(2)

return new_vel

n_sim = 1000000

sim_count = 0

frames = [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]

first_sample = []

first_sample2 = []

first_sample3 = []

final_sample = [[] for x in range(12)]

final_sample2 = [[] for x in range(12)]

final_sample3 = [[] for x in range(12)]

index_dist = [[] for x in range(12)]

while sim_count < n_sim:

r = np.random.rand()

u = np.random.rand()

while u > (1-r):

r = np.random.rand()

u = np.random.rand()

v = np.random.normal(size=3)

v = r*v/np.sqrt(norm_2(v))

first_sample.append(2*r)

# initialize the remaining coordinates according to the parametrization

v2 = np.random.normal(size=3)

v3 = -v2

scale = (np.sqrt(norm_2(v2)+norm_2(v3)))/np.sqrt(2)

v2 = v2/scale

v3 = v3/scale

vel = [np.sqrt(2)*v,beta(v)*v2 - v/np.sqrt(2), beta(v)*v3 - v/np.sqrt(2)]

first_sample2.append(norm_2(vel[1]))

first_sample3.append(norm_2(vel[2]))

t = 0

for i in range(len(frames)):

T = frames[i]

while t<T:

tau, index = sample_jump_time(vel)

index_dist[i].append(index)

t += tau

vel = conditional_jump_to_uniform3(vel, index)

final_sample[i].append(np.sqrt(norm_2(vel[0])))

final_sample2[i].append(np.sqrt(norm_2(vel[1])))

final_sample3[i].append(np.sqrt(norm_2(vel[2])))

sim_count += 1

entropy = []

first_dist = [x/2 for x in first_sample]

counts, bins, bars = plt.hist(first_dist, bins=100, density=True)

ent = sum([counts[i]*np.log(counts[i]/(5.09295817894*(bins[i])**2 *

(1-bins[i]**2)**(1/2)))/100 for i in range(1,len(counts))])

entropy.append(ent)

x = np.linspace(0,1)

plt.plot(x, 5.09295817894*x**2 * (1-x**2)**(1/2))

plt.show()

for i in range(12):

final_dist = [x/np.sqrt(2) for x in final_sample[i]]

counts, bins, bars = plt.hist(final_dist, bins=100, density=True)

ent = [counts[i]/100 * np.log( counts[i] /

((bins[i])**2*(1-bins[i]**2)**(1/2)*5.09295817894))

for i in range(1,len(counts))]

total_ent = sum(ent)

entropy.append(total_ent)

plt.plot(x, 5.09295817894*x**2 * (1-x**2)**(1/2))

plt.show()

frames0 = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]

plt.scatter(frames0, y = np.log(entropy))

## References

  * [1] M. Kac., Foundations of kinetic theory, Proc. 3rd Berkley symp. Math. Stat. Prob., J. Neyman, ed. Univ. of California, vol 3, (1956) 171-197
  * [2] E. A. Carlen, M. C. Carvalho and M. Loss, Spectral gaps for reversible Markov processes with chaotic invariant measures: the Kac process with hard sphere collisions in three dimensions, arXiv:1812.03874 [math-ph]
  * [3] E. A. Carlen, J. S. Geronimo and M. Loss, Determination of the spectral gap in the Kac model for physical momentum and energy conserving collisions, arXiv:0705.3729 [math-ph]
  * [4] E. A. Carlen, J. S. Geronimo and M. Loss, On the Markov sequence problem for Jacobi polynomials, Advances in Mathematics, Volume 226, Issue 4, (2011) 3426-3466.
  * [5] E. Janvresse, Spectral Gap for Kac’s model of Boltzmann Equation, Annals. of Prob., Volume 29, (2001) 288-304.
  * [6] C. Villani, Cercignani’s Conjecture is Sometimes True and Always Almost True, Comm. Math. Phys., Volume 234, (2003) 455–490.
  * [7] F. Barthe, D. Cordero-Erausquin and B. Maurey, Entropy of spherical marginals and related inequalities, Jour. de Math., Pures et Appliquées, Volume 86, (2006) 89–99.
  * [8] E. A. Carlen, E. H. Lieb, and M. Loss, A sharp analog of Youngs inequality on SN and related entropy inequalities, J. Geom. Anal., Volume 14, (2004) 487–520.
  * [9] P. Nevai, T. Erdélyi and A. Magnus, Generalized Jacobi weights, Christoffel functions and Jacobi polynomials, SIAM Jour. Math. Anal., Volume 25, (1994) 602–614.
  * [10] G. Szegö, Orthogonal Polynomials, American Mathematical Society Colloquium Publications, ,Volume 23, AMS, Providence R.I. 193


# Measurement of forward photon production cross-section in pp collisions at
$\sqrt{s}$ = 510 GeV with RHICf detector

O. Adriani  E. Berti  L. Bonechi  R. D’Alessandro  Y. Goto  B. Hong  Y. Itow
K. Kasahara  J. H. Lee  T. Ljubicic  M. H. Kim  H. Menjo
<EMAIL_ADDRESS>I. Nakagawa  A. Ogawa  K. Ohashi  R. Pak  T. Sako
N. Sakurai  K. Sato  R. Seidl  K. Tanida  S. Torii  A. Tricomi

###### Abstract

This study reported the inclusive differential production cross-section of
photons in six pseudorapidity regions: $6.1\,<\eta\,<6.5$, $6.5\,<\eta\,<7.0$,
$7.0\,<\eta\,<7.5$, $7.5\,<\eta\,<8.0$, $8.0\,<\eta\,<8.5$, and $\eta\,>8.5$,
measured through the RHICf experiment with pp collisions at $\sqrt{s}$ = 510
GeV conducted in June 2017. In addition, the cross-sections in the three
regions of the $x_{\mathrm{F}}$-$p_{\mathrm{T}}$ phase space coverage that are
same as those of the LHCf results at $\sqrt{s}$ = 7 and 13 TeV were obtained
and compared. Considering the uncertainties, the results were observed to be
consistent with both the Feynman scaling law and the model predictions of
EPOS-LHC, QGSJET-II-04, Sibyll 2.3d, and DPMjet-III 2019.1, although certain
models exhibited weak collision energy dependencies.

###### keywords:

RHIC , UHECR , Hadronic interaction

###### PACS:

0000 , 1111

###### MSC:

0000 , 1111

††journal: Physics Letters B

[a]organization=Department of Physics and Astronomy, University of Florence,
city=Sesto Florentino, postcode=I-50019, country=Italy [b]organization=INFN
Section of Florence, city=Sesto Fiorentino, postcode=I-50019, country=Italy
[c]organization=RIKEN Nishina Center for Accelerator-Based Science, city=Wako,
postcode=351-0198, state=Saitama, country=Japan [d]organization=RIKEN BNL
Research Center, Brookhaven National Laboratory, city=New York,
postcode=11973-5000, country=USA [e]organization=Korea University, city=Seoul,
postcode=02841, country=Korea [f]organization=Institute for Space-Earth
Environmental Research, Nagoya University, city=Nagoya, postcode=464-8602,
state=Aichi, country=Japan [g]organization=Kobayashi-Maskawa Institute for the
Origin of Particles and the Universe, Nagoya University, city=Nagoya,
postcode=464-8602, state=Aichi, country=Japan [h]organization=Shibaura
Institute of Technology, city=Minuma-ku, postcode=337-8570, state=Saitama,
country=Japan [i]organization=Brookhaven National Laboratory, city=Upton,
postcode=11973-5000, state=New York, country=USA [j]organization=Institute for
Cosmic Ray Research, University of Tokyo, city=Kashiwa, postcode=277-8582,
state=Chiba, country=Japan [k]organization=Tokushima University,
city=Tokushima, postcode=770-8051, state=Tokushima, country=Japan
[l]organization=Advanced Science Research Center, Japan Atomic Energy Agency,
city=Tokai-mura, postcode=319-1195, state=Ibaraki, country=Japan
[m]organization=RISE, Waseda University, city=Shinjuku, postcode=162-0044,
state=Tokyo, country=Japan [n]organization=Department of Physics and
Astronomy, University of Catania, city=Catania, postcode=I-95123,
country=Italy [o]organization=INFN Section of Catania, city=Catania,
postcode=I-95123, country=Italy [p]organization=CSFNSM, city=Catania,
postcode=I-95123, country=Italy

##  1 Introduction

A precise understanding of hadronic interactions over a wide energy range is
among the most important parameters for solving the problem of ultra-high-
energy cosmic rays (UHECRs). The Pierre Auger Observatory and Telescope Array
observe UHECRs with energies above $10^{18}$ eV using the extensive air shower
technique, and they have reported their observation results with high-
statistics data [1, 2, 3]. However, their observation results of chemical
composition, which is an important parameter for studying the acceleration
mechanism of UHECRs at the sources, strongly depend on the selection of the
hadronic interaction model used in their air shower simulation. The
uncertainty of the models is primarily a result of the lack of calibration
data from the accelerator experiments. Particles produced in the very forward
region of interactions carry most of the projectile particle energy, and are
crucial to air shower development. The LHC forward (LHCf) experiment measured
the production cross-sections of very forward photons, $\pi^{0}$, and neutrons
at LHC [4, 5, 6, 7, 8, 9]. The maximum center-of-mass collision energy of
$\sqrt{s}$ = 13 TeV corresponds to $10^{17}$ eV in cosmic ray interactions. In
addition to their measurement, experimental data at lower collision energies
are necessary because an air shower is caused by complex interactions with
various interaction energies ranging from GeV to the UHECR energy.

The RHIC forward (RHICf) experiment was designed to perform a measurement
similar to the LHCf at RHIC located in Brookhaven National Laboratory and
investigate the collision-energy scaling law of forward particle production
proposed in Ref. [10], so-called the Feynman scaling, via comparisons of the
LHCf results. Whereas the LHCf reported the scaling of forward $\pi^{0}$
production between $\sqrt{s}$ 2.76 and 7 TeV [7], the scaling can be tested by
a LHCf–RHICf comparison in more than 10 times wider energy range than that.
This study analyzed the data acquired with the RHICf detector in June 2017
with $pp$ collisions at $\sqrt{s}$ = 510 GeV and presented the measured
inclusive differential production cross-section of photons in the six
pseudorapidity regions covering $\eta$ $>$ 6.1 as a function of Feynman $x$
(hereafter $x_{\mathrm{F}}$111It is calculated as $x_{\mathrm{F}}$=
$E_{\gamma}/E_{\mathrm{beam}}$, where $E_{\gamma}$ and $E_{\mathrm{beam}}$ are
the energies of photons and beams (255 GeV). ). Furthermore, the results in
three additional regions covering the same $x_{\mathrm{F}}$-$p_{\mathrm{T}}$
phase spaces as the LHCf result at $\sqrt{s}$ = 7 and 13 TeV were obtained and
subsequently compared to test the Feynman scaling.

##  2 Detector and Data

The RHICf detector comprises two compact sampling and positioning calorimeters
composed of tungsten, 16 GSO scintillator plates, and 4 XY horoscopes of GSO
bars. Two calorimeter tower have diamond-shaped acceptances transverse to the
beam directions of 20 mm $\times$ 20 mm and 40 mm $\times$ 40 mm, which are
referred to as the TS and TL, respectively. It is the former LHCf-Arm1
detector [11] and its performance has been well studied via beam tests using
electron and proton beams at CERN-SPS [12]. The energy and position
resolutions of photons with energies greater than 30 GeV are better than 5%
and 0.2 mm, respectively, [13].

The detector was installed 18 m west of the Solenoidal Tracker at the RHIC
(STAR) interaction point (IP) during the low-luminosity operation period from
24th to 27th June 2017. Because of the magnetic field of the dipole magnet
located between the IP and detector, only neutral particles, photons, and
neurons hit the detector. Further, the acceptance at the location was
restricted because of the shadow of the beam pipe structure located in front
of the detector, which corresponded to a pseudorapidity coverage greater than
6.1. The detector was moved vertically to cover full acceptance, and the data
acquisition was performed in three detector positions: BOTTOM ($\Delta y$ =
-47.4 mm), MIDDLE ($\Delta y$ = 0 mm), and TOP ($\Delta y$ = +24.04 mm), where
$\Delta y$ corresponds to the relative vertical position to MIDDLE.

In this analysis, the data acquired with the two trigger modes, Shower and
High-EM [13] was used. The Shower trigger was the baseline trigger designed to
detect both electromagnetic and hadronic showers induced by photons and
neutrons, and a large prescale factor of 8–30 was applied. The trigger
efficiency was 100% for photons with energies greater than 20 GeV. Further,
High-EM triggers were generated only for events with high-energy deposits in
the 4th layer of the GSO scintillator to enrich high-energy photon events in
the sample. The prescale factor for the High-EM trigger was approximately 10
times smaller than that for the Shower trigger. In addition, the trigger
efficiencies of the High-EM trigger were 96% for photons with energies over
100 GeV.

The beam energy was 255 GeV, and the proton beams were radially polarized. The
direction of beam polarization was rotated by 90 degrees from vertical
polarization around the beam direction. The luminosity during the operation
period was measured at $L$ = $(0.7-1.5)\,\times
10^{31}\mathrm{cm^{-2}{s^{-1}}}$ using the coincidence rate of signals from
the STAR-ZDC detectors [14] located on both the east and west sides of the IP.
The recorded integral luminosities of the Shower (High-EM) trigger sample were
4.5 (27.6), 12.0 (96.0), and 20.2 (120.0) $\mathrm{nb^{-1}}$ for the BOTTOM,
MIDDLE, and TOP detector positions, respectively.

##  3 Analysis

###  3.1 Event reconstruction and selection

For the analysis, the event reconstruction algorithm reported in Ref. [6] was
employed; further details can be found in Ref. [13]. The criteria in this
analysis were reoptimized for the RHICf configuration using a detector
simulation implemented with the GEANT4 library. Events that satisfied the
criteria of particle identification (PID) for photons and a single hit were
selected, corresponding to only one particle hit in a calorimeter tower.

The energy scale of the calorimeters was calibrated using the factor estimated
in a study using $\pi^{0}$ events where photon pairs were observed in the
calorimeters. Thereafter, the invariant mass was calculated from the energies
and hit positions of a photon pair, assuming a decay vertex at the IP. The
factor was calculated to shift the peak of the mass distribution to the rest
mass of $\pi^{0}$.

Analysis regions were defined for the production cross-section measurement
through the division of the pseudorapidity coverage at the location from 6.1
to infinity by six; $6.1\,<\eta\,<6.5$, $6.5\,<\eta\,<7.0$, $7.0\,<\eta\,<7.5$
$7.5\,<\eta\,<8.0$, $8.0\,<\eta\,<8.5$, and $\eta\,>8.5$. The region is a
sector-shape area defined by the range of the distance from the beam center222
The beam center was defined as the projection of the beam direction at the IP
on the detector surface, which was measured using data [13]. $R$ and the
azimuthal angle $\phi$, whereas the region for $\eta$ $>$ 8.5 is a disk shape
with $R$ $<$ 7.25 mm. In addition to these regions, three analysis regions
were defined to test the Feynman scaling law. Each region covered the same
$x_{\mathrm{F}}$-$p_{\mathrm{T}}$ phase space as one of LHCf measurements;
$p_{\mathrm{T}}$ $<$ 0.124 $x_{\mathrm{F}}$, 0.87 $x_{\mathrm{F}}$ $<$
$p_{\mathrm{T}}$ $<$ 1.04 $x_{\mathrm{F}}$ for $\eta\,>\,10.94$ and
$8.81\,<\,\eta\,<\,8.99$ at $\sqrt{s}$ = 7 TeV [5], respectively, and
$p_{\mathrm{T}}$ $<$ 0.23 $x_{\mathrm{F}}$ for $\eta\,>\,10.94$ at $\sqrt{s}$
= 13 TeV [6]. The $R$ and $\phi$ ranges of each region and the data set
(detector position and calorimeter tower) used for the analysis regions are
listed in Tab. 1. However, region $7.0\,<\eta\,<7.5$ was not covered by a
single dataset. Therefore, subregions corresponding to $7.00\,<\eta\,<7.25$
and $7.25\,<\eta\,<7.50$ were defined, which were combined to obtain the final
result.

The measurement performed discovered a large transverse single-spin asymmetry
of forward $\pi^{0}$ production using the same dataset [15]. Owing to the
$\pi^{0}$ asymmetry inducing dependency of the photon spectrum on the spin
direction of the incoming beams to the detector, the events were categorized
with the spin direction inward or outward of the bunch associated with each
event. The difference of the integral luminosity between the spin samples was
only 1–4%.

Table 1:  Analysis regions and data sets used for this analysis.  Region | R [mm] |  $\phi$ [degrees] | Dataset
---|---|---|---
$\eta\,>8.5$ |  $<$ 7.3 | 0 - 360 | BOTTOM-TL
$8.0\,<\eta\,<8.5$ | 7.3 - 12.0 | 50 - 130, 230 - 310 | BOTTOM-TL
$7.5\,<\eta\,<8.0$ | 12.0 - 19.7 | 75 - 105, 255 - 285 | BOTTOM-TL
$7.0\,<\eta\,<7.5$ |  |  |
$7.25\,<\eta\,<7.50$ | 19.7 - 25.3 | 75 - 105 | TOP-TS
$7.00\,<\eta\,<7.25$ | 25.3 - 32.5 | 80 - 100 | MIDDLE-TL
$6.5\,<\eta\,<7.0$ | 32.5 - 53.6 | 80 - 100 | MIDDLE-TL
$6.1\,<\eta\,<6.5$ | 53.6 - 79.9 | 80 - 100 | TOP-TL
$p_{\mathrm{T}}$$<$ 0.124 $x_{\mathrm{F}}$ |  $<$ 8.7 | 45 - 135, 225 - 315 | MIDDLE-TL
$p_{\mathrm{T}}$$<$ 0.23 $x_{\mathrm{F}}$ |  $<$ 16.1 | 70 - 110, 250 - 290 | MIDDLE-TL
0.87 $x_{\mathrm{F}}$$<$ $p_{\mathrm{T}}$$<$ 1.04 $x_{\mathrm{F}}$ | 60.7 - 72.8 | 80 - 100 | TOP-TL

###  3.2 Corrections

The final result of the inclusive differential production cross-section of
photons $\mathrm{d}\,\sigma_{\gamma}/\mathrm{d}\,x_{\mathrm{F}}$ was obtained
considering the number of events passing the event selection
$N^{\mathrm{single-\gamma}}_{spin}$ with a bin-by-bin correction procedure.

$\frac{\mathrm{d}\,\sigma_{\gamma}}{\mathrm{d}\,x_{\mathrm{F}}}=C^{\mathrm{geo.}}C^{c\tau}C^{\mathrm{MC}}\sum_{spin}\frac{C^{\mathrm{PID}}_{spin}(1-R^{BKG}_{spin})}{2\,{L_{spin}}\,}\dot{\frac{\Delta
N^{\mathrm{single-\gamma}}_{spin}}{\Delta x_{\mathrm{F}}}},$ (1)

where $C^{\mathrm{geo.}}$, $C^{c\tau}$, $C^{\mathrm{MC}}$,
$C^{\mathrm{PID}}_{spin}$, and $R^{BKG}_{spin}$ are the correction factors
described below, and $L_{spin}$ is the recorded integral luminosity of the
beams in the spin direction. Subsequently, the PID purity correction
$C^{\mathrm{PID}}_{spin}$ and the beam-gas background $R^{BKG}_{spin}$ were
estimated for each spin sample, whereas the other correction factors were
common. The difference of photon cross-sections between the two spin samples
was 20% at maximum, which originated from the $\pi^{0}$ asymmetry.
Consequently, upon combining the two spin samples, the effect of spin
asymmetry was expected to be negligible in the final result even considering
the systematic uncertainty of the integral luminosity.

  * 1.

Beam-gas background $R^{BKG}_{spin}$;
The background photons produced at collisions between the beam protons and
residual gas in the beam pipe were estimated using the events associated with
the non-colliding bunches of the incoming beam to the detector. Further, the
background fraction $R^{BKG}_{spin}$ was calculated for each spin-direction
sample by considering the difference in the total beam intensities between the
colliding and non-colliding bunches [13]. However, owing to the limited
statistics of non-colliding samples and little energy dependency, the
correction was applied as a photon energy-independent factor. The estimated
factors were 1.5 – 2.5%. In addition, the observed energy dependency was
considered a systematic uncertainty.

  * 2.

PID purity correction $C^{\mathrm{PID}}_{spin}$;
In this analysis, particle identification was performed by introducing a PID
estimator $L_{90\%}$, which was defined as the longitudinal depth at which the
energy deposit integral measured by the sampling layers reached 90% of the
total energy deposition. Subsequently, the purity of the photons in the
selected events for each energy bin was estimated from the template-fit
results of the $L_{90\%}$ distribution with MC distributions for photons and
hadrons [13]. The correction factor $C^{\mathrm{PID}}_{spin}$ (= 1/purity) was
typically 0.9 – 1.0.

  * 3.

MC based correction factor $C^{\mathrm{MC}}$;
This factor was introduced as an overall correction of several contributions
related to event reconstruction: inefficiencies of the trigger, PID selection
and single-hit selection, misreconstruction of multi-hit events (two particles
hit in a tower) as single-hit events, and recovery of photon yield in the
multi-hit events rejected by the single-hit selection. The total of these
contributions was estimated using the detector simulation with $pp$ event
generation by the interaction model of EPOS-LHC or QGSJET-II-04 as

$C^{\mathrm{MC}}=N^{\mathrm{\gamma}}_{\mathrm{MC,true}}/N^{\mathrm{single-\gamma}}_{\mathrm{MC,rec}},$
(2)

where $N^{\mathrm{\gamma}}_{\mathrm{MC,true}}$ is the number of photons
hitting the analysis region and $N^{\mathrm{single-\gamma}}_{\mathrm{MC,rec}}$
is the number of events reconstructed as single photons with the reconstructed
hit position in the region. It is noted that hadron contamination was not
considered in $N^{\mathrm{single-\gamma}}_{\mathrm{MC,rec}}$ because it should
be corrected using PID correction. Thus, the average of the results obtained
using the two interaction models was used as the correction factor. The
typical value of the factor was 1.1 – 1.2. The difference from one mainly
originates from the inefficiency of PID selection ( 0.1) and the recovery of
photon yield in multi-hit events (0.02 – 0.08). Moreover, the interaction-
model dependency of the factor is primarily caused by the difference in multi-
hit events, and it was assigned as systematic uncertainty.

  * 4.

Correction of long-life particle contribution $C^{c\tau}$;
The dominant source of photons measured by the detector is the $\pi^{0}$decay
produced at collisions, while decays of long-life particles such as $K^{0}$
and $\Lambda$ along the beam pipe located between the IP and the detector
contribute to the photon yield. The contribution of long-life particles was
estimated using an MC simulation by comparing the photon flux at the IP with
that after transportation from the IP to the detector location. It was removed
via the application of the correction factor. The contribution was relatively
large in low $x_{\mathrm{F}}$ regions, typically 8% of the total yield at
$x_{\mathrm{F}}$ approximately 0.1, whereas it was less than 1% in the high
$x_{\mathrm{F}}$ region of more than 0.5. The correction factors were
calculated as the average of the results obtained using the four interaction
models: EPOS-LHC [16], QGSJET-II-04 [17], Sibyll 2.3d [18], and DPMjet-III
2019.1 [19]. Furthermore, the systematic uncertainty of the factor was
estimated as the maximum deviation from the average of each model result.

  * 5.

Geometrical correction factor $C^{\mathrm{geo.}}$;
A correction factor was introduced to correct the limitation of the $\phi$
range, and it was calculated as $C^{\mathrm{geo.}}=360^{\circ}/\Delta\phi$.

###  3.3 Systematic uncertainties

The following contributions were considered as systematic uncertainties of the
measured production cross-section, in addition to those related to the
correction factors $R^{BKG}$, $C^{\mathrm{MC}}$, and $C^{c\tau}$ described in
Sec. 3.2). Fig. 1 shows the estimated systematic uncertainties for the
pseudorapidity region $8.0\,<\eta\,<8.5$ as a typical example.

Figure 1: Systematic uncertainties of the photon production cross-section
measurement for the pseudorapidity region $8.0\,<\eta\,<8.5$. The colored and
dashed lines indicate the estimated systematic uncertainties after
normalization based on the mean values of the experimental data. The black
line indicates the total systematic uncertainties, calculated as the quadratic
summations of all uncertainties.

  * 1.

Energy scale
Three components were considered as the systematic uncertainty of the energy
scale. The first was the stability of the energy scale during the operation,
which was monitored using the peak position of the $\pi^{0}$ mass on the
reconstructed mass distribution, and it was stable within $\pm 1\%$ [13]. The
second was the nonlinearity. It was observed that the peak mass value of
$\pi^{0}$ was shifted with increasing $\pi^{0}$ energy corresponding to the
total energy of the photon pair [13]. Assuming that the shift originates only
from the nonlinearity of the energy scale conservatively, the corrected photon
energy $E^{\prime}$ was obtained as
$E^{\prime}=(1.03-0.025\,(E/100\mathrm{GeV}))E$, where $E$ is the original
value of the reconstructed energy. The last factor is the nonuniformity of the
energy scale, which was evaluated from the consistency test of the spectra
obtained by dividing the $\phi$ range. The impact of each energy-scale
uncertainty on the final result was estimated by repeating the analysis by
artificially changing the energy scale. Owing to the steep slope of the
spectrum, large systematic uncertainties of over 50% were induced in the
highest energy bin.

  * 2.

PID selection
The experimental distribution of the PID estimator $L_{90\%}$ was not
perfectly reproduced by the template distributions obtained from the MC. The
systematic effect of the difference on the final result through the purity and
efficiency estimation of PID was estimated using a method similar to that used
in Ref. [5, 8]; Another template-fitting method that allowed artificial
displacement and widening of the template distributions was performed.
Consequently, the difference in the results between the methods was assigned
as systematic uncertainty.

  * 3.

Beam center
The effect of the uncertainty of the beam-center determination was evaluated
via repeating the analysis procedure by shifting the beam center by $\pm$1 mm
vertically or horizontally. The maximum deviation in these results from the
original results was assigned as the systematic uncertainty.

  * 4.

Luminosity
The uncertainty of the integral luminosity is $\pm 5\%$, which is due to the
uncertainty of the conversion factor from the ZDC signal rate to the absolute
luminosity determined from the Vernier scans using the method similar to that
described in Ref. [20].

##  4 Results

###  4.1 differential production cross-section

Fig. 2 presents the differential production cross-section of photons as a
function of $x_{\mathrm{F}}$,
$\mathrm{d}\,\sigma_{\gamma}/\mathrm{d}\,x_{\mathrm{F}}$, for the six
pseudorapidity intervals measured by the RHICf detector. The data points in
the $x_{\mathrm{F}}$ range below 0.45 were obtained from the Shower trigger
samples, while those in the higher $x_{\mathrm{F}}$ region were obtained from
the High-EM trigger samples because of large statistics. The colored lines in
Fig. 2 show the prediction of the interaction models, EPOS-LHC, QGSJET-II-04,
Sibyll 2.3d, DPMjet-III 2019.1, that were calculated using the event generator
CRMC version 1.8.0 [21]. Fig. 3 presents the ratios of these predictions to
the experimental results. EPOS-LHC and DPMjet-III 2019.1 were consistent with
the experimental result in the low $x_{\mathrm{F}}$ region below 0.6 and 0.3,
respectively, for any pseudorapidity regions, although it predicted larger
flux than the data in the higher $x_{\mathrm{F}}$ region. QGSJET-II-04 and
Sibyll 2.3d well reproduced the shape of $x_{\mathrm{F}}$ spectra in the
highest pseudorapidity regions; however, the slope of their spectra became
progressively softer and harder, respectively, than that of the data in the
lower pseudorapidity regions. They can be interpreted as the difference in
$p_{\mathrm{T}}$ distributions between the data and the model. Thus, these
features of the data-model comparisons are very similar to those at the LHC
energies reported by LHCf [4, 5, 6], thereby suggesting that the source of the
difference is common to wide collision energies.

Figure 2: Inclusive photon production cross-section measured by the RHICf
detector. Each figure represents that for one of the pseudorapidity regions:
$\eta$ $>$ 8.5, 8.0 $<$ $\eta$ $<$ 8.5, 7.5 $<$ $\eta$ $<$ 8.0, 7.0 $<$ $\eta$
$<$ 7.5, 6.5 $<$ $\eta$ $<$ 7.0, and 6.1 $<$ $\eta$ $<$ 6.5. The bars and
hatched areas correspond to the statistical and systematic uncertainties,
respectively. Colored lines indicate the MC predictions of DPMjet-III 2019.1,
EPOS-LHC, QGSJET-II-04, and Sibyll 2.3d.

Figure 3: Ratio of inclusive photon production cross-sections predicted by
hadronic interaction models to the experimental result. The bars and hatched
areas around one correspond to the normalized statistical and systematic
uncertainties of the data, respectively.

.

###  4.2 Collision-energy scaling

The measured cross-sections were compared with the results of LHCf at $pp$
collisions at $\sqrt{s}$ = 7 and 13 TeV to test the collision energy
dependency of forward photon production. Fig. 4 shows the results obtained for
the same coverage of $x_{\mathrm{F}}$-$p_{\mathrm{T}}$ phase space as that of
the LHCf results. In this comparison, an additional normalization factor
$1/\sigma_{\mathrm{inela}}$ was introduced, where $\sigma_{\mathrm{inela}}$
was the inelastic cross section. For normalization of the LHCf results, the
values of $\sigma_{\mathrm{inela}}$ = $(72.9\,\pm\,1.5)$ mb and
$(79.5\,\pm\,1.8)$ mb reported in Ref. [22, 23] were used333A slightly lower
value $\sigma_{\mathrm{inela}}$ = 71.5 mb was used in Ref. [5]. Their results
were rescaled using this ratio.. The $pp$ inelastic cross section at
$\sqrt{s}$ = 510 GeV was estimated to be $(48.3\,\pm\,1.1)$ mb using the fit
results of total and elastic cross-sections by the COMPETE [24] and TOTEM [23]
collaborations, respectively. The RHICf result was consistent with the LHCf
results within the uncertainties, and the ratios are shown in Fig. 5. The
values of the ratios were slightly lower than those in the middle panel, which
can be interpreted as the difference in analysis methods: corrections of
photon yield in the rejected multi-hit events and photons from long-life
particle decays were not performed in the LHCf 7 TeV result, whereas they were
performed in both the RHICf and LHCf 13 TeV results. The correction increased
the ratio by 0 – 15% in $x_{\mathrm{F}}$ below 0.6.

All four models reproduced the ratio obtained, as shown in Fig. 5. Certain
models predicted a weak $x_{\mathrm{F}}$ dependency of the ratio; however, the
dependency could not be confirmed owing to the large uncertainties in the
results.

Figure 4: Comparison of the measured cross-section with the results of LHCf at
$\sqrt{s}$ = 7 and 13 TeV [5, 6]. The two left panels show the comparison with
the LHCf result at $\sqrt{s}$ = 7 TeV and the pseudorapidity region of
$\eta\,>\,10.94$ (left) and $8.81\,<\,\eta\,<\,8.99$ (middle). The right panel
shows a comparison between $\sqrt{s}$ = 13 TeV and $\eta\,>\,10.94$. The bars
and hatched areas correspond to statistical uncertainties and quadratic
summation of statistics and systematic uncertainties, respectively.

.

Figure 5: Ratios of the RHICf results to the LHCf results [5, 6]. The error
bars represent the uncertainties calculated as a quadratic summation of the
uncertainties in these results. Colored lines indicate MC predictions.

.

##  5 Summary

This study measured the inclusive differential production cross-section of
photons in the very forward region $\eta\,>\,6.1$ with $pp$ collisions at
$\sqrt{s}$ = 510 GeV. Through comparisons of the obtained results with the
LHCf, the Feynman scaling law of forward photon production in a wide energy
range from 510 GeV to 13 TeV was confirmed within the uncertainties. However,
the precision is not sufficient to discuss the weak $x_{\mathrm{F}}$
dependency predicted by certain models, and future studies must focus on
reducing the uncertainties together with LHCf collaboration. Moreover, tests
of the Feynman scaling law with other neutral hadrons $\pi^{0}$and neutrons
are an interesting concept and will be addressed in future publications. They
will contribute to improving the predictions of models with intermediate and
even higher collision energies than LHC.

## Acknowledgements

We thank the staff at BNL, STAR Collaboration, and PHENIX Collaboration for
supporting the experiments. Further, we acknowledge the essential
contributions of STAR members toward the successful operation of RHICf. We are
also grateful to BNL Collider-Accelerator Department for smooth commissioning
and operation of beams, and providing beam parameters used in this analysis
such as luminosity. This work was partly supported by the U.S.–Japan Science
and Technology Cooperation Program in High-Energy Physics, JSPS KAKENHI (No.
JP26247037 and No. JP18H01227), the joint research program of the Institute
for Cosmic Ray Research (ICRR), University of Tokyo, and the National Research
Foundation of Korea (No. 2016R1A2B2008505 and No. 2018R1A5A1025563), and
“UNICT 2020-22 Linea 2” program, University of Catania.

##  Appendix A Cross-section tables

Table 2:  Differential photon production cross-section $\mathrm{d}\,\sigma_{\gamma}/\mathrm{d}\,x_{\mathrm{F}}$ [nb] for each $x_{\mathrm{F}}$bin in the pseudorapidity ranges: $\eta\,>8.5$, $8.0\,<\eta\,<8.5$ and $7.5\,<\eta\,<8.0$. Upper and lower total uncertainties are also reported.  | $\mathrm{d}\,\sigma_{\gamma}/\mathrm{d}\,x_{\mathrm{F}}$ [nb]
---|---
Energy [GeV] | $\eta\,>8.5$ | $8.0\,<\eta\,<8.5$ | $7.5\,<\eta\,<8.0$
0.10-0.15 | $(2.51^{+0.18}_{-0.18})\times 10^{5}$ | $(4.30^{+0.30}_{-0.29})\times 10^{5}$ | $(1.18^{+0.08}_{-0.08})\times 10^{6}$
0.15-0.20 | $(1.93^{+0.13}_{-0.13})\times 10^{5}$ | $(3.29^{+0.21}_{-0.21})\times 10^{5}$ | $(9.00^{+0.56}_{-0.56})\times 10^{5}$
0.20-0.25 | $(1.49^{+0.10}_{-0.11})\times 10^{5}$ | $(2.55^{+0.16}_{-0.19})\times 10^{5}$ | $(6.76^{+0.42}_{-0.44})\times 10^{5}$
0.25-0.30 | $(1.09^{+0.07}_{-0.10})\times 10^{5}$ | $(1.89^{+0.12}_{-0.17})\times 10^{5}$ | $(4.83^{+0.30}_{-0.35})\times 10^{5}$
0.30-0.35 | $(7.52^{+0.47}_{-0.85})\times 10^{4}$ | $(1.28^{+0.08}_{-0.12})\times 10^{5}$ | $(3.51^{+0.22}_{-0.27})\times 10^{5}$
0.35-0.40 | $(5.35^{+0.34}_{-0.59})\times 10^{4}$ | $(8.97^{+0.58}_{-0.74})\times 10^{4}$ | $(2.38^{+0.16}_{-0.17})\times 10^{5}$
0.40-0.45 | $(3.76^{+0.25}_{-0.39})\times 10^{4}$ | $(6.31^{+0.42}_{-0.50})\times 10^{4}$ | $(1.70^{+0.12}_{-0.13})\times 10^{5}$
0.45-0.50 | $(2.78^{+0.19}_{-0.24})\times 10^{4}$ | $(4.61^{+0.31}_{-0.34})\times 10^{4}$ | $(1.17^{+0.09}_{-0.09})\times 10^{5}$
0.50-0.55 | $(1.78^{+0.14}_{-0.18})\times 10^{4}$ | $(3.02^{+0.23}_{-0.24})\times 10^{4}$ | $(7.97^{+0.87}_{-0.91})\times 10^{4}$
0.55-0.60 | $(1.08^{+0.08}_{-0.10})\times 10^{4}$ | $(1.83^{+0.15}_{-0.17})\times 10^{4}$ | $(4.85^{+0.66}_{-0.66})\times 10^{4}$
0.60-0.65 | $(6.02^{+0.64}_{-0.65})\times 10^{3}$ | $(1.04^{+0.10}_{-0.12})\times 10^{4}$ | $(2.65^{+0.36}_{-0.39})\times 10^{4}$
0.65-0.70 | $(3.07^{+0.44}_{-0.44})\times 10^{3}$ | $(5.46^{+0.72}_{-0.92})\times 10^{3}$ | $(1.41^{+0.23}_{-0.27})\times 10^{4}$
0.70-0.75 | $(1.74^{+0.30}_{-0.34})\times 10^{3}$ | $(2.85^{+0.70}_{-0.68})\times 10^{3}$ | $(4.51^{+1.01}_{-1.15})\times 10^{3}$
0.75-0.80 | $(8.07^{+2.35}_{-2.22})\times 10^{2}$ | $(1.17^{+0.44}_{-0.49})\times 10^{3}$ | $(2.10^{+0.73}_{-0.76})\times 10^{3}$
0.80-0.85 | $(2.32^{+1.19}_{-0.88})\times 10^{2}$ | $(2.79^{+1.45}_{-2.01})\times 10^{2}$ | $(5.58^{+3.00}_{-2.91})\times 10^{2}$
0.85-0.90 | $(8.05^{+5.24}_{-5.01})\times 10^{1}$ | $(8.62^{+5.80}_{-8.74})\times 10^{1}$ | $(2.03^{+1.51}_{-1.62})\times 10^{2}$
0.90-0.95 | $(1.01^{+0.92}_{-1.08})\times 10^{1}$ | $(1.76^{+1.70}_{-2.43})\times 10^{1}$ | $(2.91^{+2.94}_{-3.50})\times 10^{1}$
Table 3:  Differential photon production cross-section $\mathrm{d}\,\sigma_{\gamma}/\mathrm{d}\,x_{\mathrm{F}}$ [nb] for each $x_{\mathrm{F}}$bin in the pseudorapidity ranges: $7.0\,<\eta\,<7.5$, $6.5\,<\eta\,<7.0$ and $6.1\,<\eta\,<6.5$.  | $\mathrm{d}\,\sigma_{\gamma}/\mathrm{d}\,x_{\mathrm{F}}$ [nb]
---|---
Energy [GeV] | $7.0\,<\eta\,<7.5$ | $6.5\,<\eta\,<7.0$ | $6.1\,<\eta\,<6.5$
0.10-0.15 | $(3.32^{+0.22}_{-0.21})\times 10^{6}$ | $(8.24^{+0.51}_{-0.51})\times 10^{6}$ | $(1.34^{+0.08}_{-0.08})\times 10^{7}$
0.15-0.20 | $(2.46^{+0.15}_{-0.14})\times 10^{6}$ | $(5.91^{+0.35}_{-0.33})\times 10^{6}$ | $(8.81^{+0.50}_{-0.49})\times 10^{6}$
0.20-0.25 | $(1.82^{+0.11}_{-0.11})\times 10^{6}$ | $(4.20^{+0.26}_{-0.23})\times 10^{6}$ | $(5.69^{+0.35}_{-0.31})\times 10^{6}$
0.25-0.30 | $(1.31^{+0.08}_{-0.08})\times 10^{6}$ | $(2.79^{+0.18}_{-0.16})\times 10^{6}$ | $(3.54^{+0.23}_{-0.20})\times 10^{6}$
0.30-0.35 | $(9.05^{+0.62}_{-0.57})\times 10^{5}$ | $(1.85^{+0.13}_{-0.11})\times 10^{6}$ | $(2.21^{+0.15}_{-0.13})\times 10^{6}$
0.35-0.40 | $(6.17^{+0.45}_{-0.41})\times 10^{5}$ | $(1.22^{+0.09}_{-0.08})\times 10^{6}$ | $(1.34^{+0.09}_{-0.08})\times 10^{6}$
0.40-0.45 | $(4.07^{+0.32}_{-0.29})\times 10^{5}$ | $(8.02^{+0.61}_{-0.55})\times 10^{5}$ | $(8.16^{+0.60}_{-0.54})\times 10^{5}$
0.45-0.50 | $(2.75^{+0.22}_{-0.19})\times 10^{5}$ | $(4.91^{+0.39}_{-0.37})\times 10^{5}$ | $(4.66^{+0.35}_{-0.33})\times 10^{5}$
0.50-0.55 | $(1.78^{+0.18}_{-0.14})\times 10^{5}$ | $(3.02^{+0.25}_{-0.25})\times 10^{5}$ | $(2.58^{+0.21}_{-0.20})\times 10^{5}$
0.55-0.60 | $(1.06^{+0.12}_{-0.10})\times 10^{5}$ | $(1.68^{+0.16}_{-0.17})\times 10^{5}$ | $(1.37^{+0.12}_{-0.12})\times 10^{5}$
0.60-0.65 | $(5.62^{+0.73}_{-0.66})\times 10^{4}$ | $(8.68^{+1.03}_{-1.15})\times 10^{4}$ | $(6.71^{+0.68}_{-0.78})\times 10^{4}$
0.65-0.70 | $(2.87^{+0.50}_{-0.47})\times 10^{4}$ | $(4.28^{+0.62}_{-0.75})\times 10^{4}$ | $(3.18^{+0.38}_{-0.49})\times 10^{4}$
0.70-0.75 | $(1.25^{+0.32}_{-0.31})\times 10^{4}$ | $(1.99^{+0.34}_{-0.45})\times 10^{4}$ | $(1.43^{+0.22}_{-0.30})\times 10^{4}$
0.75-0.80 | $(5.87^{+2.05}_{-2.12})\times 10^{3}$ | $(8.47^{+1.74}_{-2.58})\times 10^{3}$ | $(5.83^{+1.13}_{-1.66})\times 10^{3}$
0.80-0.85 | $(2.98^{+1.40}_{-1.57})\times 10^{3}$ | $(2.68^{+0.67}_{-1.10})\times 10^{3}$ | $(1.65^{+0.40}_{-0.65})\times 10^{3}$
0.85-0.90 | $(1.05^{+0.69}_{-0.77})\times 10^{3}$ | $(7.23^{+3.32}_{-4.89})\times 10^{2}$ | $(4.85^{+2.08}_{-2.83})\times 10^{2}$
0.90-0.95 | $(4.50^{+4.24}_{-4.69})\times 10^{1}$ | $(8.99^{+7.97}_{-10.41})\times 10^{1}$ | $(8.58^{+7.69}_{-8.62})\times 10^{1}$

## References

  * [1] The energy spectrum of cosmic rays beyond the turn-down around $10^{17}$ ev as measured with the surface detector of the pierre auger observatory, The European Physical Journal C 81 (11) (2021) 966.  doi:10.1140/epjc/s10052-021-09700-w.
URL https://doi.org/10.1140/epjc/s10052-021-09700-w

  * [2] A. Aab, et al., Depth of maximum of air-shower profiles at the Pierre Auger Observatory. I. Measurements at energies above 1017.8 eV, Phys. Rev. D 90 (12) (2014) 122005.  arXiv:1409.4809, doi:10.1103/PhysRevD.90.122005.
  * [3] R. U. Abbasi, et al., Depth of Ultra High Energy Cosmic Ray Induced Air Shower Maxima Measured by the Telescope Array Black Rock and Long Ridge FADC Fluorescence Detectors and Surface Array in Hybrid Mode, Astroparticle Physics 858 (2018) 76.  doi:10.3847/1538-4357/aabad7.
  * [4] O. Adriani, et al., Measurement of zero degree inclusive photon energy spectra for $\sqrt{s}$ = 900GeV proton-proton collisions at LHC, Physics Letters B 715 (2012) 298–303.  arXiv:1207.7183, doi:10.1016/j.physletb.2012.07.065.
  * [5] O. Adriani, et al., Measurement of zero degree single photon energy spectra for $\sqrt{s}$ = 7 TeV proton-proton collisions at LHC, Physics Letters B 703 (2011) 128–134.  arXiv:1104.5294, doi:10.1016/j.physletb.2011.07.077.
  * [6] O. Adriani, et al., Measurement of forward photon production cross-section in proton–proton collisions at $\sqrt{s}$ = 13 TeV with the LHCf detector, Phys. Lett. B780 (2018) 233–239.  arXiv:1703.07678, doi:10.1016/j.physletb.2017.12.050.
  * [7] O. Adriani, others, LHCf Collaboration, Measurements of longitudinal and transverse momentum distributions for neutral pions in the forward-rapidity region with the LHCf detector, Phys. Rev. D 94 (3) (2016) 032007.  arXiv:1507.08764, doi:10.1103/PhysRevD.94.032007.
  * [8] O. Adriani, et al., Measurement of very forward neutron energy spectra for 7 TeV proton-proton collisions at the Large Hadron Collider, Physics Letters B 750 (2015) 360–366.  arXiv:1503.03505, doi:10.1016/j.physletb.2015.09.041.
  * [9] O. Adriani, et al., Measurement of energy flow, cross section and average inelasticity of forward neutrons generated in $\mathrm{\sqrt{s}=13TeV}$ proton-proton collisions with the LHCf Arm2 detector, JHEP 07 (2020) 016.  arXiv:2003.02192, doi:10.1007/JHEP07(2020)016.
  * [10] R. P. Feynman, Very high-energy collisions of hadrons, Phys. Rev. Lett. 23 (1969) 1415–1417.  doi:10.1103/PhysRevLett.23.1415.
URL https://link.aps.org/doi/10.1103/PhysRevLett.23.1415

  * [11] O. Adriani, et al., Technical design report of the LHCf experiment (2 2006).
URL https://cds.cern.ch/record/926196

  * [12] Y. Makino, et al., Performance study for the photon measurements of the upgraded LHCf calorimeters with Gd2SiO5 (GSO) scintillators, Journal of Instrumentation 12 (2017) P03023.  doi:10.1088/1748-0221/12/03/P03023.
  * [13] O. Adriani, et al., Performance of RHICf detector during operation in 2017, Journal of Instrumentation 16 (10) (2021) P10027.  arXiv:2108.00634, doi:10.1088/1748-0221/16/10/P10027.
  * [14] C. Adler, A. Denisov, E. Garcia, M. J. Murray, H. Strobele, S. N. White, The RHIC zero degree calorimeter, Nucl. Instrum. Meth. A 470 (2001) 488–499.
  * [15] M. H. Kim, et al., Transverse single-spin asymmetry for very forward neutral pion production in polarized $p+p$ collisions at $\sqrt{s}$ = 510 gev, Phys. Rev. Lett. 124 (2020) 252501.  doi:10.1103/PhysRevLett.124.252501.
  * [16] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, K. Werner, Epos lhc: Test of collective hadronization with data measured at the cern large hadron collider, Phys. Rev. C 92 (2015) 034906.  doi:10.1103/PhysRevC.92.034906.
URL https://link.aps.org/doi/10.1103/PhysRevC.92.034906

  * [17] S. Ostapchenko, Monte carlo treatment of hadronic interactions in enhanced pomeron scheme: Qgsjet-ii model, Phys. Rev. D 83 (2011) 014018.  doi:10.1103/PhysRevD.83.014018.
URL http://link.aps.org/doi/10.1103/PhysRevD.83.014018

  * [18] F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, T. Stanev, Hadronic interaction model sibyll 2.3d and extensive air showers, Phys. Rev. D 102 (2020) 063002.  doi:10.1103/PhysRevD.102.063002.
URL https://link.aps.org/doi/10.1103/PhysRevD.102.063002

  * [19] A. Fedynitch, R. Engel, Revision of the high energy hadronic interaction models phojet/dpmjet-iii, in: Proceedings of the 14th International Conference on Nuclear Reaction Mechanisms, NRM 2015, 15.06 - 19.06., Varenna, Italy, 2015, pp. 301–309, 51.03.03; LK 01.
  * [20] Effective zdc cross section measurement during a dedicated star elastic scattering physics store in run15 (11 2019).  doi:10.2172/1580251.
URL https://www.osti.gov/biblio/1580251

  * [21] R. Ulrich, T. Pierog, C. Baus, Cosmic Ray Monte Carlo Package, CRMC (Feb. 2021).  doi:10.5281/zenodo.4558706.
  * [22] G. Antchev, et al., Luminosity-independent measurements of total, elastic and inelastic cross-sections at $\sqrt{s}=7$ TeV, EPL 101 (2) (2013) 21004.  doi:10.1209/0295-5075/101/21004.
  * [23] G. Antchev, et al., First measurement of elastic, inelastic and total cross-section at $\sqrt{s}=13$ TeV by TOTEM and overview of cross-section data at LHC energies, Eur. Phys. J. C 79 (2017) 103. 16 p.  arXiv:1712.06153, doi:10.1140/epjc/s10052-019-6567-0.
  * [24] J. R. Cudell, V. V. Ezhela, P. Gauron, K. Kang, Y. V. Kuyanov, S. B. Lugovsky, E. Martynov, B. Nicolescu, E. A. Razuvaev, N. P. Tkachenko, Benchmarks for the forward observables at rhic, the tevatron-run ii, and the lhc, Phys. Rev. Lett. 89 (2002) 201801.  doi:10.1103/PhysRevLett.89.201801.
URL https://link.aps.org/doi/10.1103/PhysRevLett.89.201801