\section{Introduction}
The relativistic study of the radial oscillations and of the related instability of a star have begun with the works of Chandrasekhar \cite{1,2}.
The radial oscillations of a star are not related to the emission of gravitational waves, that are intrinsically quadrupolars. The emission of gravitational waves is related to the so called non-radial oscillations.
The theory of the non-radial oscillations of a spherical star started with the work in \cite{3} and improved in \cite{4,5}. In
\cite{6} Chandrasekhar presented a complete relativistic theory of the non-radial oscillations of a black hole as a problem of resonant scattering. There \cite{6}, an incoming gravitational wave perturbates the black hole. The incident wave is in part absorbed and in part reflected by the
event horizon: the quasi-normal frequencies are obtained by imposing no outgoing waves at spatial infinity.
The necessary mathematical formalism to
expand a generic perturbation in a tensorial basis of $10$ spherical harmonics can be found in
\cite{0,01}. This technology has been applied by Chandrasekhar and Ferrari \cite{7,8} and by Ferrari and collaborators (see for example \cite{10,11})
to study the non-radial oscillations of a star.
In particular, the
equations governing the perturbations of the gravitational field ({\it axial}) outside the star reduce \cite{0,6} to a one-dimensional Schrodinger-like wave equation named Regge-Wheeler equation, while the ones governing the perturbations of the matter fluid outside the star
({\it polar}) reduce to the Zerilli \cite{01} equation. The presence of a Schrodinger-like equation for stationary states in the theory of non-radial oscillations is an intriguing fact that can be the starting point of some reasonings concerning a quantum description of a gravitational waves in terms of gravitons.

As pointed in \cite{14}, thanks to the intuitions of the founding fathers of quantum mechanics,
the theory of free electromagnetic waves have played a fundamental role
for the quantization of the free electromagnetic field and for the Planck derivation of the black body radiation.

For the complexity of the
general relativity equations, the same has not
happened for the gravitational waves. Moreover,
an electromagnetic radiation can be confined in thermal equilibrium within a cavity (black body radiation), but the same for gravitons is a complicated task. Confining a gravitational radiation in a box is thus a complicated but necessary task in order to explore fundamental
quantum properties of the gravitational waves. In \cite{15} S. W. Hawking, in order to study the stability of a black hole thermodynamically,
confined a black hole in a fancy box. In \cite{16} it is noticed that a graviton cannot be trapped in a box. Nevertheless, an embedding
within an anti de Sitter spacetime can make the job, although in a idealized way. In \cite{17,18}, L. Smolin concludes that cannot exist a
realistic substance capable to absorb a gravitational radiation, while D. Garfinkle and R.M. Wald \cite{19} presented a counter-example of
\cite{17,18} criticized in \cite{20}. Finally, in \cite{21} T. Padmanabhan and T.P. Singh presented
an interesting example to confine a linearized gravitational radiation
in thermal equilibrium within a cavity. The first goal of this paper is to obtain a physically motivated formula for the spectrum
of gravitons trapped in a box, after that
we study the statistical mechanics of gravitons in a box ('hard wall')
and some possible applications to the confining gravitons procedure. In particular, we are able to explain the
semiclassical Bekenstein-Hawking
black hole formula $S_{BH}=K_B A_h/L_P^2$ (with $K_B$ the Boltzmann constant, $L_P$ the Planck length and $A_h$ the proper area of
the event horizon) for a Schwarzschild black hole
starting from the discrete spectrum of gravitons confined within the black hole.

In section 2 we rewrite the Regge-Wheeler and Zerilli equations in the special case of a perturbed Minkowskian spacetime.
In section 3 we obtain a discrete spectrum for gravitons in a spherical box,
while section 4 is devoted to the statistical description of 'trapped' gravitons. In section 5 we analyze the thermodynamics of trapped
gravitons. In section 6 we apply the machinery of sections above to obtain the semi-classical black hole entropy formula
in a simple and clear way, while in section 7 the case with summation over the legendre index $m$ is done.
Finally, section 8 is devoted to some conclusions and final remarks.

\section{The equations in the vacuum}

The first goal of this paper is to obtain a physically motivated expression for the spectrum of gravitons trapped in a
spherical box. It is customary (see \cite{E1} and references therein) in perturbative quantum field theory,
to depict graviton field $h_{\mu\nu}$ as a perturbation about the flat metric $diag(1,-1,-1,-1)$.
To this purpose, it is sufficient to consider a gravitational wave traveling in a Minkowski
(perturbed) spacetime. Hence,
we consider perturbations of the Minkowski spacetime (we initially use geometrized unity with $G=c=1$)
\begin{equation}
ds^2=dt^2-r^2\sin^2\theta\;d{\phi}^2-dr^2-r^2 d{\theta}^2.
\label{1}
\end{equation}
As usual, we indicate with $g_{ik}^{(0)}$ the unperturbed Minkowskian metric given by (\ref{1}) and with $h_{ik}$ a small perturbation with
$|h_{ik}|<<|g_{ik}^{(0)}|$. The perturbed metric is
\begin{equation}
g_{ik}=g_{ik}^{(0)}+h_{ik}.
\label{2}
\end{equation}
Since of the spherical symmetry of the 'bare' unperturbed metric $g_{ik}^{(0)}$, the perturbating term $h_{ik}$ can be expanded in a basis of spherical tensorial harmonics \cite{0,01} depending on the real spherical harmonics
$Y_{\ell m}(\theta,\phi), \ell\in\mathbb{N}
, m\in\mathbb{Z}
, m\in [-\ell, +\ell]$. In this frame, the polar perturbations
$h_{ik}^{p}$ are the ones dependent on harmonics that under parity operator look like ${(-1)}^{\ell}$, while the axial perturbations
$h_{ik}^{a}$ look like ${(-1)}^{\ell+1}$.\\
In the diagonal gauge \cite{10} the perturbed metric becomes, for the axial case:
\begin{equation}
h_{ik}^a=
\begin{pmatrix}
(t) & (\phi) & (r) & (\theta)\\
0 & h_0(t,r)\sin\theta\;Y_{\ell m,\theta} & 0 & -h_0(t,r)\frac{1}{\sin\theta}Y_{\ell m,\phi}\\
h_0(t,r)\sin\theta\;Y_{\ell m,\theta} & 0 & h_1(t,r)\sin\theta\;Y_{\ell m,\theta} &
0\\
0 & h_1(t,r)\sin\theta\;Y_{\ell m,\theta} & 0 & -h_1(t,r)\frac{1}{\sin\theta}\;Y_{\ell m,\phi}\\
-h_0(t,r)\frac{1}{\sin\theta}Y_{\ell m,\phi} & 0 & -h_1(t,r)\frac{1}{\sin\theta}\;Y_{\ell m,\phi} &
0
\end{pmatrix},
\label{7}
\end{equation}
and for the polar case:
\begin{equation}
h_{ik}^p=
\begin{pmatrix}
(t) & (\phi) & (r) & (\theta)\\
2N(t,r) Y_{\ell m} & 0 & 0 & 0\\
0 & -2r^2\sin^2\theta\;H_{11}(t,r,\theta,\phi) & 0 & -r^2 V(t,r) X_{\ell m}\\
0 & 0 & -2L(t,r) Y_{\ell m} & 0\\
0 & -r^2 V(t,r) X_{\ell m} & 0 & -2 r^2 H_{33}(t,r,\theta,\phi)
\end{pmatrix},
\label{8}
\end{equation}
where
\begin{eqnarray}
X_{\ell m}(\theta,\phi) &=& 2Y_{\ell m,\theta,\phi}-2\cot\theta\;Y_{\ell m,\phi}\label{8b}\\
W_{\ell m}(\theta,\phi) &=& Y_{\ell m,\theta,\theta}-\cot\theta\;Y_{\ell m,\theta}-
\frac{1}{\sin^2\theta}\;Y_{\ell m,\phi,\phi}\nonumber\\
H_{11}(t,r,\theta,\phi) &=&TY_{\ell m}+\frac{V}{\sin^2\theta}\;Y_{\ell m,\phi,\phi}+
V\cot\theta\;Y_{\ell m,\theta}\nonumber\\
H_{33}(t,r,\theta,\phi) &=&TY_{\ell m}+VY_{\ell m,\theta,\theta}.
\end{eqnarray}
As usual, the static nature of the unperturbed metric (\ref{1}) allows to assume that the perturbation functions all have the time dependence
$e^{-\imath\omega t}$ \cite{6} (with $\omega$ a constant): this is equivalent to a Fourier analysis of the metric coefficient with
frequency $\omega$. Since we are considering vacuum perturbations, we have
$\delta T_{ik}=0$.

 \subsection{Axial perturbations with $\ell\geq 2$}

As well known \cite{0,6}, after defining the Regge-Wheeler function
$Z_{\ell m}^{(a)}$ by $h_{1\ell m}=rZ_{\ell m}^{(a)}$, the axial perturbations are driven by
the Regge-Wheeler\footnote{We named the (\ref{23}) Regge-Wheeler equation for simplicity, although it refers to
	the Schwarzschild case.} equation \cite{14}:
\begin{equation}
Z_{\ell m,r,r}^{(a)}+{\omega}^2 Z_{\ell m}^{(a)}=\frac{\ell(\ell+1)}{r^2}Z_{\ell m}^{(a)}.
\label{23}
\end{equation}
Equation (\ref{23}) reminds the Schrodinger equation for stationary states with vanishing potential $V(r)$.

\subsection{Polar perturbations with $\ell\geq 2$}
By following the same manipulations of \cite{14}
and after defining $Z_{\ell m}^{(p)}=rL_{\ell m}$, the
solutions of the linearized Einstein's equations become:
\begin{eqnarray}
& & qV_{\ell m}(r,\omega)=-L_{\ell m}(r,\omega)-\frac{1}{r}\int_0^r L_{\ell m}(r^{\prime},\omega)dr^{\prime},\label{f1}\\
& & N_{\ell m}(r,\omega)=L_{\ell m}(r,\omega)+2r\int_0^r\frac{L_{\ell m}(r^{\prime},\omega)}{{r^{\prime}}^2}dr^{\prime},
\label{f2}\\
& & Z_{\ell m,r,r}^{(p)}-\frac{\ell(\ell+1)}{r^2}Z_{\ell m}^{(p)}+{\omega}^2 Z_{\ell m}^{(p)}=0.\label{f3}
\end{eqnarray}
To obtain the equation (\ref{f1}) and (\ref{f2}) we have imposed the regularity of $L_{\ell m}$ for $r\rightarrow 0$.

\section{Gravitons in a box}

Starting from this section, we restore the constants $c$ and $G$.
In the literature of statistical black hole entropy, we can often see gravitons inside the horizon depicted as an
ensemble of $n$ non-interacting massless particles
with (continuum) angular frequency \cite{E2} $\omega$ or by considering discrete ad hoc expressions for $\omega$ like
$\omega\sim n$ or inspired \cite{E3,CC1,CC2,CC3}
by quasi-normal black hole frequencies. However, although the main proposals present in the literature are physically possible, a derivation from the onset
\footnote{A convincing quantum gravity theory is still not at our disposal.}
of the (discrete) expected frequency of trapped gravitons is still lacking.
Hence, our primary goal in this section is to
obtain a physically reasonable ('phenomenological') formula for the frequency of trapped gravitons motivated by the techniques of perturbative quantum field
theory (see \cite{E1} and references therein), i.e. perturbation around a Minkowski spacetime.

To start with,
as well known, \cite{0,01,6,14} equations (\ref{23}) and (\ref{f3}) governing the axial and polar perturbations respectively look like the radial part of the Schrodinger equation for free particles and for stationary states. The only difference is that to the metric ('wave') function
$L_{\ell m}$ are associated, as the quantum mechanical case, the Legendre polynomials $Y_{\ell m}(\theta,\phi)$, while for the axial case to
$Z_{\ell m}^{(a)}$ are associated the Legendre polynomials given by $Y_{\ell m,\theta}$ and $Y_{\ell m,\phi}$ (see \cite{7}). This means that the analogy with the Schrodinger equation is more stringent with the polar case. Hence, in the following we only consider equation (\ref{f3}),
although the following considerations concerning the (\ref{f3}) obviously still hold for the (\ref{23}).

As well known, a particle with rest mass $m$,
potential $V(r)$ and energy $E$ satisfies the spherical Schrodinger equation with a wave
function ${\psi}_E(r,\theta,\phi)$ that can be split as
${\psi}_E(r,\theta,\phi)=R(r)Y_{\ell m}(\theta,\phi)$. The angular part $Y_{\ell m}$ satisfies the eigenvalues equation
\begin{equation}
-\left[\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}
\left(\sin\theta\frac{\partial}{\partial\theta}\right)+\frac{1}{{\sin}^2\theta}
\frac{{\partial}^2}{\partial{\phi}^2}\right]Y_{\ell m}=\ell(\ell+1)Y_{\ell m},
\label{42}
\end{equation}
while for the radial part, after defining $u_E(r)=rR(r)$, we have
\begin{equation}
-\frac{{\hbar}^2}{2m}u_{E,r,r}+\left[V(r)+\frac{\ell(\ell+1){\hbar}^2}{2m r^2}\right]u_E=E u_E,\;\;u_E(r=0)=0.
\label{43}
\end{equation}
Consider now the free case by setting $V(r)=0$. The (\ref{43}) becomes
\begin{equation}
-{\hbar}^2 u_{E,r,r}+\frac{\ell(\ell+1){\hbar}^2}{r^2} u_E=P^2 u_E.
\label{kk}
\end{equation}
In the (\ref{43}) the Planck constant arises thanks to the usual form of the impulse
operator $\mathbf{P}$ as $\mathbf{P}=-\imath\hbar\nabla$. To make the analogy between the (\ref{43}) and the (\ref{f3})
(and (\ref{23})) we can write the (\ref{43}) as
\begin{equation}
-u_{E,r,r}+\frac{\ell(\ell+1)}{r^2}u_E=k^2u_E,\;\;\;k=\sqrt{\frac{2mE}{{\hbar}^2}}.
\label{44}
\end{equation}
Since for the (\ref{44}) $E=\frac{P^2}{2m}$, we have $k=\frac{|P|}{\hbar}$.
We can now rewrite equation (\ref{f3}) as
\begin{equation}
-Z_{\ell m,r,r}^{(p)}+\frac{\ell(\ell+1)}{r^2}Z_{\ell m}^{(p)}=\frac{{\omega}^2}{c^2}Z_{\ell m}^{(p)}.
\label{45}
\end{equation}
Thanks to the (\ref{45}) and (\ref{44}) we can write $|P|=\frac{\hbar\omega}{c}$. This means that the equation governing the polar
\footnote{For the axial perturbation the analogy is only with the radial eqution (\ref{23}).}
perturbations can be interpreted as a wave function for free massless gravitons with energy $E=\hbar\omega$ and momentum
$|P|=\frac{E}{c}$.
As a consequence, the way we have to take evident the quantum nature of gravitons is to confine gravitons in a finite
(spherical in our context) box. As preliminarily discussed in the introduction, this is not a simple task. Various attempts can be found in the
literature (see for example \cite{15,16,17,18,19,20,21}) with opposite answers. In this paper we use the pragmatic view present in
\cite{15} avoiding the issues of a practical realization. Quite
reasonably \cite{21}, we could think to a gravitational radiation in thermal equilibrium with the (hard) wall of a spherical box at a given temperature $T$. Hence, we have a potential $V(r)$ that is zero inside the box and practically infinity at the boundary
$R$.

To start with, we must specify the solution of equation (\ref{45}). As well known \cite{14}, the regular solution for
$L_{\ell m}(r,\omega)$ can be obtained in terms of regular Bessel functions $j_{\ell}(kr)$ with $k=\frac{\omega}{c}$.
Thus we have $L_{\ell m}(r,\omega)=A_{\ell}\;j_{\ell}(kr),A_{\ell}\in\mathbb{R}$. The regular Bessel functions for $r\rightarrow 0$ look like
$j_{\ell}(kr)\sim {{(kr)}^{\ell}}/(2\ell+1)!!$, while for large values of $kr$ we have
\begin{equation}
j_{\ell}(kr)\simeq\frac{1}{kr}\cos\left[kr-\frac{(\ell+1)\pi}{2}\right].
\label{51}
\end{equation}
It should be noticed that the (\ref{51}) is a good approximation also for $kr > \ell$.
To confine our gravitational radiation within a spherical box of radial radius $R$, we must impose the boundary condition
$L_{\ell m}(R,\omega)=0$: this condition is a Dirichlet boundary condition and in fact the confining box
is a Dirichlet boundary condition for the metric function $Z_{\ell m}$.

As a consequence, with the help of (\ref{51}), the following discrete spectrum for the so trapped gravitons
does arise:
\begin{equation}
{\omega}_{\ell n}\simeq\frac{c}{2R}\left(2+\ell+2n\right)\pi,\,\;\ell\geq 2,\;\;n\in\mathbb{N}.
\label{52}
\end{equation}
As expected, finite size effects lead to the manifestation of a quantum behavior for gravitons. The classical behavior is restored
for $R\rightarrow\infty$ where a strictly continuum spectrum arises.
Note that the allowed frequencies look like the ones of an harmonic oscillator but with the further
dependence on the Legendre index $\ell$.
As the usual harmonic oscillator, the ground state of (\ref{52}) has a non-vanishing energy $E_{0\ell}$ given by
$E_{\ell 0}=\hbar{\omega}_{\ell 0}=\hbar\frac{c(\ell+2)\pi}{2R}$.
In the following, we will discuss how to interpret equation (\ref{52}).
Thanks to the similarity depicted above, we has been able to obtain a physically reasonable expression for the spectrum of trapped gravitons.
By forcing this analogy,
we are tempted to promote the indices $\{\ell,m,n\}$ to quantum numbers. However, note that
formula (\ref{52}) is independent from the azimuthal index $m$ and as a result the energy levels are degenerate with a degeneration factor
given by $2\ell+1$. However, it should be noted that the Regge-Wheeler and Zerilli equations (\ref{f3}) and (\ref{23}) have been obtained
by perturbing a spherically symmetric vacuum region. As a consequence, for any value of the Legendre index $\ell$ we can have
quadrupolar ($\ell=2$), sextupolar ($\ell=3$) perturbations. To any kind of perturbations we can associate a different species of gravitons, i.e.
we can have a given harmonic oscillator spectrum for any given fixed $\ell$: thus we can have
quadrupolar, sextupolar... gravitons. Hence, only the index $n$ could be promoted to a quantum index, whereas the index $\ell$ can be seen as
a species index. This interpretation of the index $\ell$ could be of interest, for example, in relation to the so named 'species problem'
\cite{s1}. In fact, 'ordinary' entropy depends on the number of species, while the black hole entropy is species-independent.
In our case, we have infinity 'species' of gravitons
\footnote{Only massless particles can stably survive inside a black hole.}
that must be summed up in order to obtain the geometric expression of the black hole entropy.
In the following, we adopt this point of view. Nevertheless, for completeness of presentation, in section 7 we study the
interpretation with all
the indices $\{\ell,m,n\}$ promoted to quantum numbers.

\section{Statistical mechanics of gravitons in a box}

Our starting point \cite{hua} is the canonical partition function $Z_{\ell}$ for the discrete spectrum (\ref{52}) with the summation over the quantum index $n$:
\begin{equation}
Z_{\ell}=\sum_{n=0}^{\infty} e^{-\beta\hbar{\omega}_{\ell n}},\;\;\;\beta=\frac{1}{K_B T}.
\label{55}
\end{equation}
First of all, we study the case  of the partition function (\ref{55}) with a given Legendre index $\ell$ .
With a trivial algebra we obtain
\begin{equation}
Z_{\ell}=\frac{e^{-\left[\frac{c\beta}{2R}(\ell+2)\pi\hbar\right]}}
{1-e^{-\left[\frac{c\beta\pi\hbar}{R}\right]}}.
\label{56}
\end{equation}
The Helmholtz free energy $F_{\ell}$ is $F_{\ell}=-K_B T\ln(Z_{\ell})$.
Hence, the numerator in the (\ref{56}) is a constant independent from
$\beta$ and for the entropy $S_{\ell}$ we obtain the one of a discrete harmonic oscillator with
angular frequency $\omega=c\pi/R$. For the mean energy $<E_{\ell}>$ we obtain
\begin{equation}
<E_{\ell}>=\frac{c}{2R}(\ell+2)\pi\hbar+
\frac{c\pi\hbar\;e^{-\left[\frac{c\beta\pi\hbar}{R}\right]}}{R\left({1-e^{-\left[\frac{c\beta\pi\hbar}{R}\right]}}\right)}.
\label{57}
\end{equation}
By denoting	with ${\omega}_g=c\pi/R$, expression (\ref{57}) becomes
\begin{equation}
<E_{\ell}>=\left(<n>+\frac{(\ell+2)}{2}\right)\hbar{\omega}_g.
\label{58}
\end{equation}
It is intriguing to note that, after formally setting $\ell=1$ in (\ref{58}), we obtain the usual formula for an ordinary three dimensional
harmonic oscillator\footnote{This is what we expect, thanks to its dipolar nature, for the electromagnetic radiation.}.
The term $\frac{(\ell+2)\hbar{\omega}_g}{2}$ in (\ref{58}) is nothing else but the non-vanishing energy of the ground state.
In the Planck-like derivation, the ground state energy is
usually set to zero: in our case this implies that the term dependent on $\ell$ is absent. In fact, this term
\footnote{This represents the zero point energy.}
is typically subtracted to $<E>$. In this way, the usual Planck distribution emerges with the mean number occupation $<n_{\omega}>$ for the level
with energy $\hbar\omega$ given by
\begin{equation}
<n_{{\omega}_g}>=\frac{1}{e^{\beta\hbar{\omega}_g}-1},
\label{59}
\end{equation}
that is the Planck distribution. It may be argued that if we consider a general gravitational perturbation, then a summation over
$\ell$ and $m$
must be performed. As stated at the end od section 3, in the following we consider the index $\ell$ as representing the species number and we leave
the case of the summation over $m$ at the appendix.
We must thus calculate the full partition function after marginalization of $\ell$. We obtain:
\begin{equation}
Z_g=\sum_{\ell=2}^{\infty}\sum_{n=0}^{\infty} e^{-\beta\hbar{\omega}_{\ell n}}=
\frac{e^{-\left(\frac{2c\pi\beta\hbar}{R}\right)}}
	{\left[1-e^{-\left(\frac{c\pi\beta\hbar}{2R}\right)}\right]}
	\frac{1}{\left[1-e^{-\left(\frac{c\pi\beta\hbar}{R}\right)}\right]}.
\label{60}
\end{equation}
As a consequence of the (\ref{60}), after calculating the Helmholtz energy $F_g=-K_B T\ln(Z_g)$, we see that the entropy is the one of two
harmonic oscillators with angular frequencies ${\omega}_1=c\pi/(2R)$ and ${\omega}_2=2{\omega}_1$. For the mean energy
$<E_g>$ we obtain in the space of frequencies ${\nu}$
\begin{equation}
<E_g>=h{\nu}_0+\frac{h{\nu}_1}{\left[e^{\beta h{\nu}_1}-1\right]}+\frac{h{\nu}_2}{\left[e^{\beta h{\nu}_2}-1\right]},\;
{\nu}_0=\frac{c}{R},\;{\nu}_1=\frac{c}{4R},\;{\nu}_2=\frac{c}{2R}.
\label{61}
\end{equation}
Note that the ground state energy $h{\nu}_0$ can be obtained from the one given in (\ref{57}), i.e.
$\frac{c}{2R}(\ell+2)\pi\hbar$, by setting $\ell=2$: this could be a manifestation of the quadrupolar nature of the gravitational radiation.
After subtracting this term we obtain
\begin{equation}
<E_g>=\frac{h{\nu}_1}{\left[e^{\beta h{\nu}_1}-1\right]}+\frac{h{\nu}_2}{\left[e^{\beta h{\nu}_2}-1\right]}.
\label{62}
\end{equation}
The (\ref{62}) is a summation of two discrete Planck distributions with proportional frequencies, i.e.
${\nu}_2=2{\nu}_1$. In this regard, gravitons seem to differentiate from Planck distributions of photons in thermal
equilibrium within a cavity. The term with ${\nu}_1$ is a consequence of the summation over the Legendre index $\ell$ and thus
can be seen as the mean contribution to $E_g$ related to the superposition of all possible kind of gravitational perturbations
(quadrupolar, sextupolar...). More generally,
a gravitational radiation in thermal equilibrium within a cavity can be seen as the superposition of $N$
harmonic oscillators with fundamental
frequencies ${\nu}_1$ and ${\nu}_2=2{\nu}_1$ with partition function $Z_T$ \cite{hua}
\begin{equation}
Z_T=Z_g^N.
\label{63}
\end{equation}
All these facts imply a difference with the ordinary electromagnetic radiation in a cavity where the term with
${\nu}_1$ is absent. To this purpose, it should be always held fixed in mind the fact that electromagnetic radiation is typically dipolar, while
the gravitational one is quadrupolar.

\section{Thermodynamics of gravitons}

First of all, thanks to the (\ref{63}), for a system of $N$ oscillators we have:
\begin{equation}
Z_T=
{\left(\frac{e^{-\left(\frac{2c\pi\beta\hbar}{R}\right)}}
{\left[1-e^{-\left(\frac{c\pi\beta\hbar}{2R}\right)}\right]}
\frac{1}{\left[1-e^{-\left(\frac{c\pi\beta\hbar}{R}\right)}\right]}\right)}^N.
\label{64}
\end{equation}
The internal energy $<E_T>=U$ can be derived from
(\ref{64}) in the usual way as $-{(\ln Z_T)},{\beta}$:
\begin{equation}
U=\frac{c\pi\hbar N}{2R\left[e^{\frac{\beta c\pi\hbar}{2R}}-1\right]}+\frac{c\pi\hbar N}{R\left[e^{\frac{\beta c\pi\hbar}{R}}-1\right]}.
\label{65}
\end{equation}
Note that in the limit for $R\rightarrow\infty$ (at fixed $T$ or more precisely for $\hbar\rightarrow 0$) we have $U\rightarrow 2N K_B T$.
By inspection of the (\ref{65}) we see that, in our discretized system, the internal energy $U$ is a function of the thermodynamic
temperature $T$ and the radius $R$. Hence, we must pay care to use Maxwell's relations involving partial derivatives with respect to the volume
$V=4/3\pi R^3$. As an example, by the first law we have $T dS=dU+PdV$. Since $U=U(R,T)$, it is not more true that, for example
$\frac{\partial S}{\partial V}=\frac{P}{T}$, with $P$ the pressure.\\
To start with, we thus calculate the entropy by using the usual definition
\begin{equation}
S=-K_B\sum_{\ell=2}^{\infty}\sum_{n=0}^{\infty}P_{n\ell}\ln(P_{\ell,n}),\;\;
P_{n\ell}=\frac{e^{-\beta E_{\ell n}}}{Z_T}.
\label{66}
\end{equation}
After a simple algebra we obtain
\begin{equation}
TS=NK_B T\ln(Z_g)+U,
\label{67}
\end{equation}
and with $F=U-TS$ we regain the usual relation $F_T=-NK_B T\ln(Z_g)$. As a consequence, since $dF_T=-SdT-PdV$ certainly we have
that $\frac{\partial F_T}{\partial T}=-S$ and $\frac{\partial F_T}{\partial V}=-P$.
Hence, for $S$  we obtain
\begin{eqnarray}
& & S=-N K_B\left[\ln\left(1-e^{-\frac{X}{2}}\right)+\ln\left(1-e^{-X}\right)\right]+\label{68}\\
& & +\frac{c\pi\hbar N e^{-\frac{X}{2}}\left[1+3 e^{-\frac{X}{2}}\right]}{2TR\left(1-e^{-X}\right)},\;\;
X=\frac{c\pi\beta\hbar}{R},\nonumber
\end{eqnarray}
and for the pressure $P$
\begin{equation}
P=\frac{c\pi N\hbar\left[e^{-\frac{X}{2}}+3e^{-X}\right]}{2R^2\left(1-e^{-X}\right)}\frac{dR}{dV},
\label{69}
\end{equation}
with $dR/dV=3^{-1}{(3/(4\pi))}^{1/3} V^{-2/3}$. It is a simple matter to verify that, thanks to the formula (\ref{65}), we regain the usual relation for a
radiation field $PV=U/3$. In the limit for $\hbar\rightarrow 0$ (classical limit) we obtain in particular
$PV=2NK_B T/3=U/3$.\\
Summarizing, in this section we have shown that the application of the (\ref{52}) for trapped gravitons
in the semiclassical limit, by considering the Legendre index as
labeling the species of gravitons, leads to well defined thermodynamic quantities together with the expected equation of state
(\ref{69}) suitable for a radiation field. In the appendix we performed the same calculations of this section but with the degeneracy factor
$(2\ell+1)$ taken into account.\\
In the next section we apply the results of this section in order to obtain the Schwarzschild
black hole entropy from a statistical point of view.

\section{An application: the Schwarzschild black hole entropy}

An important issue of modern physics is the statistical understanding of the black hole entropy $S_{BH}=K_B A_h/L_P^2$, i.e. the calculation
of the internal degrees of freedom of a black hole. A lot of attempts are present in the literature (see for example
\cite{22,23,24,25,26,27,27b,27bb} and references therein). Typically,
in many approaches black hole entropy arises from a proposed quantum gravity theory
(string theory, loop quantum gravity) and then the validity of these derivations is related to the unproven used theory.
Other approaches use excitation modes propagating near the horizon \cite{22,23,24}. In particular, in \cite{22} and
\cite{24} the authors suppose that the black hole entropy is the entanglement between internal and external (with respect to the horizon)
degrees of freedom, but the expression so obtained is ultraviolet divergent.
Generally, the main obstacle to get a statistical derivation of the Bekenstein-Hawking entropy is due to the presence of the Planck length
$L_P$, i.e. how to obtain the constant $L_P$ without using a prosed quantum gravity theory. In fact, in many attempts the proportionality
$S\sim A$ is obtained but the constant of proportionality is introduced by hand.\\
In the following, we show that, thanks to the formula (\ref{52}) and to the fact that we do not perform the continuum (classical) limit
by retaining $\hbar\neq 0$, the Bekenstein-Hawking entropy formula can be obtained in a clear and simple way. In practice, since the formula $S_{BH}=K_B A_h/L_P^2$ is semi-classical, we use
a discrete version of the canonical ensemble distribution.\\
We note that black hole is a natural arena to apply the results of sections 4 and 5 since nothing can escape from a black hole and
the event horizon can be considered as a perfect hard wall. Inspired by perturbative quantum field theory \cite{E1},
we obtained the frequency formula
(\ref{52}) by perturbating a spherical region in a Minkowski spacetime,
while the Schwarzschild geometry is not flat.
Nevertheless, many arguments can be invoked to justify the application of this section.

In the situation considered in the sections above gravitons are trapped in a spherical region and the (\ref{52}) it gives the spectrum of these
gravitons. We may suppose the interior of a black hole as composed of a gas of massless excitations provided by gravitons. These gravitons
with their frequencies (energy) spectrum (\ref{52}) generate the curvature inside the event horizon at $r=R_s$.
The condition that these so obtained gravitons generate a black hole is provided by equating the internal energy of this gas with the ADM
mass-energy $Mc^2$ of the black hole\footnote{This can be obtained by putting a sufficiently large amount $N$ of
oscillators within a given radius $R$, as shown by equation (\ref{71})}. This is a standard 'empirical' procedure that avoids
the shortcomings due to the lack of a complete quantum gravity theory.\\
Moreover, there exists evidence that
(see for example \cite{22,23,24}) excitation modes inside the Schwarzschild black hole propagates near the event horizon.
This fact maybe a consequence of the holographic principle  \cite{28,29}
since the internal degrees of freedom are in some sense 'projected'
onto the black hole event horizon, i.e. the holographic screen. It is a known fact that the metric near the horizon is 'almost'
flat (see for example \cite{a2}). Near the horizon of a Schwarzschild black hole, the metric looks like the one of
a flat Minkowski spacetime times a sphere of radius $R_s$ given by the one of the event horizon, i.e.
$R_s=2GM/c^2$ (Rindler spacetime).
An experimenter very close to the event horizon cannot distinguish between a Rindler coordinate system
in a flat spacetime with the black hole Schwarzschild spacetime
\footnote{The Hawking temperature is in fact defined by the surface gravity on the event horizon.}
and he is legitimate to adopt usual thermodynamics.
This point of view also inspired the works based on the concept of entanglement entropy
(see again the nice paper \cite{E1}).\\
Moreover, In the study of the thermodynamics of (at least) static black holes, for its volume $V_h$ \cite{a7,33} it is appropriate
(named in this case thermodynamic volume) the usual Minkowskian form $V_h=4\pi R_s^3/3$, irrespective of the fact that inside the
horizon the spatial coordinate $r$ becomes timelike and the time coordinate $t$ spacelike and the usual definition of a volume is rather problematic.
Summarizing, also the holographic principle conspires to elevate the formula (\ref{52}) as a physically reasonable expression
for the spectrum of trapped gravitons inside the event horizon with the black hole degrees of freedom stored
near the horizon. In the following, we show that this is encoded by the behavior of $N$ that looks like $A\sim R_s^2$ rather than the
volume $V$ of the event horizon.

To start with, we must specify the internal energy $U$. For a Schwarzschild black hole, as stated above,
this is done by its rest ADM energy, i.e.
$U=Mc^2$, where $M$ is the ADM black hole mass. By the usual relation $2GM/c^2=R$ we obtain $U=c^4 R/(2G)$. The first step is the determination
of the behavior of $N$ and $T_i$. This can be done by considering a suitable thermodynamic limit.
At a first look, we can use the standard set up $N\sim R^3$. With this choice
and at a classical level with with the temperature independent from the Planck constant,
formula (\ref{65}) implies that\footnote{This can be understood by the fact that
in the limit $\hbar=0$ we have $U=2NK_BT$.} $T_i\sim 1/R^2$. When the expression $T_i\sim 1/R^2$ is inserted in the entropy formula
(\ref{68}), in the limit $R\rightarrow\infty$ we obtain $S\rightarrow 0$. This result is obviously
not acceptable. Moreover, the constraint $U\sim R,\;\forall R>0$ rule out this possibility.
The only non-classical possibility leading to $S\sim A$
is that, according to the holographic principle \cite{28,29}, we have
$N\sim R^2$. Since for a black hole $U\sim R$, thanks to the (\ref{65}), for the internal temperature $T_i$ of a black hole we must have
\begin{equation}
T_i=\alpha T_h=\frac{\alpha c\hbar}{4\pi K_B R},\;\;\alpha\in (0,\infty).
\label{70}
\end{equation}
Formula (\ref{70}) means that the internal temperature is proportional to the one of the event horizon. There is not a priori reason to assume
$\alpha=1$. In fact, the temperature $T_h$ is proportional to the surface gravity parameter on the event horizon and is normalized in such a way that the norm of the timelike killing vector is the unity at spatial infinity. However, for an unlucky observer inside the black hole the outside
spatial infinity is not attainable and the usual normalization is not available. Hence, the result that $\alpha\neq 1$ is
physically reasonable, although we expect $\alpha$ of the order of unity.\\
As a consequence, thanks to the (\ref{70}) and using $U=c^4 R/(2G)$, the (\ref{65}) it gives
\begin{equation}
R=\sqrt{\pi N}L_P\sqrt{\frac{1}{\left[e^{\frac{{2\pi}^2}{\alpha}}-1\right]}+
	\frac{2}{\left[e^{\frac{4{\pi}^2}{\alpha}}-1\right]}}.
\label{71}
\end{equation}
Note the presence of $L_P$ in the (\ref{71}) that is a consequence of the presence of $\hbar$ in (\ref{65}).
Finally, with
(\ref{71}) the entropy (\ref{68}) becomes
\begin{eqnarray}
& & S=K_B Y(\alpha)\frac{A_h}{4L_P^2},\;\;A_h=4\pi R^2,\label{72}\\
& & Y(\alpha)=\frac{b}{\alpha\pi^2\left(3+e^{\frac{2\pi^2}{\alpha}}\right)}\nonumber\\
& & b=-\alpha e^{\frac{4\pi^2}{\alpha}}\ln\left(1-e^{-\frac{2\pi^2}{\alpha}}\right)-
\alpha e^{\frac{4\pi^2}{\alpha}}\ln\left(1-e^{-\frac{4\pi^2}{\alpha}}\right)+6\pi^2+\nonumber\\
& & +2\pi^2 e^{\frac{2\pi^2}{\alpha}}+\alpha \ln\left(1-e^{-\frac{2\pi^2}{\alpha}}\right)+
 \alpha\ln\left(1-e^{-\frac{4\pi^2}{\alpha}}\right).
\end{eqnarray}
It is a simple matter to see that the equation $Y(\alpha)=1$ has (remarkable)
a unique solution approximatively given by
$\alpha\simeq 2.225$. First of all, it should be noted that formula (\ref{52}) is an approximation, but a good approximation for the roots of
$j_{\ell}$. In particular, the maximum error is for the first zero that is of the order of the $3\%$. This error quickly decreases for the
remaining zeros of $j_{\ell}$ and as a result
a small correction to the value obtained for $\alpha$ arises with its 'true' value of the order of
$2$. Physically, this factor can be interpreted as due to the fact that the temperature $T_h$ is measured by an observer at spatial infinity, while
for an observer within the event horizon the outside spatial infinity is not available. Moreover, note that, thanks to
the (\ref{69}), a non-vanishing pressure is present inside the horizon due to the gravitons radiation pressure. As a consequence, since
$T_idS_{BH}=dU+PdV_h$ and $S$ and $U$ are left unchanged\footnote{Remember that for the first law
	at $R=2GM/c^2$ we have $T_hdS_{BH}=dU$}, we must have $T_i\neq T_h$. This pressure $P$ looks like $P\sim 1/R^2$, that is a further
indication that gravitons behave holographic inside the black hole. This is an interesting result that solves the problem concerning the presence
of a work term in the usual black hole first law of thermodynamics. An observer measures the Hawking flux at spatial infinity with a temperature
$T_h$ and thus is forced to associate to the black hole an entropy $S_{BH}\sim A_h/4$, a temperature $T_h$
(the temperature he/her perceives) and an internal energy given by the ADM mass.
For an internal observer the situation is different. He perceives the temperature of the gravitons gas present within the horizon
and also the pressure generated by this gas. Both agree that the entropy is $S_{BH}\sim A_h/4$, but the 'internal' point of view in a more usual way since he can see gravitons inside the black hole. Hence the usual form of the first law $T_h dS_{BH}=c^2 dM$ is a consequence of the
lack of knowledge of the nature of the degrees of freedom inside the event horizon of the external observer.

It should be noticed that our result has been obtained only by requiring that the internal energy is given by the total ADM mass-energy of the black hole. This is a physically sound assumption since the internal energy is the energy necessary to 'create' a given system and
the ADM mass is the total energy of the black hole.

As a final consideration for this section, we study the ultraviolet term $R\rightarrow 0$. In the expression (\ref{65}) we have subtracted, as
usual, the 'ultraviolet' term $U_0=2\pi c\hbar N/R$.
This term is obviously diverging for $R\rightarrow 0$ at fixed finite $N$ and it is thus expected
to play a role when the Planck length $L_P$ is reached. This term becomes dominant when $R\sim L_P$.
By inserting this term in the (\ref{65}), we see that the leading term of the (\ref{71}) for
$R\sim L_P$ is given by:
\begin{equation}
R\simeq 2\sqrt{\pi N}\;L_P.
\label{q1}
\end{equation}
At the Planck length it is expected that
the spacetime geometry is not more given by the usual one but substituted with a non-commutative
spacetime thus evidencing a quantum structure for the
spacetime, i.e. a quantum spacetime \cite{a3,a4}.
For a quantum spacetime, uncertainty relations involving coordinates come into action (see \cite{a3} and the generalization
in \cite{a5,a6}). In particular, for spherical geometries, a minimal uncertainty $\Delta R$
for the radial coordinates does arise of the order of the
Planck length \cite{a5}, i.e. $\Delta R\geq \sqrt{3} L_P$: the expression (\ref{q1}) is in agreement with this inequality. In particular,
we may think that the minimum radius $R_{min}$ can be obtained in (\ref{q1}) by taking a single ($N=1$) oscillator thus obtaining
$R_{min}\simeq 2\sqrt{\pi}\;L_P$. This suggests that a quantum spacetime in a state with spherical shape can be seen as composed of
spherical balls each filled with a single
harmonic oscillator made of a single oscillating graviton.

\section{A summation over $m$}

In section 3 the analogy between the Regge-Wheeler and Zerilli equations
with the Schrodinger one allowed us to physically motivate the
(\ref{52}) as representing the spectrum for trapped gravitons. By forcing this interpretation, we can promote the indices
$\{\ell,m,n\}$ to quantum numbers representing quantum states.
In this section we study the case where the summation over the azimuthal index $m$ is performed and as a consequence instead of the
(\ref{60}) we have
\begin{equation}
Z_g=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{m=+\ell}\sum_{n=0}^{\infty} e^{-\beta\hbar{\omega}_{\ell n}}=
\sum_{\ell=2}^{\infty}\sum_{n=0}^{\infty}(2\ell+1)e^{-\beta\hbar{\omega}_{\ell n}}.
\label{a1}
\end{equation}
To perform the summation (\ref{a1}), we can use the following usual formula that is fulfilled for a geometric series for
$x>0$;
\begin{equation}
f(x)=\frac{1}{1-e^{-x}}=\sum_{i=0}^{\infty}e^{-ix},\;\frac{df(x)}{dx}={\left(\frac{1}{1-e^{-x}}\right)}_{,x}=
-\sum_{i=1}^{\infty}ie^{-ix}.
\label{a2}
\end{equation}
With the (\ref{a2}) we obtain for $Z_g$, with $X=\frac{c\pi\beta\hbar}{R}$,
\begin{equation}
Z_g=\frac{e^{-2X}(3-2 e^{-\frac{X}{2}})}{(1-e^{-X}){(1-e^{-\frac{X}{2}})}^2}.
\label{a3}
\end{equation}
Formula (\ref{a3}) contains extra terms with respect to (\ref{60}). For the total free energy $F_T$
we obtain:
\begin{equation}
F_T=-K_BN T\left[-2X+\ln\left(3-2e^{-\frac{X}{2}}\right)-\ln\left(1-e^{-X}\right)-
2\ln\left(1-e^{-\frac{X}{2}}\right)\right].
\label{a4}
\end{equation}
For the entropy $S$ we obtain:
\begin{equation}
S=-\frac{F_T}{T}+\frac{c\pi N\hbar e^{-\frac{X}{2}}}{TR}
\frac{\left[2+4 e^{-\frac{X}{2}}-3e^{-X}\right]}{\left(3-2 e^{-\frac{X}{2}}\right)\left(1-e^{-X}\right)}.
\label{a5}
\end{equation}
First of all, we can study the internal energy $U$ (after subtracting the zero energy value):
\begin{equation}
U=\frac{c\pi N\hbar}{R}\left[-\frac{1}{\left(3e^{\frac{X}{2}}-2\right)}+\frac{1}{\left(e^X-1\right)}
+\frac{1}{\left(e^{\frac{X}{2}}-1\right)}\right].
\label{a6}
\end{equation}
Note the negative term in the (\ref{a6}). Nevertheless, we have $U>0\;\forall X>0$.
In the limit $X\rightarrow 0$ ($T\rightarrow\infty$ at fixed $R$ or $R\rightarrow\infty$ at $T$ fixed) we have
$U\rightarrow 3NK_B T$ ($2NK_BT$ for the (\ref{65})), whereas for
$X\rightarrow\infty$ we always have $U\rightarrow 0$.

Also for the (\ref{a1}) we can apply the same arguments of section 6. In particular, the (\ref{70}) is still valid and the
(\ref{71}) becomes:
\begin{equation}
R=\sqrt{\pi N}L_P\sqrt{\frac{2}{\left[e^{\frac{{2\pi}^2}{\alpha}}-1\right]}+
	\frac{2}{\left[e^{\frac{4{\pi}^2}{\alpha}}-1\right]}-\frac{2}{\left[3e^{\frac{{2\pi}^2}{\alpha}}-2\right]}}.
\label{a7}
\end{equation}
By following the same arguments of section 6, we found that also with the (\ref{a1}) we obtain the
Bekenstein-Hawking entropy only for a value, i.e. $\alpha\simeq 8.48$.

Summarizing, the summation over $m$ leads to corrections with respect to the formulas of sections 4-6, but the physical features remain substantially unchanged.

\section{Conclusions and final remarks}

In this paper, our attention focused on the analogy between the Regge-Wheeler
(\ref{23}) and Zerilli (\ref{f3}) equations with the Schrodinger one. Obviously, the key question is under what conditions
$Z^{(a,p)}$ behaves, in practice, like a quantum wavefunction. Since we described free massless gravitons, equations
(\ref{23}) and (\ref{f3}) are always 'classical'.
A reasonable way to see the quantum behavior of gravitons is to
confine the gravitational radiation in a finite spherical box of radial radius $R$ and thus impose the appropriate
(Dirichlet) boundary conditions,
i.e. $Z_{\ell m}^{(a,p)}(t,R)=0$. As a result of this condition, a discrete spectrum arises dependent on the Legendre index
$\ell$ and the integer $n$ of usual harmonic oscillators. The physical situation can be assimilated to the one of a certain gravitational
radiation \cite{21} in thermal equilibrium within a given cavity. The index $\ell$ can be as well be interpreted as representing
the species of gravitons present after perturbation. In this case,
by calculating the canonical partition function, as a result of a
summation over $\ell$ and $n,$ the so obtained distribution is a superposition of two Planck distributions
(\ref{62}) with proportional frequencies. This is an interesting result since differs from the single Planck distribution of photons in thermal equilibrium within
a cavity.\\
The other possibility is to promote all the indices (not only $n$) $\{\ell,m,n\}$ to quantum numbers and a summation
with respect to $m$ becomes necessary. This case is studied in section 7 with formulas similar to the ones obtained in sections
3-6 and with physical features substantially unchanged.

As an interesting application, we used our results to obtain the Bekenstein-Hawking entropy in a simple and physically clear way.
The only assumption has been that the internal energy $U$ of the black hole is provided by its ADM rest energy, a very reasonable assumption.
As a result, in order to have consistency with the behavior $U\sim R$ and an entropy proportional to the area of the event horizon,
we must have $N\sim R^2$ and as a consequence we obtain
for the internal temperature $T_i$ of the black hole $T_i\simeq 2.225 T_h$ ($T_i\simeq 8.48 T_h$ by summing up with respect to $m$).
These rather intriguing results have a simple physical
interpretation. The internal degrees of freedom of a black hole are provided by gravitons that can be seen as an ensemble of $N$
harmonic oscillators. Every single graviton contributes to the total entropy by a quantity proportional to the Planck area
$A_P=4\pi L_P^2$, in agreement with the holographic nature of the black hole entropy.

Concerning the approach used in this paper, some analogy can be done with the papers quoted in
\cite{E3,CC1,CC2,CC3}. There, the author uses the correspondence between Hawking radiation and black hole quasi-normal modes
that are in turn interpreted as quantum levels in a semi-classical model inspired by the Bohr approach used in 1913.
The interesting consequences of these setups are that black hole entropy is proportional to the black hole excited states and
quite remarkably represent a possible reasonable solution to the black hole information paradox.
Since we attempt to describe the interior of the black hole near the horizon, while in the papers \cite{E3,CC1,CC2,CC3}
the exterior near horizon region is considered, our approach could be seen as complementary to the one present in
the quoted above papers.\\
Moreover, note that the discrete spectrum formula (\ref{52}) for a graviton of species $\ell$
can be rewritten in the form
\begin{equation}
{\omega}_{\ell n}\frac{R}{c}\simeq \left(1+\frac{\ell}{2}+n\right)\pi
\label{af}
\end{equation}
that, apart from the term $\left(1+\frac{\ell}{2}\right)\pi$, is related to the famous
Bohr quantization rule $\oint pdq=nh$ for the coordinate $q$ and its conjugate variable $p$ in the form
\begin{equation}
\oint pdq=\left(1+\frac{\ell}{2}+n\right)\pi h,
\label{af2}
\end{equation}
with $q=R$ and $p=\hbar{\omega}_{\ell n}/c$.

As a final consideration, notice that to obtain our results we have not invoked particular quantum gravity theories,
and for this reason this makes our result in some sense 'phenomenological' and thus physically sound. A possible further
development of the argument presented in this paper is to explain in a statistical sense the recently proposed generalization for the black
hole entropy in a cosmological context \cite{30,31,32,33} and in particular the recent proposal
\cite{34} concerning the nature of the cosmological constant in terms of 'slow' gravitons.

\section*{Acknowledgements}
I would like to thank Luca Tomassini for useful discussions concerning non commutative spacetimes and the anonymous referee who
suggested the references \cite{E3,CC1,CC2,CC3} in relation with the approach used in this paper.


\section{introduction}

The problem of communication for omniscience (CO) was originally formulated in \cite{Csiszar2004}. It is assumed that there are a finite number of users in a system. Each of them observes a distinct component of a discrete multiple correlated source in private. The users are allowed to exchange their observations over public authenticated noiseless broadcast channels. The purpose is to attain \textit{omniscience}, the state that each user obtains all the components in the entire multiple source. The CO problem in \cite{Csiszar2004} is based on an \textit{asymptotic model} where the length of observation sequence is allowed to approach infinity. Whereas the finite linear source model \cite{Chan2011ITW} and packet model in the coded cooperative data exchange (CCDE) \cite{Roua2010,SprintRand2010,CourtIT2014} can be considered as the \textit{non-asymptotic model} where the number of observations is finite and the communication rates are restricted to be integral.

While, in the majority of the studies, e.g. \cite{Csiszar2004,CourtIT2014,MiloIT2015}, the omniscience is attained in a one-off manner, the concept of successive omniscience (SO) is proposed in \cite{ChanSuccessiveIT,Ding2015NetCod}. The idea is to let the omniscience be achieved in a successive manner: attain the local omniscience in a user subset before the global omniscience. It is shown in \cite{ChanSuccessiveIT} that we can attain the local omniscience in a \textit{complimentary} set, a user subset that has a multivariate mutual information (MMI) no less than the MMI in the entire system, while still keep the overall communication rates for CO minimized.
SO is also an attractive idea when we want to design a practical method for CO. Firstly, solving the local omniscience problem (in a user subset) is less complex than the global one. Secondly, in a large scale system where the global omniscience takes a long time, it is practical to let a small group of users attain the local omniscience first so that they can be treated as a super user, which reduces the dimension of the CO problem. Thirdly, when the system parameters vary, e.g., when some users in the original system move out of the communication range,\footnote{This could happen in a large scale CCDE system where the users are mobile clients.} the solution to SO is optimal up-to-the-date.

While the implementation of SO boils down to the problem of how to determine a complimentary user subset, the study in \cite{ChanSuccessiveIT} does not provide a practical solution and the iterative merging algorithm proposed in \cite{Ding2015NetCod} only applies to the CCDE systems. In this paper, we propose an efficient algorithm for searching the complimentary subset for SO in both asymptotic and non-asymptotic models. First, the necessary and sufficient condition for a user subset to be complimentary in \cite{ChanSuccessiveIT} is converted to the one that is conditioned on the value of the Dilworth truncation. Based on this condition, we propose an algorithm for searching the complimentary user subset for SO, which is called the CompSetSO algorithm. While running this algorithm still requires the value of the minimum sum-rate in both asymptotic and non-asymptotic models, we derive a suffient condition for a user subset to be complimentary so that knowing only the lower bound on the minimum sum-rate is sufficient. We show that, based on this lower bound, the CompSetSO algorithm either searches a complimentary subset or returns a global-omniscience-achievable rate vector with the minimum sum-rate in $O(|V| \cdot \SFM(|V|))$ time. Here, $V$ denotes the user set and $\SFM(|V|)$ denotes the complexity of minimizing a submodular set function that is defined on $2^V$. Finally, an example is presented to show how to attain omniscience in multi-stages of SO in a CCDE system by adopting random linear network coding (RLNC) scheme \cite{Ho2006RLNC}.

\section{System Model}
\label{sec:system}

Let $V$ with $\Card{V}>1$ be a finite set that contains the indices of all users in the system. We call $V$ the \textit{ground set}. Let $\RZ{V}=(\RZ{i}:i\in V)$ be a vector of discrete random variables indexed by $V$. For each $i\in V$, user $i$ privately observes an $n$-sequence $\RZ{i}^n$ of the random source $\RZ{i}$ that is i.i.d.\ generated according to the joint distribution $P_{\RZ{V}}$. We allow the users exchange their sources directly so as to let all users in $V$ recover the source sequence $\RZ{V}^n$. The state that each user obtains the total information in the entire multiple source is called \textit{omniscience}, and the process that users communicate with each other to attain omniscience is called \textit{communication for omniscience (CO)} \cite{Csiszar2004}. For $X \subseteq V$, the \textit{(local) omniscience} in $X$ refers to the state that all users $i \in X$ recover the sequence $\RZ{X}^n$. The \textit{global omniscience} in the ground set $V$ is a special case of the local omniscience when $X = V$.

\subsection{Minimum Sum-rate}

Consider the local omniscience problem in $X \subseteq V$. Let $\rv_X=(r_i:i\in X)$ be a rate vector indexed by $X$. We call $\rv_X$ an \textit{achievable rate vector} if the omniscience in $X$ can be attained by letting users communicate with the rates designated by $\rv_X$. Let $r$ be the function associated with $\rv_X$ such that
$$ r(C)=\sum_{i\in C} r_i, \quad \forall C \subseteq X $$
with the convention $r(\emptyset)=0$. We call $r(X)$ the \textit{sum-rate} of $\rv_X$. For $C \subseteq X$, let $H(\RZ{C})$ be the amount of randomness in $\RZ{C}$ measured by Shannon entropy and $H(\RZ{X \setminus C}|\RZ{C})=H(\RZ{X})-H(\RZ{C})$ be the conditional entropy of $\RZ{X \setminus C}$ given $\RZ{C}$. In the rest of this paper, without loss of generality, we simplify the notation $\RZ{C}$ by $C$.

It is shown in \cite{Csiszar2004} that an achievable rate vector $\rv_X$ must satisfy the Slepian-Wolf (SW) constraints: $r(C) \geq H(C|X \setminus C), \forall C \subset X$. The interpretation is: To attain the omniscience in $X$, the total amount of information sent from user set $C$ should be at least equal to what is missing in $X \setminus C$. We have the set of all achievable rate vectors being
$$ \RRCO(X)=\Set{ \rv_X\in\Real^{|X|} \colon r(C) \geq H(C|X\setminus C),\forall C \subset X }.$$

In an asymptotic model, we study the CO problem by considering the asymptotic limits as the \emph{block length} $n$ goes to infinity so that the communication rates could be real or fractional; In a non-asymptotic model, the block length $n$ is restricted to be finite and the communication rates are required to be integral. $\RACO(X)=\min\Set{r(X) \colon \rv_X\in \RRCO(X)}$ and $\RNCO(X)=\min\Set{r(X) \colon \rv_X\in \RRCO(X) \cap \Z^{|X|} }$ are the minimum sum-rates for attaining omniscience in $X$ in the asymptotic and non-asymptotic models, respectively. $\RRACO(X)=\Set{\rv_X\in \RRCO(X) \colon r(X)=\RACO(X)} $ and $\RRNCO(X)=\Set{\rv_X \in \RRCO(X) \cap \Z^{|X|} \colon r(X)=\RNCO(X)}$ are the corresponding optimal rate vector sets for the asymptotic and non-asymptotic models, respectively.

For $X \subseteq V$, let $\Pi(X)$ be the set of all partitions of $X$ and $\Pi'(X) = \Pi(X) \setminus \Set{X}$. It is shown in \cite{Csiszar2004,ChanMMI,Ding2016IT} that
    \begin{equation} \label{eq:MinSumRate}
            \RACO(X) = \max_{\Pat \in \Pi'(X)} \sum_{C \in \Pat} \frac{H(X) - H(C)}{|\Pat|-1}
    \end{equation}
and $\RNCO(X) = \lceil \RACO(X) \rceil$.
The maximization problem \eqref{eq:MinSumRate} can be solved and an optimal rate vector in $\RRACO(X)$ or $\RRNCO(X)$ can be determined in $O(|X|^2 \cdot \SFM(|X|))$ time \cite{Ding2016IT}.\footnote{$O(\SFM(|X|))$ is the complexity of solving the problem $\min\Set{f(C) \colon C \subseteq X}$ for submodular set function $f$, which is strongly polynomial \cite{Fujishige2005}.}

\begin{example} \label{ex:aux}
Consider the system where $V=\{1,\dotsc,5\}$ and each user observes respectively
    \begin{equation}
        \begin{aligned}
            \RZ{1} & = (\RW{a},\RW{b},\RW{c},\RW{d},\RW{f},\RW{g},\RW{i},\RW{j}),  \\
            \RZ{2} & = (\RW{a},\RW{b},\RW{c},\RW{f},\RW{i},\RW{j}),\\
            \RZ{3} & = (\RW{e},\RW{f},\RW{h},\RW{i}),   \\
            \RZ{4} & = (\RW{b},\RW{c},\RW{e},\RW{j}), \\
            \RZ{5} & = (\RW{b},\RW{c},\RW{d},\RW{h},\RW{i}),
        \end{aligned}  \nonumber
    \end{equation}
where $\RW{j}$ is an independent uniformly distributed random bit. In the corresponding CCDE system, each $\RW{j}$ represents a packet and $\RZ{i}$ denotes the packets received by mobile client $i$ from a base station. The mobile clients in $V$ transmit linear combinations of $\RZ{i}$s by some network coding scheme, e.g., \cite{Roua2010}, over noiseless peer-to-peer channels in order to recover all packets in $\RZ{V}$. In this system, $\RACO(V)= \frac{13}{2}$ and $\RNCO(V) = 7$. We have $(\frac{9}{2},0,\frac{1}{2},\frac{1}{2},1) \in \RRACO(V)$ and $(5,0,1,1,0) \in \RRNCO(V)$ being one of the optimal rate vectors in the asymptotic and non-asymptotic models, respectively.\footnote{The optimal rate vector is not unique, i.e., $\RRACO(V)$ and $\RRNCO(V)$ are not singleton, in general.}
\end{example}

\section{Successive Omniscience}

The idea of successive omniscience (SO) is proposed in \cite{ChanSuccessiveIT,Ding2015NetCod}, which allows the omniscience to be achieved in a successive manner: first attain the local omniscience in a user subset $X$; then solve the global omniscience problem in $V$ by assuming that all the users $i \in V \setminus X$ have obtained the information in the communications for achieving the local omniscience in $X$.

\subsection{Complimentary User Subset}

Let $X \subseteq V$ such that $|X| > 1$ be a non-singleton user subset. $X$ is called a \textit{complimentary} subset (for SO) if the local omniscience in $X$ can be achieved first without increasing the minimum sum-rate, $\RACO(V)$ and $\RNCO(V)$ in the asymptotic and non-asymptotic models, respectively, for the global omniscience in $V$\cite{ChanSuccessiveIT}.
To be more specific, take the asymptotic model for example. For a non-singleton subset $X \subseteq V$, let $\rv_X \in \RRACO(X)$ be an optimal rate vector for attaining the local omniscience in $X$. Let $\rv'_V = (r'_i \colon i \in V )$ where $r'_i = r_i$ if $i \in X$ and $r'_i=0$ otherwise. If there exists a rate vector $\rv''_V \in \RealP^{|V|}$ such that $\rv'_V + \rv''_V \in \RRACO(V)$, then local omniscience in $X$ is complimentary. Similarly, in the non-asymptotic model, let $\rv'_V $ be constructed in the same way by $\rv_X \in \RRNCO(X)$ for $X \subseteq V$. $X$ is complimentary if there exists a rate vector $\rv''_V \in \ZP^{|V|}$ such that $\rv'_V + \rv''_V \in \RRNCO(V)$.

\begin{example}\label{ex:auxSO}
    Consider the asymptotic model for the system in Example~\ref{ex:aux}. The minimum sum-rate for achieving the local omniscience in $\Set{1,2}$ is $\RACO(\Set{1,2}) = 2$ and $\rv_{\Set{1,2}} = (r_1,r_2) = (2, 0) \in \RRACO(\Set{1,2})$ is an optimal rate vector. In this case, we have $\rv'_V = (2,0,0,0,0)$. After the users transmit $\rv'_V$ for achieving the local omniscience in $\Set{1,2}$, we can let them transmit $\rv''_V = (\frac{5}{2},0,\frac{1}{2},\frac{1}{2},1)$, which attains the global omniscience in $V = \Set{1,\dotsc,5}$. Therefore, $\Set{1,2}$ is a complimentary user subset. Here, $\rv'_V + \rv''_V = (\frac{9}{2},0,\frac{1}{2},\frac{1}{2},1) \in \RRACO(V)$ is an optimal rate vector as shown in Example~\ref{ex:aux}. In other words, the optimal rate vector $(\frac{9}{2},0,\frac{1}{2},\frac{1}{2},1)$ for the omniscience in $V = \Set{1,\dotsc,5}$ can be implemented in an SO manner so that the local omniscience in $\Set{1,2}$ can be achieved before the global omniscience.

    For the non-asymptotic model, $\Set{1,2}$ is also a complimentary user subset since SO can be done by rate vectors $\rv'_V = (2,0,0,0,0)$ and $\rv''_V = (3,0,1,1,0)$, which first achieve local omniscience in $\Set{1,2}$ and then global omniscience in $V = \Set{1,\dotsc,5}$. Here, $\rv'_V + \rv''_V = (5,0,1,1,0) \in \RRNCO(V)$ is an optimal rate vector as shown in Example~\ref{ex:aux}.
\end{example}

\subsubsection{necessary and sufficient condition}

The necessary and sufficient condition for a user subset $X$ to be complimentary is derived in \cite{ChanSuccessiveIT} for both asymptotic and non-asymptotic models, which is stated as follows.

\begin{theorem}[necessary and sufficient condition {\cite[Theorems 4.2 and 5.2]{ChanSuccessiveIT}}]\label{theo:SO}
    In an asymptotic model, a non-singleton user subset $X \subset V$ is complimentary if and only if $H(V) - H(X) + \RACO(X) \leq \RACO(V)$; In a non-asymptotic model, a non-singleton user subset $X \subset V$ is complimentary for SO if and only if $H(V) - H(X) + \RNCO(X) \leq \RNCO(V)$.\footnote{Let $I(X)$ denote the multivariate mutual information (MMI) in $X$. In {\cite[Theorems 4.2 and 5.2]{ChanSuccessiveIT}}, the necessary and sufficient condition for $X$ to be complimentary is $I(X) \geq I(V)$ and $\lfloor I(X) \rfloor \geq \lfloor I(V) \rfloor$ for the asymptotic and non-asymptotic models. They can be converted to the conditions in Theorem~\ref{theo:SO} via the dual relationships: $\RACO(V) = H(V)- I(V) $ and $\RNCO(V) = H(V)- \lfloor I(V) \rfloor$ \cite{ChanMMI,Court2011}.
    Also, since the ground set $V$ is always a complimentary subset, we restrict our attention to the non-singleton proper subsets $X$ of $V$ that are complimentary.} \hfill\IEEEQED
\end{theorem}

Here, $H(V) - H(X)$ is the amount of information that is missing in user subset $X$, the omniscience of which only relies on the transmissions from the users in $V \setminus X$. If we let the users in $X$ attain local omniscience with the minimum sum-rate $\RACO(X)$, the users in $V \setminus X$ are required to transmit at least $H(V) - H(X)$ for attaining the global omniscience. Then, the total number of transmissions is no less than $H(V) - H(X) + \RACO(X)$. If $H(V) - H(X) + \RACO(X)>\RACO(V)$, the global omniscience is not achievable by the minimum sum-rate $\RACO(V)$ if we allow the users in $X$ to attain the local omniscience first. The condition $H(V) - H(X) + \RNCO(X) \leq \RNCO(V)$ for the non-asymptotic model in Theorem~\ref{theo:SO} can be interpreted in the same way.

However, Theorem~\ref{theo:SO} cannot be directly applied for determining a complimentary subset since the power set $2^V$ is exponentially large in $|V|$. In the following context, we convert Theorem~\ref{theo:SO} to the conditions on the Dilworth truncation and propose a polynomial time algorithm for searching a complimentary user subset for SO.

For $0 \leq \alpha \leq H(V)$, let
$$ \FuHash{\alpha}(X) = \begin{cases}
                       0, & \mbox{if $X = \emptyset$} \\
                       \alpha - H(V) + H(X), & \mbox{otherwise}.
                     \end{cases}$$
$\FuHashHat{\alpha}(X) = \min_{\Pat \in \Pi(X)} \sum_{C \in \Pat} \FuHash{\alpha}(C), \forall X \subseteq V$ is the \textit{Dilworth truncation} of $\FuHash{\alpha}$ \cite{Dilworth1944}.

\begin{corollary} \label{coro:SOEQV1}
    In an asymptotic model, a non-singleton user subset $X \subset V$ is complimentary for SO if and only if $\FuHash{\RACO(V)}(X) = \FuHashHat{\RACO(V)}(X)$; In a non-asymptotic model, a non-singleton user subset $X \subset V$ is complimentary for SO if and only if $\FuHash{\RNCO(V)}(X) = \FuHashHat{\RNCO(V)}(X)$. \hfill \IEEEQED
\end{corollary}

The proof of Corollary~\ref{coro:SOEQV1} is in Appendix~\ref{app:coro:SOEQV1}

\begin{example} \label{ex:auxSOCompSet}
    For the system in Example~\ref{ex:aux}, we have
     \begin{equation}
        \begin{aligned}
            &\Set{X \subset V \colon |X| > 1, H(V) - H(X) + \RACO(X) \leq \RACO(V)} \\
            &\qquad = \Set{X \subset V \colon |X|>1, \FuHash{\RACO(V)}(X) = \FuHashHat{\RACO(V)}(X) } \\
            &\qquad = \big\{ \Set{1,2},\Set{1,5},\Set{1,2,5},\Set{1,3,4,5} \big\},
        \end{aligned} \nonumber
     \end{equation}
     being all complimentary subsets in the asymptotic model and
     \begin{equation}
        \begin{aligned}
            &\Set{X \subset V \colon |X| > 1, H(V) - H(X) + \RNCO(X) \leq \RNCO(V)} \\
            & \qquad = \Set{ X \subset V \colon |X|>1, \FuHash{\RNCO(V)}(X) = \FuHashHat{\RNCO(V)}(X) } \\
            & \qquad = \big\{ \Set{1,2},\Set{1,4},\Set{1,5},\Set{2,4},\Set{2,5},\Set{1,2,4},\\
            & \qquad\qquad \Set{1,2,5},\Set{1,3,4},\Set{1,3,5},\Set{1,4,5},\Set{2,3,4},\\
            & \qquad\qquad \Set{2,3,5},\Set{2,4,5},\Set{1,2,3,4},\Set{1,2,3,5},\\
            & \qquad\qquad \Set{1,2,4,5},\Set{1,3,4,5},\Set{2,3,4,5} \big\},
        \end{aligned} \nonumber
     \end{equation}
     being all complimentary subsets in the non-asymptotic model.
\end{example}

Then, the task reduces to finding a subset $X$ such that $\FuHash{\alpha}(X) = \FuHashHat{\alpha}(X)$, where $\alpha=\RACO(V)$ and $\alpha=\RNCO(V)$ for the asymptotic and non-asymptotic models, respectively. It can be converted to a minimization problem
               \begin{equation} \label{eq:SFMAlgo}
                    \min_{X \subseteq V_i \colon i \in X} \Set{ \FuHash{\alpha}(X) - r(X) } ,
                \end{equation}
where $i \in V$ and $V_i = \Set{1,\dotsc,i}$, based on which, we propose an algorithm for searching for the complimentary subset for SO (CompSetSO) in Algorithm~\ref{algo:CompSetSO}.

        \begin{algorithm} [t]
           \caption{complimentary subset for SO (CompSetSO)}
	       \label{algo:CompSetSO}
	       \small
	       \SetAlgoLined
	       \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
	       \SetKwFor{For}{for}{do}{endfor}
            \SetKwRepeat{Repeat}{repeat}{until}
            \SetKwIF{If}{ElseIf}{Else}{if}{then}{else if}{else}{endif}
	       \BlankLine
           \Input{the ground set $V$, an oracle that returns the value of $H(X)$ for a given $X \subseteq V$ and $\alpha$, which is determined based on Theorem~\ref{theo:algo:CompSetSO} or Theorem~\ref{theo:algo:CompSetSOExt}.}
	       \Output{$\hat{X}$ which is a complimentary user subset for SO}
	       \BlankLine
            $ r_{1} \leftarrow \FuHash{\alpha}(\Set{1})$ and $r_i \leftarrow \alpha - H(V), \forall i \in V \setminus \Set{1}$\;
            \For{$i=2$ \emph{\KwTo} $|V|$}{

                \eIf{there exists a non-singleton minimizer $\hat{X}$ of $\min_{X \subseteq V_i \colon i \in X} \Set{ \FuHash{\alpha}(X) - r(X) }$ such that $\hat{X} \subset V$}{
                    terminate iteration and return $\hat{X}$;}{

                    $ r_i \leftarrow r_i + \min_{X \subseteq V_i \colon i \in X} \Set{ \FuHash{\alpha}(X) - r(X)}$\; }
            }
	   \end{algorithm}

\begin{theorem} \label{theo:algo:CompSetSO}
    For the CompSetSO algorithm in Algorithm~\ref{algo:CompSetSO}, the output $\hat{X}$ is a complimentary user subset for the asymptotic and non-asymptotic models if the input $\alpha = \RACO(V)$ and $\alpha = \RNCO(V)$, respectively; If there is no output, there does not exist a complimentary subset for SO. \hfill\IEEEQED
\end{theorem}

The proof of Theorem~\ref{theo:algo:CompSetSO} is in Appendix~\ref{app:theo:algo:CompSetSO}.

\begin{example}  \label{ex:auxCompSetSOMinSumRate}
    We apply the CompSetSO algorithm to the asymptotic model of the system in Example~\ref{ex:aux} by inputting $\alpha=\RACO(V) = \frac{13}{2}$. It can be shown that for $i = 2$, $\hat{X} = \Set{1,2}$ is returned as a complimentary subset. For the non-asymptotic model, by inputting $\alpha = \RNCO(V) = 7$, we get $\hat{X} = \Set{1,2}$ returned as a complimentary subset for $i = 2$.
\end{example}

\subsubsection{Sufficient Condition}

In Theorem~\ref{theo:algo:CompSetSO}, knowing the value of the minimum sum-rate, $\RACO(V)$ or $\RNCO(V)$, is a prerequisite. However, if we obtain the value of $\RACO(V)$ or $\RNCO(V)$, say, by the modified decomposition algorithm in \cite{Ding2016IT} or the deterministic algorithms in \cite{CourtIT2014,MiloIT2015}, we necessarily know an optimal rate vector in $\RRACO(V)$ or $\RRNCO(V)$ for the global omniscience, in which case solving a local omniscience problem in a user subset may not be necessary. It is also not consistent with the advantage of SO that the local omniscience problem is less complex than the global one.
So, the question is: Can we find a complimentary user subset for SO without knowing the minimum sum-rate? The answer is yes. In this section, we derive the sufficient conditions for a subset to be complimentary such that the value of $\alpha$ in the CompSetSO algorithm can be relaxed from the exact value of $\RACO(V)$ or $\RNCO(V)$ to a lower bound on $\RACO(V)$ or $\RNCO(V)$ that can be obtained in $O(|V|)$ time.

\begin{lemma}\label{lemma:SOEQV2}
    In an asymptotic model, a non-singleton user subset $X \subset V$ is complimentary if $\FuHash{\alpha}(X) = \FuHashHat{\alpha}(X)$ for $\alpha = \sum_{i \in V} \frac{H(X) - H(\Set{i})}{|V| - 1}$; In a non-asymptotic model, a non-singleton user subset $X \subset V$ is complimentary if $\FuHash{\alpha}(X) = \FuHashHat{\alpha}(X)$ for $\alpha = \big\lceil \sum_{i \in V} \frac{H(X) - H(\Set{i})}{|V| - 1} \big\rceil$. \hfill \IEEEQED
\end{lemma}
The proof of Lemma~\ref{lemma:SOEQV2} is in Appendix~\ref{app:lemma:SOEQV2}. The values of $\alpha$ in Lemma~\ref{lemma:SOEQV2} are the lower bounds on $\RACO(V)$ and $\RNCO(V)$ for asymptotic and non-asymptotic models, respectively \cite[Theorem IV.1]{Ding2016IT}.

\begin{theorem} \label{theo:algo:CompSetSOExt}
    The CompSetSO algorithm in Algorithm~\ref{algo:CompSetSO} returns a complimentary user subset $\hat{X}$ for the asymptotic and non-asymptotic models if the input $\alpha = \sum_{i \in V} \frac{H(V) - H(\Set{i})}{|V| - 1}$ and $\alpha = \big\lceil \sum_{i \in V} \frac{H(V) - H(\Set{i})}{|V| - 1} \big\rceil$, respectively; If there is no output, there does not exist complimentary subset for SO.  \hfill\IEEEQED
\end{theorem}

The proof of Theorem~\ref{theo:algo:CompSetSOExt} is in Appendix~\ref{app:theo:algo:CompSetSOExt}. The case when there is no complimentary subset in both Theorems~\ref{theo:algo:CompSetSO} and \ref{theo:algo:CompSetSOExt} is given in Appendix~\ref{app:NoCompSetSO}.

\begin{example} \label{ex:auxCompSetSOLB}
    For the asymptotic model of the system in Example~\ref{ex:aux}, by inputting $\alpha = \sum_{i \in V} \frac{H(X) - H(\Set{i})}{|V| - 1} = \frac{23}{4}$ in the CompSetSO algorithm, we get output $\hat{X} = \Set{1,2}$ as a complimentary subset. By inputting $\alpha = \big\lceil \sum_{i \in V} \frac{H(X) - H(\Set{i})}{|V| - 1} \big\rceil = 6$ in the CompSetSO algorithm, we get output $\hat{X} = \Set{1,2}$ as a complimentary subset for the non-asymptotic model.
\end{example}

\subsubsection{Complexity}

In the CompSetSO algorithm, \eqref{eq:SFMAlgo} is a submodular function minimization (SFM) problem,\footnote{The submodularity of \eqref{eq:SFMAlgo} is proved in \cite{Ding2016IT} based on the submodularity of the entropy function $H$.} which can be solved in $O(\SFM(|V|))$ time. Let $\alpha$ be initiated according to Theorem~\ref{theo:algo:CompSetSOExt}. Then, in an asymptotic or non-asymptotic model, if there exists a complimentary user subset $\hat{X}$ for SO, it can be found in $O(|V| \cdot \SFM(|V|))$ time. Consider the situation when there is no complimentary user subset for SO. The CompSetSO algorithm completes in $O(|V| \cdot \SFM(|V|))$ time without any output. But, in this case, according to \cite[Theorem V.1]{Ding2016IT}, $\rv_V$ is updated to an optimal rate vector in $\RRACO(V)$ and $\RRNCO(V)$ for the asymptotic and non-asymptotic models, respectively. In summary, the CompSetSO algorithm determines either a complimentary user subset for SO or an optimal rate vector in $O(|V| \cdot \SFM(|V|))$ time.

\section{Multi-stage Successive Omniscience}

We show an example of multi-stage SO in a CCDE system.

\begin{example} \label{ex:auxMultiSO}
    Consider the system in Example~\ref{ex:aux} as a CCDE system where the linear combinations of packets are transmitted by random linear network coding (RLNC) scheme\cite{Ho2006RLNC}. For example, if $r_i = 2$, then user $i$ broadcasts the linear coding $S_i = \gamma^{\intercal} \Rz{i} = \sum_{j \in \Rz{i}} \gamma_j \RW{j} $ twice. Here, $\gamma = (\gamma_j : \RW{j} \in \Rz{i})$ and, at each broadcast, $\gamma_j$ is randomly chosen from a Galois field $\F_q$ with $q > H(V) \cdot |V|$. If $r_i = \frac{3}{2}$, then each packet $\RW{j} \in \RZ{i}$ is broken into two chunks: $\RW{j}^{(1)}$ and $\RW{j}^{(2)}$ and user $i$ broadcast $S_i = \sum_{j \in \RZ{i}} (\gamma_j^{(1)} \RW{j}^{(1)} + \gamma_j^{(2)} \RW{j}^{(2)})$ for six times. At each broadcast, each $\gamma_j^{(1)}$ and $\gamma_j^{(2)}$ are randomly chosen from a Galois field $\F_q$ with $q > 2 \cdot H(V) \cdot |V|$.

    For the asymptotic model, the global omniscience can be achieved by three stages of SO: By transmitting the rates $\rv'_V = (2,0,0,0,0)$, $\rv''_V = (2,0,0,0,1)$ and $\rv'''_V = (\frac{1}{2},0,\frac{1}{2},\frac{1}{2},0)$, the omniscience is achieved in $\Set{1,2}$, $\Set{1,2,5}$ and $\Set{1,2,3,4,5}$ in sequence. Here, we have $\rv'_V + \rv''_V + \rv'''_V = (\frac{9}{2},0,\frac{1}{2},\frac{1}{2},1) \in \RRACO(V)$, which also means that the optimal rate vector $(\frac{9}{2},0,\frac{1}{2},\frac{1}{2},1)$ can be implemented by three stages of SO. We have shown how to find the complimentary user subset $\Set{1,2}$ by the CompSetSO algorithm in Examples~\ref{ex:auxCompSetSOMinSumRate} and \ref{ex:auxCompSetSOLB}. The complimentary subset $\Set{1,2,5}$ is determined as follows.

    Since local omniscience in $\Set{1,2}$ is attained after the transmission of $\rv'_V$, we can treat $\Set{1,2}$ as a super user and assign a user index $12'$. For each user $i \in V \setminus \Set{1,2}$, we assign a new index $i'$ with $\RZ{i'} = \RZ{i} \cup \Gamma$, where $\Gamma$ contains all the transmissions that is received by user $i$ in the first stage of SO, i.e., all the broadcasts for achieving the local omniscience in $\Set{1,2}$. For the super user $12'$, we have $\RZ{12'} = \RZ{1} \cup \RZ{2}$. We construct the ground set of the new system as $V' = \Set{12',3',4',5'}$ with the observations $\RZ{i'}$ for all $i \in V'$. By applying the CompSetSO algorithm to the new system, we get $\Set{12',5'}$ as a complimentary subset, which corresponds to $\Set{1,2,5}$ in the original system.

    In the same way, for the non-asymptotic model, it can be shown that the global omniscience can be achieved by three stages of SO: By transmitting the rates $\rv'_V = (2,0,0,0,0)$, $\rv''_V = (3,0,0,1,0)$ and $\rv'''_V = (0,0,1,0,0)$, the omniscience is achieved in $\Set{1,2}$, $\Set{1,2,4}$ and $\Set{1,\dotsc,5}$ in sequence. And, $\rv'_V + \rv''_V + \rv'''_V = (5,0,1,1,0) \in \RRNCO(V)$, which also means that the optimal rate vector $(5,0,1,1,0)$ can be implemented by three stages of SO.
\end{example}

\section{Conclusion}
We studied the problem of how to efficiently search a complimentary user subset so that the omniscience of a discrete multiple random source among a set of users can be attained in a successive manner. Based on the existing necessary and sufficient condition for a user subset to be complimentary, we proposed a CompSetSO algorithm, which searches a complimentary subset for SO in both asymptotic and non-asymptotic models. We showed that inputting a lower bound, instead of the exact value, of the minimum sum-rate is sufficient for the CompSetSO algorithm to return either a complimentary subset or an optimal rate vector in $O(|V| \cdot \SFM(|V|))$ time. The CompSetSO algorithm can be implemented recursively so that the omniscience can be attained in multi-stages of SO.

For the future research work, it is worth studying how to implement the SO more efficiently than the CompSetSO algorithm. Also, for the CCDE problem, it would be of interest if the multi-stage SO can be implemented by network coding schemes other than RLNC.

\appendices

\section{Proof of Corollary~\ref{coro:SOEQV1}} \label{app:coro:SOEQV1}
    Based on Theorem~\ref{theo:SO}, we have $ \RACO(X) \leq \RACO(V) - H(V) + H(X) = \FuHash{\RACO(V)}(X) $ being the necessary and sufficient condition for $X$ to be complimentary in the asymptotic model.
    On the other hand, $\RACO(X) \geq \sum_{C \in \Pat} \frac{H(X) - H(C)}{|\Pat| - 1}, \forall \Pat \in \Pi'(X)$. Then,
    we have inequality $\sum_{C \in \Pat} \frac{H(X) - H(C)}{|\Pat| - 1} \leq \RACO(V) - H(V) + H(X), \forall \Pat \in \Pi'(X)$, which is equivalent to $\FuHash{\RACO(V)}(X) \leq \sum_{C \in \Pat} \FuHash{\RACO(V)}(C), \forall \Pat \in \Pi'(X)$, i.e., $\FuHash{\RACO(V)}(X) = \FuHashHat{\RACO(V)}(X)$. In the same way, we can prove that $X$ is complimentary in the non-asymptotic model if and only if $\FuHash{\RNCO(V)}(X) = \FuHashHat{\RNCO(V)}(X)$. \hfill\IEEEQED

\section{Proof of Theorem~\ref{theo:algo:CompSetSO}} \label{app:theo:algo:CompSetSO}

    For $\alpha = \RACO(V)$, let $\hat{X}$ be the user subset returned by the CompSetSO algorithm, i.e., we find a non-singleton user subset $\hat{X} \subset V$ that minimizes $\min\Set{ \FuHash{\alpha}(X) - r(X) \colon \hat{i} \in X \subseteq V_{\hat{i}}}$, for some $\hat{i} \in \Set{1,\dotsc,|V|}$. It also means that we have not found any non-singleton proper subset of $V$ that minimizes $\min\Set{ \FuHash{\alpha}(X) - r(X) \colon i \in X \subseteq V_i}$, i.e., $\Set{i}$ is the only minimizer, for all $i \in \Set{1,\dotsc,\hat{i}-1}$. In each iteration of the CompSetSO algorithm, we have $\rv_V \in P(\FuHash{\alpha},\leq)$, where $P(\FuHash{\alpha},\leq) = \Set{\rv_V \in \Real^{|V|} \colon r(X) \leq \FuHash{\alpha}(X)}$ is the polyhedron of $\FuHash{\alpha}$\cite{Ding2016IT}. Then, for all $X \subseteq V_{\hat{i}-1}$, we have $\sum_{i \in X} \FuHash{\alpha}(\Set{i}) = r(X) \leq \FuHash{\alpha}(X)$, i.e., $\FuHashHat{\alpha}(X) = \sum_{i \in X} \FuHash{\alpha}(\Set{i}) \leq \sum_{C \in \Pat} \FuHash{\alpha}(C)$ for all $\Pat \in \Pi(X)$.

    On the other hand, since $\hat{X} \subset V$ is a non-singleton minimizer of $\min\Set{ \FuHash{\alpha}(X) - r(X) \colon \hat{i} \in X \subseteq V_{\hat{i}}}$, we have $\FuHash{\alpha}(\hat{X}) - r(\hat{X}) \leq \FuHash{\alpha}(C) - r(C), \forall C \subseteq \hat{X} \colon \hat{i} \in C$. Then,
    \begin{equation}
        \begin{aligned}
            \FuHash{\alpha}(\hat{X}) & \leq \FuHash{\alpha}(C) + r(\hat{X} \setminus C)  \\
                                     & = \FuHash{\alpha}(C) + \sum_{i \in \hat{X} \setminus C} \FuHash{\alpha}(\Set{i}) \\ & \leq \FuHash{\alpha}(C) + \sum_{C' \in \Pat} \FuHash{\alpha}(C'),
        \end{aligned} \nonumber
    \end{equation}
    for all $C \subseteq \hat{X}$ such that $\hat{i} \in C$ and all $\Pat \in \Pi(\hat{X})$. Therefore, $\FuHash{\alpha}(\hat{X}) \leq \sum_{C \in \Pat} \FuHash{\alpha}(C)$ for all $\Pat \in \Pi(\hat{X})$. Then, $\FuHash{\alpha}(\hat{X}) = \FuHashHat{\alpha}(\hat{X})$. According to Corollary~\ref{coro:SOEQV1}. $\hat{X}$ is a complimentary user subset in the asymptotic model. In the same way, we can prove the statement for the non-asymptotic model when $\alpha=\RNCO(V)$. \hfill\IEEEQED

\section{Proof of Lemma~\ref{lemma:SOEQV2}} \label{app:lemma:SOEQV2}
    We have $0 \leq \alpha = \sum_{i \in V} \frac{H(X) - H(\Set{i})}{|V| - 1} \leq \RACO(V)$. If $\FuHash{\alpha}(X) = \FuHashHat{\alpha}(X)$, then
    $$ \sum_{C \in \Pat} \FuHash{\RACO(V)}(C) - \FuHash{\RACO(V)}(X) \geq \sum_{C \in \Pat} \FuHash{\alpha}(C) - \FuHash{\alpha}(X) \geq 0,$$
    for all $\Pat \in \Pi(X)$, i.e., $\FuHash{\RACO(V)}(X) = \FuHashHat{\RACO(V)}(X)$. According to Corollary~\ref{coro:SOEQV1}, $X$ is a complimentary subset in the asymptotic model. In the same way, we can prove that $X \subset V$ such that $|X| > 1$ is complimentary in the non-asymptotic model if $\FuHash{\alpha}(X) = \FuHashHat{\alpha}(X)$ for $\alpha = \big\lceil \sum_{i \in V} \frac{H(X) - H(\Set{i})}{|V| - 1} \big\rceil$.

\section{Proof of Theorem~\ref{theo:algo:CompSetSOExt}} \label{app:theo:algo:CompSetSOExt}
    According to Lemma~\ref{lemma:SOEQV2} and by using the same way as in the proof in Theorem~\ref{theo:algo:CompSetSO}, we can show that the output $\hat{X}$ is complimentary. For $\alpha = \sum_{i \in V} \frac{H(V) - H(\Set{i})}{|V| - 1}$, we have $\sum_{i \in V} H(\Set{i}) = |V| H(V) - (|V| - 1) \alpha$. If there is no output of the CompSetSO algorithm, it means
    \begin{equation}
        \begin{aligned}
            & \quad \sum_{i \in V} \FuHash{\alpha}(\Set{i}) - \sum_{C \in \Pat} \FuHash{\alpha}(C) \\
            & = \sum_{i \in V} \big( \alpha - H(V) + H(\Set{i}) \big) - \sum_{C \in \Pat} \big( \alpha - H(V) + H(C) \big) \\
            & = (|V| - |\Pat|)\big( \alpha - H(V) \big) + \sum_{i \in V} H(\Set{i}) - \sum_{C \in \Pat} H(C) \\
            & = \sum_{C \in \Pat} \big( H(V) - H(C) \big) -(|\Pat| - 1) \alpha < 0, \quad \forall \Pat \in \Pi'(V),
        \end{aligned} \nonumber
    \end{equation}
    which is equivalent to
    $$ \alpha > \sum_{C \in \Pat} \frac{H(V) - H(C)}{|\Pat| - 1}, \quad \forall \Pat \in \Pi'(V), $$
    i.e., $\alpha = \RACO(V)$. Also, we have
    \begin{equation}
        \begin{aligned}
            & \quad \FuHash{\alpha}(X) - \sum_{i \in X} \FuHash{\alpha}(\Set{i}) \\
            & = \alpha - H(V) + H(X) - \sum_{i \in X} \big( \alpha - H(V) + H(\Set{i}) \big) \\
            & = H(X) - \sum_{i \in X} H(\Set{i}) - (|X| - 1) \big( \alpha - H(V) \big)  \\
            & = \sum_{i \in X} \big( H(X) - H(\Set{i}) \big) - (|X| - 1) \big( \alpha - H(V) + H(X) \big)\\
            & > 0, \qquad\quad \forall X \subset V \colon |X| > 1,
        \end{aligned}\nonumber
    \end{equation}
    which is equivalent to
    \begin{equation}
        \begin{aligned}
            \RACO(V) & = \sum_{i \in V} \frac{H(X) - H(\Set{i})}{|X| - 1} \\
                     & > \RACO(V) - H(V) + H(X), \ \forall X \subset V \colon |X| > 1,
        \end{aligned}\nonumber
    \end{equation}
    i.e., there is no complimentary user subset for SO in the asymptotic model. In the same way, we can prove the statement for the non-asymptotic model. \hfill\IEEEQED

\section{Nonexistence of Complimentary Subset} \label{app:NoCompSetSO}
In a system, either asymptotic or non-asymptotic model, if there does not exist a complimentary subset for SO, it means the omniscience can not be attained with the minimum sum-rate in a manner such that the local omniscience in some non-singleton subset $X \subset V$ can be attained first. According to Corollary~\ref{coro:SOEQV1}, it happens when the $\Set{\Set{i} \colon i \in V}$ and $\Set{\Set{V}}$ are the only minimizers of $\FuHashHat{\alpha}(V) = \min_{\Pat \in \Pi(V)} \sum_{C \in \Pat} \FuHash{\alpha}(C)$, i.e., when $\Set{\Set{i} \colon i \in V}$ is the only maximizer of $\RACO(V) = \max_{\Pat \in \Pi'(V)} \sum_{C \in \Pat} \frac{H(V) - H(C)}{|\Pat|-1}$. This is the case when the components in $\RZ{V}$ are mutually independent. But, it is not necessary that a system with mutually independent $\RZ{V}$ does not have a complimentary subset.

\begin{example}
Consider the system where $V=\{1,2,3\}$ and each user observes respectively
    \begin{equation}
        \begin{aligned}
            \RZ{1} & = (\RW{a},\RW{b}),  \\
            \RZ{2} & = (\RW{b},\RW{c}),\\
            \RZ{3} & = (\RW{a},\RW{c}),
        \end{aligned}  \nonumber
    \end{equation}
where $\RW{j}$ is an independent uniformly distributed random bit. It can be shown that $\RACO(V) = \frac{3}{2}$, $\RRACO(V) = \Set{(\frac{1}{2},\frac{1}{2},\frac{1}{2})}$ and all the components in $\RZ{V}$ are mutually independent. In this case, $\Set{\Set{1},\Set{2},\Set{3}}$ and $\Set{\Set{1,2,3}}$ are the minimizers of $\min_{\Pat \in \Pi(V)} \sum_{C \in \Pat} \FuHash{\alpha}(C)$ and $\Set{\Set{1},\Set{2},\Set{3}}$ is the only maximizer of $\RACO(V) = \max_{\Pat \in \Pi'(V)} \sum_{C \in \Pat} \frac{H(V) - H(C)}{|\Pat|-1}$. In this system, there is no complimentary user subset for SO, i.e., the optimal rate vector $(\frac{1}{2},\frac{1}{2},\frac{1}{2})$ cannot be implemented in a way such that the local omniscience in some non-singleton proper subset of $V$ is attained before the global omniscience in $V$.

Consider the system where $V=\{1,2,3\}$ and each user observes respectively
    \begin{equation}
        \begin{aligned}
            \RZ{1} & = (\RW{a}),  \\
            \RZ{2} & = (\RW{b}),\\
            \RZ{3} & = (\RW{c}),
        \end{aligned}  \nonumber
    \end{equation}
where $\RW{j}$ is an independent uniformly distributed random bit. It can be shown that $\RACO(V) = 3$, $\RRACO(V) = \Set{(1,1,1)}$ and all the components in $\RZ{V}$ are mutually independent. In this case, all partitions $\Pat \in \Pi(V)$ are the minimizers of $\min_{\Pat \in \Pi(V)} \sum_{C \in \Pat} \FuHash{\alpha}(C)$ and all partitions in $\Pat \in \Pi'(V)$ are the maximizers of $\RACO(V) = \max_{\Pat \in \Pi'(V)} \sum_{C \in \Pat} \frac{H(V) - H(C)}{|\Pat|-1}$. In this system, all non-singleton subset $X \subset V$ are complimentary.
\end{example}

Based on the proof of Theorem~\ref{theo:algo:CompSetSOExt} in Appendix~\ref{app:theo:algo:CompSetSOExt}, when the CompSetSO algorithm does not return a complimentary user set for an asymptotic and non-asymptotic models when $\alpha = \sum_{i \in V} \frac{H(V) - H(\Set{i})}{|V| - 1}$ and $\alpha = \big\lceil \sum_{i \in V} \frac{H(V) - H(\Set{i})}{|V| - 1} \big\rceil$, respectively, we have $\alpha = \RACO(V)$ and $\alpha= \RNCO(V)$, which necessarily means that the omniscience cannot be attained in a successive manner.