text:
\section{Introduction}
A breakthrough in our understanding of high-energy nuclear collisions
is the recognition~\cite{Alver:2006wh,Alver:2010gr}
that quantum fluctuations in the wavefunctions of
projectile and target, followed by hydrodynamic expansion, result in
unique long-range azimuthal correlations between outgoing particles.
The importance of these fluctuations was pointed out in the context of
detailed analyses of elliptic flow in nucleus-nucleus
collisions~\cite{Alver:2006wh,Miller:2003kd}.
It was later realized that fluctuations produce triangular
flow~\cite{Alver:2010gr}, which has subsequently been measured in nucleus-nucleus
collisions at RHIC~\cite{Adare:2011tg,Adamczyk:2013waa} and LHC~\cite{ALICE:2011ab,Chatrchyan:2012wg,ATLAS:2012at}.
Recently, fluctuations were predicted to generate significant
anisotropic flow in proton-nucleus
collisions~\cite{Bozek:2011if}, which quantitatively
explains~\cite{Bozek:2012gr}  the long-range correlations
observed by LHC experiments~\cite{CMS:2012qk,Abelev:2012ola,Aad:2012gla}.

Recently, the ATLAS and CMS experiments reported the observation of a
nonzero  4-particle cumulant of azimuthal correlations, dubbed
$v_2\{4\}$, in proton-nucleus
collisions~\cite{Aad:2013fja,Chatrchyan:2013nka}.
The occurrence of a large $v_2\{4\}$ in proton-nucleus collisions is
not fully understood, even though it is borne out by hydrodynamic
calculations with fluctuating initial conditions~\cite{Bozek:2013uha}.
Such higher-order cumulants were originally
introduced~\cite{Borghini:2000sa,Borghini:2001vi}
to measure elliptic flow in the reaction plane of non-central
nucleus-nucleus collisions, and isolate it from other, ``nonflow''
correlations.
It turns out that the simplest fluctuations one can think of, namely,
Gaussian fluctuations, do not contribute to
$v_2\{4\}$~\cite{Voloshin:2007pc}.
Since flow in proton-nucleus collisions is thought to originate from
fluctuations in the initial geometry,
one naively expects $v_2\{4\}\sim 0$, {\it even if\/}
there is collective flow in the system.

In this paper, we argue that the values observed for $v_2\{4\}$ in p-Pb collisions
are naturally explained by non-Gaussian fluctuations, which are expected for small systems.
Our explanation differs from that recently put forward by
Bzdak et al.~\cite{Bzdak:2013rya} that it is due to symmetry breaking
(see Eq.~(\ref{bg}) and discussion below).
As Bzdak et al., we assume that anisotropic flow $v_n$ scales like the
corresponding initial-state anisotropy $\varepsilon_n$ on an
event-by-event basis. This is known to be a very good approximation in
ideal~\cite{Qiu:2011iv} and viscous~\cite{Niemi:2012aj} hydrodynamics.
Thus flow fluctuations directly reflect $\varepsilon_n$ fluctuations.
Now, $\varepsilon_n$ is bounded by unity by definition.
On the other hand, Gaussian fluctuations are not bounded, which is the
reason why they fail to model small systems.
We propose a simple alternative to the Gaussian parameterization
which naturally satisfies the constraint $\varepsilon_n<1$. We show
that it provides an excellent fit to all Monte-Carlo
calculations.

\section{Distribution of the initial anisotropy}
In each event, the anisotropy in harmonic $n$ is defined (for $n=2,3$) by~\cite{Teaney:2010vd}
\begin{eqnarray}
\label{defepsilon}
\varepsilon_{n,x} &\equiv& -\frac{\int r^n \cos(n\phi)\rho(r,\phi)rdrd\phi}
{\int r^n\rho(r,\phi)rdrd\phi}\cr
\varepsilon_{n,y} &\equiv& -\frac{\int r^n \sin(n\phi)\rho(r,\phi)rdrd\phi}
{\int r^n\rho(r,\phi)rdrd\phi},
\end{eqnarray}
where $\rho(r,\phi)$ is the initial transverse density profile near midrapidity
 in a centered polar coordinate system.

\begin{figure}
 \includegraphics[width=\linewidth]{his_pA}
 \caption{ (Color online)
Histogram of the distribution of
$\varepsilon_2$ obtained in a Monte-Carlo Glauber simulation of a p-Pb
collision at LHC, and fits using
Eqs.~(\ref{gaussian})-(\ref{powerlaw}). }
\label{fig:histogram}
\end{figure}

Fig.~\ref{fig:histogram} displays the histogram of the distribution of
$\varepsilon_2$ in a p-Pb collision at 5.02~TeV obtained in a Monte-Carlo Glauber
calculation~\cite{Miller:2007ri}.
We use the PHOBOS implementation~\cite{Alver:2008aq} with a Gaussian
wounding profile~\cite{Alvioli:2009ab,Rybczynski:2011wv}.
We assume that the initial density $\rho(r,\phi)$ is a sum of
Gaussians of width $\sigma_0=0.4$~fm, centered around each participant
nucleon with a normalization that
fluctuates~\cite{Rybczynski:2013yba}.
These fluctuations, which increase anisotropies~\cite{Dumitru:2012yr},
are modeled as in Ref.~\cite{Bzdak:2013rya}.
We have selected events with number of participants $14\le N\le 16$,
corresponding to typical values in a central p-Pb collision.

We now compare different parameterizations of this distribution, which
we use to fit our numerical results.
The first is an isotropic two-dimensional Gaussian (we drop the
subscript $n$ for simplicity):
\begin{equation}
\label{gaussian}
P(\varepsilon)=\frac{2\varepsilon}{\sigma^2}\exp\left(-\frac{\varepsilon^2}{\sigma^2}\right),
\end{equation}
where $\varepsilon\equiv\sqrt{\varepsilon_x^2+\varepsilon_y^2}$ and
the distribution is normalized:
$\int_0^{\infty}P(\varepsilon)d\varepsilon=1$.
This form is motivated by the central limit theorem,
assuming that the eccentricity solely originates from event-by-event
fluctuations, and neglecting fluctuations in the denominator.
Note that this distribution does not strictly satisfy  the constraint
$\varepsilon<1$, which follows from the definition
(\ref{defepsilon}).
When fitting our Monte-Carlo results, we have therefore multiplied
Eq.~(\ref{gaussian})
by a constant to ensure normalization between 0 and 1.
The rms $\varepsilon$ has been fitted to that of the Monte-Carlo
simulation.
Fig.~\ref{fig:histogram} shows that Eq.~(\ref{gaussian}) gives a
reasonable approximation to our Monte-Carlo results, but not a good
fit.

Bzdak et al.~\cite{Bzdak:2013rya} have proposed to replace Eq.~(\ref{gaussian})
by a ``Bessel-Gaussian'':
\begin{equation}
\label{bg}
P(\varepsilon)=\frac{2\varepsilon}{\sigma^2}
I_0\left(\frac{2\varepsilon\bar\varepsilon}{\sigma^2}\right)
\exp\left(-\frac{\varepsilon^2+\bar\varepsilon^2}{\sigma^2}\right).
\end{equation}
This parameterization introduces an additional free parameter $\bar\varepsilon$,
corresponding to the mean eccentricity in the reaction plane in nucleus-nucleus
collisions~\cite{Voloshin:2007pc}.
It reduces to (\ref{gaussian}) if  $\bar\varepsilon=0$.
A nonzero value of $\bar\varepsilon$ is however difficult to justify for a
symmetric system in which anisotropies are solely created by fluctuations.
In Fig.~\ref{fig:histogram}, $\bar\varepsilon$ and $\sigma$ have been
chosen so that the first even moments $\langle\varepsilon^2\rangle$
and $\langle\varepsilon^4\rangle$ match exactly the Monte-Carlo
results, as suggested in~\cite{Bzdak:2013rya}.
The quality of the fit is not much improved compared to the Gaussian
distribution, even though there is an additional free parameter.
Note that the Bessel-Gaussian, like the
Gaussian, does not take into account the constraint $\varepsilon<1$.

We now introduce the one-parameter power law distribution:
\begin{equation}
\label{powerlaw}
P(\varepsilon)=2\alpha\varepsilon(1-\varepsilon^2)^{\alpha-1},
\end{equation}
where $\alpha>0$.
Eq.~(\ref{powerlaw}) reduces to Eq.~(\ref{gaussian}) for
$\alpha\gg 1$, with $\sigma^2\equiv 1/\alpha$.
The main advantage of Eq.~(\ref{powerlaw}) over previous
parameterizations is that the support of  $P(\varepsilon)$ is the
unit disc: it satisfies for all $\alpha>0$ the
normalization $\int_0^1 P(\varepsilon)d\varepsilon=1$.
In the limit $\alpha\to 0^+$,   $P(\varepsilon)\simeq\delta(\varepsilon-1)$.

Eq.~(\ref{powerlaw}) is the {\it
  exact\/}~\cite{Ollitrault:1992bk}\footnote{See Eq.~(3.10)
  of~\cite{Ollitrault:1992bk}. What is derived there is the
  distribution of anisotropy in momentum space, but the algebra is
  identical for the distribution of eccentricity.} distribution of
$\varepsilon_2$ for $N$ identical pointlike sources with a
2-dimensional isotropic Gaussian  distribution, with $\alpha=(N-1)/2$,
if one ignores the recentering correction.
In a more realistic situation, Eq.~(\ref{powerlaw}) is no longer
exact. We adjust $\alpha$ to match the rms $\varepsilon$ from the
Monte-Carlo calculation.
Fig.~\ref{fig:histogram} shows that Eq.~(\ref{powerlaw}) (with
$\alpha\simeq 5.64$) agrees much
better with Monte-Carlo results than Gaussian and
Bessel-Gaussian distributions.

\section{Cumulants}
Cumulants of the distribution of $\varepsilon$ are derived from a
generating function, which is the logarithm of the two-dimensional
Fourier transform of the distribution of $(\varepsilon_x,\varepsilon_y)$:
\begin{equation}
\label{defg}
G(k_x,k_y)\equiv \ln \langle \exp(ik_x \varepsilon_x+ik_y\varepsilon_y)\rangle,
\end{equation}
where angular brackets denote an expectation value over the ensemble of events.
If the system has azimuthal symmetry, by integrating over the relative
azimuthal angle of ${\bf k}$ and $\boldsymbol{\varepsilon}$, one obtains
\begin{equation}
\label{defgbessel}
G(k)=\ln \langle J_0(k\varepsilon)\rangle,
\end{equation}
where $k\equiv\sqrt{k_x^2+k_y^2}$ and
$\varepsilon\equiv\sqrt{\varepsilon_x^2+\varepsilon_y^2}$.
The cumulant to a given order $n$, $\varepsilon\{n\}$, is obtained by
expanding Eq.~(\ref{defgbessel}) to order $k^n$,
and identifying with the expansion of $\ln J_0(k\varepsilon\{n\})$ to
the same order.
This uniquely defines $\varepsilon\{n\}$ for all {\it even\/} $n$.
One thus obtains~\cite{Miller:2003kd}
$\varepsilon\{2\}^2=\langle\varepsilon^2\rangle$,
$\varepsilon\{4\}^4=2\langle\varepsilon^2\rangle^2-\langle\varepsilon^4\rangle$.
Expressions of $\varepsilon\{6\}$ and $\varepsilon\{8\}$ are given
in~\cite{Bzdak:2013rya}.

\begin{table}[ht]
\caption{Values of the first eccentricity cumulants
for the Gaussian (\ref{gaussian}), Bessel-Gaussian (\ref{bg}) and
power law
 (\ref{powerlaw}) distributions. }
\label{table:models}
\smallskip
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
&Gauss&BG&Power\\
\hline
$\varepsilon\{2\}$&$\sigma$&$\sqrt{\sigma^2+\bar\varepsilon^2}$&$\frac{\displaystyle
  1}{\displaystyle\sqrt{1+\alpha}}$\\
\hline
$\varepsilon\{4\}$&0&$\bar\varepsilon$&
$\left[\frac{\displaystyle 2}{\displaystyle (1+\alpha)^2(2+\alpha)}\right]^{1/4}$\\
\hline
$\varepsilon\{6\}$&0&$\bar\varepsilon$&
$\left[\frac{\displaystyle 6}{\displaystyle (1+\alpha)^3(2+\alpha)(3+\alpha)}\right]^{1/6}$\\
\hline
$\varepsilon\{8\}$&0&$\bar\varepsilon$&
$\left[\frac{\displaystyle 48\left(1+\frac{5\alpha}{11}\right)}
{\displaystyle (1+\alpha)^4(2+\alpha)^2(3+\alpha)(4+\alpha)}\right]^{1/8}$\\
\hline
\end{tabular}
\end{center}
\end{table}
Expressions of the first four cumulants are listed in
Table~\ref{table:models}.
For the power law distribution (\ref{powerlaw}), these results are
obtained by expanding the generating function (\ref{defgbessel}):
\begin{eqnarray}
\label{gpowerlaw}
G(k)&=&\ln\left[\int_0^1
J_0(k\varepsilon)P(\varepsilon)d\varepsilon\right]\cr
&=&\ln\left[\frac{2^{\alpha}\alpha!}{k^{\alpha}}J_{\alpha}(k)\right].
\end{eqnarray}
General results have been obtained previously
in the case of $N$ pointlike sources and in the large $N$ limit for
$\varepsilon_2\{2\}$~\cite{Bhalerao:2006tp} and
$\varepsilon_2\{4\}$~\cite{Alver:2008zza}.
Our results derived from Eq.~(\ref{powerlaw}) are exact for a Gaussian distribution
of sources and therefore agree with these general results for $N\gg 1$.
Similar results have also been derived for $\varepsilon_3\{2\}$ and
$\varepsilon_3\{4\}$~\cite{Bhalerao:2011bp}, but not for
cumulants of order 6 or higher.

\begin{figure}
 \includegraphics[width=\linewidth]{e2}
 \caption{ (Color online) Cumulants of the eccentricity distribution
   as a function of the number of participants $N$
   for the power law distribution (\ref{powerlaw}), where we have
set   $\alpha=(N-2)/2$.
}
\label{fig:gaussian}
\end{figure}

Fig.~\ref{fig:gaussian} displays the cumulants $\varepsilon\{2\}$ to
$\varepsilon\{8\}$ as a function of $N$, as
predicted by Eq.~(\ref{powerlaw}) for pointlike
sources.\footnote{Here, we
  assume that the recentering correction
  effectively reduces   by one unit the number of independent
  sources.
  We thus replace $N$ by $N-1$ in the exact result of
  Ref.~\cite{Ollitrault:1992bk}.}
These results are similar to those obtained in full Monte-Carlo Glauber
calculations~\cite{Bzdak:2013rya}.
In the limit  $N\gg 1$, the power law distribution yields
$\varepsilon\{k\}\propto N^{(1-k)/k}$.
It thus predicts
a strong ordering $\varepsilon\{8\}\ll \varepsilon\{6\}\ll \varepsilon\{4\}\ll \varepsilon\{2\}\ll 1$,
unlike the Bessel-Gaussian which predicts
$\varepsilon\{4\}=\varepsilon\{6\}=\varepsilon\{8\}$.
For fixed $N$, however, the cumulant expansion
quickly converges, as illustrated in Fig.~\ref{fig:gaussian}.
In practice, for typical values of $N$ in p-Pb collisions, one
observes $\varepsilon\{4\}\simeq \varepsilon\{6\}\simeq \varepsilon\{8\}$, in
agreement with numerical findings of Bzdak et
al.~\cite{Bzdak:2013rya}.
This rapid convergence can be traced back to the fact that the
generating function $G(k)$ in Eq.~(\ref{gpowerlaw}) has a singularity
at the first zero of $J_\alpha(k)$, denoted by $j_{\alpha 1}$.
This causes the cumulant expansion to quickly converge to the
value~\cite{Bhalerao:2003xf}
\begin{equation}
\label{lyz}
\varepsilon\{\infty\}=\frac{j_{01}}{j_{\alpha 1}}.
\end{equation}
This asymptotic limit is also plotted in Fig.~\ref{fig:gaussian}. It
is hardly distinguishable from $\varepsilon\{6\}$ and $\varepsilon\{8\}$
for these values of $N$.

\begin{figure}
 \includegraphics[width=\linewidth]{universal}
 \caption{ (Color online) $\varepsilon\{4\}$ versus
   $\varepsilon\{2\}$. The dashed line in both panels is
   Eq.~(\ref{eps4vseps2}). Left: p-Pb collisions. ``Full'' refers to
   Gaussian sources associated with each participant, and
   fluctuations in the weights of each source. ``Pointlike'' refers to
   pointlike identical sources. DIPSY results for p-p collisions are replotted
   from~\cite{Avsar:2010rf}. Right: Pb-Pb collisions. The dotted
   line is $\varepsilon\{4\}=\varepsilon\{2\}$, corresponding to a
   nonzero mean eccentricity, and negligible fluctuations.}
\label{fig:e2e4}
\end{figure}

\section{Testing universality}
The power law distribution (\ref{powerlaw}) predicts the following
parameter-free relation between the first two cumulants:
\begin{equation}
\label{eps4vseps2}
\varepsilon\{4\}=\varepsilon\{2\}^{3/2}\left(\frac{2}{1+\varepsilon\{2\}^2}\right)^{1/4}.
\end{equation}
This relation can be used to test the universality of the
distribution (\ref{powerlaw}).
For p-Pb collisions at 5.02~TeV, we run two different types of
Monte-Carlo Glauber
calculations: a full Monte-Carlo identical to that of
Fig.~\ref{fig:histogram}, and a second one where fluctuations and
smearing are switched off (identical pointlike sources).
We calculate $\varepsilon_2$ and $\varepsilon_3$ for each event.
Events are then binned according to the number of participants $N$,
mimicking a centrality selection.
For p-p collisions at 7~TeV, we use published
results~\cite{Avsar:2010rf} obtained with the event
generator DIPSY~\cite{Flensburg:2011wx}, which are binned according
to multiplicity.
Results are shown in Fig.~\ref{fig:e2e4} (left).
Each symbol of a given type corresponds to a different bin.
All Monte-Carlo results are in very good agreement with Eq.~(\ref{eps4vseps2}).
A closer look at the results show that the ``full'' Monte-Carlo
Glauber calculations are above the line by $\sim 0.015$ (for both
$\varepsilon_2$ and $\varepsilon_3$), the
``pointlike'' results for $\varepsilon_3$ by $\sim 0.005$, and the
``pointlike'' results for $\varepsilon_2$ (where our result is exact, up to the
recentering correction) by $\sim 0.002$.
DIPSY results are above the line by $\sim 0.01$.

For Pb-Pb collisions at 2.76~TeV (Fig.~\ref{fig:e2e4} right), we use the
results obtained in Ref.~\cite{Bhalerao:2011yg} using the Monte-Carlo
Glauber~\cite{Alver:2008aq} and Monte-Carlo KLN~\cite{Drescher:2007ax}
models. These results are in 5\% centrality bins.
For $\varepsilon_3$, both models are in very good agreement with
Eq.~(\ref{eps4vseps2}) (within $0.01$ or so). Note that Pb-Pb
collisions probe this relation closer to the origin, in the large $N$
limit where more general results are
available~\cite{Bhalerao:2011bp}. These general  results predict
$\varepsilon\{4\}\propto\varepsilon\{2\}^{3/2}$ for $N\to\infty$, but
with a proportionality constant that depends on the density
profile. Our results show that it is in practice very close to
the value predicted by Eq.~(\ref{eps4vseps2}), namely, $2^{1/4}$.

Monte-Carlo results for $\varepsilon_2$ in Pb-Pb differ from
Eq.~(\ref{eps4vseps2}). This
is expected, since $\varepsilon_2$ in mid-central Pb-Pb collisions is
mostly driven by the almond shape of the overlap area between colliding
nuclei~\cite{Ollitrault:1992bk}, not by fluctuations.
In the limiting case where fluctuations are negligible,
$\varepsilon_2\{4\}=\varepsilon_2\{2\}$.
Our results show that fluctuations dominate only for the most central
and most peripheral bins.

We conclude that the power law distribution (\ref{powerlaw}) is
a very good approximation to the distribution of fluctuation-driven
eccentricities, irrespective of the details of the model.
This could be checked explicitly with other initial-state
models~\cite{Dumitru:2012yr,Schenke:2012hg}.

\begin{figure}
 \includegraphics[width=\linewidth]{predictions}
 \caption{ (Color online) Predictions of the model for ratios of
   higher order cumulants and $\varepsilon\{2\}$ as a function of the measured
   $v\{4\}/v\{2\}$. Typical values for $v_3$ in
   Pb-Pb~\cite{ALICE:2011ab,Chatrchyan:2013kba}  and $v_2$ in p-Pb
   collisions~\cite{Chatrchyan:2013nka} are indicated by arrows. }
\label{fig:predictions}
\end{figure}
\section{Applications}
We now discuss applications of our result.
The distribution of $\varepsilon_n$ is completely
determined by the parameter $\alpha$ in Eq.~(\ref{powerlaw}). This
parameter can be obtained directly from experimental data.
Assuming that anisotropic flow is proportional to eccentricity in the
corresponding harmonic, $v_n\propto\varepsilon_n$, which is proven to
be a very good approximation for $n=2,3$~\cite{Niemi:2012aj}, one
obtains
\begin{equation}
\frac{v\{4\}}{v\{2\}}=\frac{\varepsilon\{4\}}{\varepsilon\{2\}}=
\left(\frac{2}{2+\alpha}\right)^{1/4}.
\end{equation}
The first equality has already been checked against Monte-Carlo models
and experimental data~\cite{Bhalerao:2011ry,Chatrchyan:2013kba}. The
second equality directly relates the parameter $\alpha$ in
Eq.~(\ref{powerlaw}) to the measured ratio $v\{4\}/v\{2\}$.

This in turn gives a prediction for ratios of higher-order flow
cumulants, which scale like the corresponding ratios of eccentricity
cumulants. These predictions are displayed in
Fig.~\ref{fig:predictions}.
One can also directly obtain the rms eccentricity $\varepsilon\{2\}$,
which is a property of the initial state.

The ratio  $v_3\{4\}/v_3\{2\}$ in Pb-Pb is close to
0.5 in mid-central collisions~\cite{ALICE:2011ab,Chatrchyan:2013kba}.
We thus predict $v_3\{6\}/v_3\{4\}\simeq 0.84$ and
$v_3\{8\}/v_3\{6\}\simeq 0.94$ in the same centrality.
We also obtain $\varepsilon_3\{2\}\simeq
0.17$, which is a typical prediction from Monte-Carlo models
in the 10\%-20\% or 20\%-30\% centrality
range~\cite{Retinskaya:2013gca}.

Similarly, the ratio $v_2\{4\}/v_2\{2\}\sim 0.7$ measured in p-Pb
collisions~\cite{Aad:2013fja,Chatrchyan:2013nka} implies
$v_2\{6\}/v_2\{4\}\simeq 0.93$ and
$v_2\{8\}/v_2\{6\}\simeq 0.98$, that is, almost degenerate
higher-order cumulants.
We obtain $\varepsilon_2\{2\}\simeq 0.37$, in agreement with
Monte-Carlo Glauber models~\cite{Bzdak:2013rya}.

\section{Conclusions}
We have proposed a new parameterization of the distribution
of the initial anisotropy $\varepsilon_n$ in proton-proton,
proton-nucleus and nucleus-nucleus
collisions which, unlike previous parameterizations, takes into
account the condition $\varepsilon_n<1$.
This new parameterization is found in good agreement with results
of Monte-Carlo simulations when
$\varepsilon_n$ is created by fluctuations of the initial geometry.
Our results explain the observation, in these Monte-Carlo models, that
cumulants of the distribution of $\varepsilon_n$ quickly converge as
the order increases. This is because the Fourier transform of the
distribution of $\varepsilon_n$ has a zero at
a finite value of the conjugate variable $k$. This, in turn, is a
consequence of the fact that the probability distribution of
$\varepsilon_n$ has compact support (that is, $\varepsilon_n<1$).

The consequence of this universality is that while the rms
$\varepsilon_n$ is strongly model-dependent~\cite{Retinskaya:2013gca},
the probability distribution of $\varepsilon_n$ is fully determined
once the rms value is known --- in particular, the magnitudes of
higher-order cumulants such as $\varepsilon_n\{4\}$.
Assuming that anisotropic flow $v_n$ is proportional to
$\varepsilon_n$ in every event, we have predicted the values of
$v_3\{6\}$ and $v_3\{8\}$ in Pb-Pb collisions, and the values of
$v_2\{6\}$ and $v_2\{8\}$ in p-Pb collisions.

If future experimental data confirm our prediction, these results will
strongly support the picture that the long-range correlations observed
in proton-nucleus and nucleus-nucleus collisions are due to
anisotropic flow, which is itself proportional to the anisotropy in
the initial state. This picture, furthermore, will be confirmed
irrespective of the details of the initial-state model.

\begin{acknowledgments}
JYO thanks Art Poskanzer for pointing out, back in 2009, that Bessel-Gaussian fits
to Monte-Carlo Glauber calculations fail because
they miss the constraint $\varepsilon_2<1$,
Larry McLerran for discussing Ref.~\cite{Bzdak:2013rya} prior to
publication, Christoffer Flensburg for sending DIPSY results,
Ante Bilandzic and Wojciech Broniowski for useful discussions,
and Jean-Paul Blaizot and Raju Venugopalan for comments on the manuscript.
We thank the Yukawa Institute for Theoretical Physics, Kyoto University.
Discussions during the YITP workshop YITP-T-13-05 on ``New Frontiers in QCD'' were useful to complete this work.
LY is funded  by the European Research Council under the
Advanced Investigator Grant ERC-AD-267258.
\end{acknowledgments}


text:
\section{\textcolor{black}{\normalsize Introduction}}

An intermediate state is a quantum state which reduces to two or
more distinguishably different states (normally, distinguishable in
terms of photon number distribution) in different limits. In 1985,
such a state was first time introduced by Stoler \emph{et al.}
{[}\ref{the:D-Stoler,-B}]. To be precise, they introduced Binomial
state (BS) as a state which is intermediate between the most
nonclassical number state $|n\rangle$ and the most classical
coherent state $|\alpha\rangle$. They defined BS as \begin{equation}
\begin{array}{lr}
|p,M\rangle=\sum_{n=0}^{M} & B_{n}^{M}\end{array}|n\rangle=\sum_{n=0}^{M}\sqrt{^{M}C_{n}p^{n}(1-p)^{M-n}}|n\rangle\,\,\,\:0\leq p\leq1.\label{eq:binomial1}\end{equation}
This state%
\footnote{The state is named as binomial state because the photon
number distribution associated with this state
$\left(i.e.\,|B_{n}^{M}|^{2}\right)$ is
simply a binomial distribution.%
} is called intermediate state as it reduces to number state in the
limit $p\rightarrow0$ and $p\rightarrow1$ (as
$|0,M\rangle=|0\rangle$ and $|1,M\rangle=|M\rangle$) and in the
limit of $M\rightarrow\infty,\, p\rightarrow1$, where $\alpha$ is a
real constant, it reduces to a coherent state with real amplitude.
Since the introduction of BS as an intermediate state, it was always
been of interest to quantum optics, nonlinear optics, atomic physics
and molecular physics community. Consequently, different properties
of binomial states have been studied
{[}\ref{the:VIdiella-Barraco}-\ref{the:OEBS}]. In these studies it
has been observed that the nonclassical phenomena (such as,
antibunching, squeezing and higher order squeezing) can be seen in
BS. This trend of search for nonclassicality in Binomial state,
continued in nineties. In one hand, several versions of generalized
BS have been proposed
{[}\ref{the:hong-chenfu-genralized-BS}-\ref{the:Hong-yi-Fan-generalized-BS}]
and in the other hand, people went beyond binomial states and
proposed several other form of intermediate states (such as,
\textcolor{black}{excited binomial state {[}\ref{the:EBS}], odd
excited binomial state {[}\ref{the:OEBS}], hypergeometric state
{[}\ref{the:hypergeomtric}], negative hypergeometric state
{[}\ref{the:negativehyper geometric}], reciprocal binomial state
{[}\ref{the:reciprocalBS}], shadowed state
{[}\ref{the:CTLee(shaddowed)}], shadowed negative binomial state
{[}\ref{the:ranjani and ctlee}] and photon added coherent state
{[}\ref{the:agarwal-photonadded}] etc.).} The studies in the
nineties were mainly limited to theoretical predictions but the
recent developments in the experimental techniques made it possible
to verify some of those theoretical predictions. For example, we can
note that, as early as in 1991 Agarwal and Tara
{[}\ref{the:agarwal-photonadded}] had introduced photon added
coherent state as \begin{equation} |\alpha,m\rangle=\frac{a^{\dagger
m}|\alpha\rangle}{\langle a^{m}a^{\dagger
m}\rangle},\label{eq:photonadded1}\end{equation} (where $m$ is an
integer and $|\alpha\rangle$ is coherent state) but the experimental
generation of the state has happened only in recent past when
Zavatta, Viciani and Bellini {[}\ref{the:photonadded-experiment}]
succeed to produce it in 2004. It is easy to observe that
(\ref{eq:photonadded1}) represents an intermediate state, since it
reduces to coherent state in the limit $m\rightarrow0$ and to number
state in the limit $\alpha\rightarrow0$. \textcolor{black}{This
state can be viewed as a coherent state in which additional $m$
photon are added. The photon number distribution of all the above
mentioned states are different but all these states belong to a
common family of states called intermediate state.} It has also been
found that most of these intermediate states show antibunching,
squeezing, higher order squeezing, subpoissonian photon statistics
etc. but higher order antibunching has been reported only in
shadowed negative binomial state {[}\ref{the:ranjani and ctlee}].
Inspired by these observations, many schemes to generate
intermediate states have been proposed in recent past
{[}\ref{the:photonadded-experiment}-\ref{the:R LOFranco2}].

\textcolor{black}{The reason behind the study of nonclassical properties
of intermediate states lies in the fact that the most of the interesting
recent developments in quantum optics have arisen through the nonclassical
properties of the radiation field only. For example, antibunching
and squeezing, which do not have any classical analogue {[}\ref{elements of quantum optics}-\ref{hbt}],
have extensively been studied in last thirty years. But the majority
of these studies are focused on lowest order nonclassical effects.
Higher order extensions of these nonclassical states have only been
introduced in recent past {[}\ref{hong}-\ref{lee2}].  Among these
higher order nonclassical effects, higher order squeezing has already
been studied in detail {[}\ref{hong}, \ref{hong2}, \ref{hillery},
\ref{giri}] but the higher order antibunching (HOA) is not yet studied
rigourously. }

\textcolor{black}{The idea of HOA was introduced by Lee in a
pioneering paper {[}\ref{lee1}] in 1990, since then it has been
predicted in two photon coherent state {[}\ref{lee1}], shadowed
negative binomial state {[}\ref{the:ranjani and ctlee}], trio
coherent state {[}\ref{ba an}] and in the interaction of intense
laser beam with an inversion symmetric third order nonlinear medium
{[}\ref{the:martin}]. From the fact that in first 15 years after its
introduction, HOA was reported only in some particular cases, HOA
appeared to be a very rare phenomenon. But recently we have shown
that the HOA is not a rare phenomenon {[}\ref{the:pathak-jphysB}]
and it can be seen in simple optical processes} like six wave mixing
process, four wave mixing process and second harmonic
generation\textcolor{black}{. To establish that further, here we
have shown the existence of HOA in different intermediate states,
namely, binomial state, reciprocal binomial state, photon added
coherent state, hypergeometric state, Roy-Roy generalized binomial
state and negative binomial state. }

\textcolor{black}{The present work is motivated by the recent
experimental observation of intermediate state
{[}\ref{the:photonadded-experiment}], theoretical observation of
possibility of observing HOA in some simple optical systems
{[}\ref{the:pathak-jphysB}] and the fact that the intermediate
states, which frequently show different kind of nonclassicality,
form a big family of quantum state. But till now HOA has been
predicted only in one member  (Shadowed negative binomial state) of
such a big family of quantum states {[}\ref{the:ranjani and ctlee}].
Motivated by these facts the present work aims to study the
possibility of HOA in all the popularly known intermediate states.
The theoretical predictions of the present study can be
experimentally verified with the help of various intermediate state
generation schemes and homodyne experiment, since the criteria for
HOA appears in terms of factorial moment, which can be measured by
using homodyne photon counting experiments
{[}\ref{the:Bachor-H-A,}-\ref{the:singh2}]. In the next section we
have briefly described the criteria of HOA. In section 3 we have
shown that the HOA of any arbitrary order can be seen in BS. In
section 4, Roy-Roy generalized binomial state
{[}\ref{the:broy&proy-generalized-bs}] is studied and existence of
HOA is predicted. Calculational details and methodology have been
shown only in algebraically simple cases which are described in
section 3 and 4. Section 5 is divided in several subsections and we
have followed the same procedure and have studied the possibilities
of observing HOA in different intermediate states namely, reciprocal
binomial state, negative binomial state, hypergeometric states and
photon added coherent state. One subsection is dedicated for the
discussion of one particular intermediate state. Finally section 6
is dedicated to conclusions. }

\section{\textcolor{black}{\normalsize Criteria of HOA}}
\textcolor{black}{The criterion of HOA is expressed in terms of
higher order factorial moments of number operator. There exist
several criterion for the same which are essentially equivalent.
Initially, using the negativity of P function {[}\ref{elements of
quantum optics}], Lee introduced the criterion for HOA as }

\textcolor{black}{\begin{equation} R(l,m)=\frac{\left\langle
N_{x}^{(l+1)}\right\rangle \left\langle N_{x}^{(m-1)}\right\rangle
}{\left\langle N_{x}^{(l)}\right\rangle \left\langle
N_{x}^{(m)}\right\rangle }-1<0,\label{eq:ho3}\end{equation} where
$N$ is the usual number operator, $\left\langle N^{(i)}\right\rangle
=\left\langle N(N-1)...(N-i+1)\right\rangle $ is the $ith$ factorial
moment of number operator, $\left\langle \right\rangle $ denotes the
quantum average, $l$ and $m$ are integers satisfying the conditions
$1\leq m\leq l$ and the subscript $x$ denotes a particular mode. Ba
An {[}\ref{ba an}] choose $m=1$ and reduced the criterion of $l$th
order antibunching to \begin{equation} A_{x,l}=\frac{\left\langle
N_{x}^{(l+1)}\right\rangle }{\left\langle N_{x}^{(l)}\right\rangle
\left\langle N_{x}\right\rangle }-1<0\label{eq:bhuta1}\end{equation}
or, \begin{equation} \left\langle N_{x}^{(l+1)}\right\rangle
<\left\langle N_{x}^{(l)}\right\rangle \left\langle
N_{x}\right\rangle .\label{eq:ba an (cond)}\end{equation}
Physically, a state which is antibunched in $l$th order has to be
antibunched in $(l-1)th$ order. Therefore, we can further simplify
(\ref{eq:ba an (cond)}) as \begin{equation} \left\langle
N_{x}^{(l+1)}\right\rangle <\left\langle N_{x}^{(l)}\right\rangle
\left\langle N_{x}\right\rangle <\left\langle
N_{x}^{(l-1)}\right\rangle \left\langle N_{x}\right\rangle
^{2}<\left\langle N_{x}^{(l-2)}\right\rangle \left\langle
N_{x}\right\rangle ^{3}<...<\left\langle N_{x}\right\rangle
^{l+1}\label{eq:ineq}\end{equation} and obtain the condition for
$l-th$ order antibunching as \begin{equation} d(l)=\left\langle
N_{x}^{(l+1)}\right\rangle -\left\langle N_{x}\right\rangle
^{l+1}<0.\label{eq:ho21}\end{equation} This simplified criterion
(\ref{eq:ho21}) coincides exactly with the physical criterion of HOA
introduced by Pathak and Garica {[}\ref{the:martin}] and the
criterion of Erenso, Vyas and Singh {[}\ref{the:Erenso}], recently
Vogel has reported a class of nonclassicality conditions based on
higher order factorial moments {[}\ref{the:vogel}]. All these
criteria essentially lead to same kind of nonclassicality which
belong to the class of strong nonclassicality according to the
classification scheme of Arvind} \textcolor{black}{\emph{et al}}
\textcolor{black}{{[}\ref{the:Arvind}]. Here we can note that
$d(l)=0$ and $d(l)>0$ corresponds to higher order coherence and
higher order bunching (many photon bunching) respectively.}
Actually, $\left\langle a^{\dagger l}a^{l}\right\rangle
=\left\langle N^{(l)}\right\rangle $ is a measure of the probability
of observing $l$ photons of the same mode at a particular point in
space time coordinate. Therefore the physical meaning of
inequalities (\ref{eq:ineq}) is that the probability of detection of
single photon pulse is greater than that of a two photon in a bunch
and that is greater than the probability of detection of three
photon in a bunch and so on. This is exactly the characteristic that
is required in a probabilistic single photon source used in quantum
cryptography. In other words all the probabilistic single photon
sources used in quantum cryptography should satisfy the criteria
(\ref{eq:ho21}) of HOA {[}\ref{the:IJP}].

\section{Binomial State}

Binomial state is originally defined as (\ref{eq:binomial1}), from
which it is straight forward to see that \begin{equation}
\begin{array}{lcl}
a|p,M\rangle & = & \sum_{n=0}^{M}\left\{ \frac{M!}{(n-1)!(M-n)!}p^{n}(1-p)^{M-n}\right\} ^{\frac{1}{2}}|n-1>\\
 & = & \sum_{l=0}^{M-1}\left\{ \frac{M(M-1)!}{l!(M-1-l)!}p^{l+1}(1-p)^{M-1-l}\right\} ^{\frac{1}{2}}|l>\,\,\,\,\,\,\,\,(assuming\, n-1=l)\\
 & = & [Mp]^{\frac{1}{2}}\sum_{l=0}^{M-1}\left\{ \frac{(M-1)!}{l!(M-1-l)!}p^{l}(1-p)^{M-1-l}\right\} ^{\frac{1}{2}}|l>\\
 & = & [Mp]^{\frac{1}{2}}|p,M-1\rangle.\end{array}\label{eq:eigen1}\end{equation}
 Similarly, we can write, \begin{equation}
\begin{array}{lcl}
a^{2}|p,M\rangle & = & [M(M-1)p^{2}]^{\frac{1}{2}}|p,M-2\rangle\\
a^{3}|p,M\rangle & = & [M(M-1)(M-2)p^{3}]^{\frac{1}{2}}|p,M-3\rangle\\
\vdots & \vdots & \vdots\\
a^{l}|p,M\rangle & = & [M(M-1)....(M-l+1)p^{l}]^{\frac{1}{2}}|p,M-l\rangle\\
 & = & \left[\frac{M!}{(M-l)!}p^{l}\right]^{\frac{1}{2}}|M-l,p\rangle.\end{array}\label{eq:eigen2}\end{equation}
Therefore, \begin{equation} \langle M,p|a^{\dagger l}=\langle
M-l,p|\left[\frac{M!}{(M-l)!}p^{l}\right]^{\frac{1}{2}}\label{eq:eigen3}\end{equation}
and consequently, \begin{equation} \langle
M,p|n^{(l)}|p,M\rangle=\langle M,p|a^{\dagger
l}a^{l}|p,M\rangle=\left[\frac{M!}{(M-l)!}p^{l}\right].\label{eq:expectation1}\end{equation}
Now substituting (\ref{eq:expectation1}) in equation (\ref{eq:ho21})
we obtain the condition for $lth$ order antibunching as
\begin{equation}
d(l)=\left[\frac{M!}{(M-l-1)!}p^{l+1}\right]-[Mp]^{l+1}<0\label{eq:dl-binomial}\end{equation}
or, \begin{equation}
(M-1)(M-2)....(M-l)<M^{l}\label{eq:hoasatisfied}\end{equation} which
is always satisfied for any $M>l$ and both $M$ and $l$ are positive
(since every term in left is $<M$). As $M$ is the number of photons
present in the field and $d(l)$ is a measure of correlation among
$(l+1)$ photons, therefore $M\geq(l+1)$ or $M>l$. Consequently, a
binomial state always shows HOA and the highest possible order of
antibunching that can be seen in a binomial state is equal to $M-1$,
where $M$ is the number of photon present in the field. From
(\ref{eq:dl-binomial}) it is straight forward to see that the number
state is always higher order antibunched and in the other extreme
limit (when $p\rightarrow1,\, M\rightarrow\infty$ and the BS reduces
to coherent state) $d(l)=0,$ which is consistent with the physical
expectation.

\section{Generalized Binomial State}

We have already mentioned that there are different form of
generalized binomial states
{[}\ref{the:hong-chenfu-genralized-BS}-\ref{the:Hong-yi-Fan-generalized-BS}].
For the present study we have chosen generalized binomial state
introduced by Roy and Roy {[}\ref{the:broy&proy-generalized-bs}].
They have introduced the generalized binomial state (GBS) as
\begin{equation}
|N,\alpha,\beta\rangle=\sum_{n=o}^{N}\sqrt{\omega(n,N,\alpha,\beta)}|n\rangle\label{eq:gen-bino1}\end{equation}
where, \begin{equation}
\omega(n,N,\alpha,\beta)=\frac{N!}{(\alpha+\beta+2)_{N}}\frac{(\alpha+1)_{n}(\beta+1)_{N-n}}{n!(N-n)!}\label{eq:gen-bino2}\end{equation}
 with $\alpha,\beta>-1$, $n=0,1,....,N$, and \begin{equation}
\begin{array}{lcr}
(a)_{0}=1 & \,\,\,\,\,\,\, &
(a)_{n}=a(a+1)....(a+n-1)\end{array}.\label{eq:gen-bino3}\end{equation}
This intermediate state reduces to vacuum state, number state,
coherent state, binomial state and negative binomial state in
different limits of $\alpha$, $\beta$ and $N$. In order to obtain an
analytic expression of $d(l)$ for this particular generalized
binomial state we need to prove following useful identity:

\textbf{Identity1:} \begin{equation}
a(a+1)_{n}=(a)_{n+1}\label{eq:identity 1}\end{equation}

Proof: Using (\ref{eq:gen-bino3}) we can write \[
(a+1)_{n}=(a+1)...(a+n)=\frac{a(a+1)...(a+n)}{a}=\frac{(a)_{n+1}}{a}.\]
Therefore, \[ a(a+1)_{n}=(a)_{n+1}.\]

Now it is easy to see that the above identity (\ref{eq:identity 1})
yields the following useful relations: \begin{equation}
(\alpha+1)_{l+1}=(\alpha+1)(\alpha+2)_{l}\label{eq:relation1}\end{equation}
 and \begin{equation}
(\alpha+\beta+2)_{N}=(\alpha+\beta+2)(\alpha+\beta+3)_{N-1}=(\alpha+\beta+2)(\alpha+2+\beta+1)_{N-1}.\label{eq:relation2}\end{equation}
Using (\ref{eq:gen-bino1}) and (\ref{eq:gen-bino2}) we can obtain
\begin{equation}
\begin{array}{lcl}
a|N,\alpha,\beta\rangle & = & \sum_{n=0}^{N}\left\{ \frac{N!}{(\alpha+\beta+2)_{N}}\frac{(\alpha+1)_{n}}{(n-1)!}\frac{(\beta+1)_{N-n}}{(N-n)!}\right\} ^{\frac{1}{2}}|n-1\rangle\\
 & = & \sum_{l=0}^{N-1}\left\{ \frac{N(N-1)!}{(\alpha+\beta+2)_{N}}\frac{(\alpha+1)_{l+1}}{l!}\frac{(\beta+1)_{N-1-l}}{(N-1-l)!}\right\} ^{\frac{1}{2}}|l\rangle\end{array},\label{eq:relation 3}\end{equation}
where $n=l-1$ has been used. Now we can apply (\ref{eq:relation1})
and (\ref{eq:relation2}) on (\ref{eq:relation 3}) to obtain \begin{equation}
\begin{array}{lcl}
a|N,\alpha,\beta\rangle & = & \left\{ \frac{N(\alpha+1)}{(\alpha+\beta+2)}\right\} ^{\frac{1}{2}}\sum_{l=0}^{N-1}\left\{ \frac{(N-1)!(\alpha+2)_{l}(\beta+1)_{N-1-l}}{(\alpha+2+\beta+1)_{N-1}l!(N-1-l)!}\right\} ^{\frac{1}{2}}|l\rangle\\
 & = & \left\{ \frac{N(\alpha+1)}{(\alpha+\beta+2)}\right\} ^{\frac{1}{2}}\sum_{n=0}^{N-1}\sqrt{\omega(n,N-1,\alpha+1,\beta)}|n\rangle,\end{array}\label{eq:gb1}\end{equation}
 where dummy variable $l$ is replaced by $n.$ Therefore, \[
\begin{array}{lcl}
\langle N,\alpha,\beta|a^{\dagger}a|N,\alpha,\beta\rangle & = & \frac{N(\alpha+1)}{(\alpha+\beta+2)}\\
\langle N,\alpha,\beta|a^{\dagger2}a^{2}|N,\alpha,\beta\rangle & = & \frac{N(N-1)(\alpha+1)(\alpha+2)}{(\alpha+\beta+2)(\alpha+\beta+3)}\\
\vdots & \vdots & \vdots\\
\langle N,\alpha,\beta|a^{\dagger l}a^{l}|N,\alpha,\beta\rangle & = & \frac{\left[N(N-1).....(N-l+1)\right]\left[(\alpha+1)(\alpha+2).....(\alpha+l)\right]}{(\alpha+\beta+2)(\alpha+\beta+3).....(\alpha+\beta+l+1)}\\
 & = & \frac{N!(\alpha+l)!(\alpha+\beta+1)!}{(N-l)!\alpha!(\alpha+\beta+l+1)!}\end{array}\]
and \begin{equation}
\begin{array}{lcl}
d_{GBS}(l) & = & \frac{N!(\alpha+l+1)!(\alpha+\beta+1)!}{(N-l-1)!\alpha!(\alpha+\beta+l+2)!}-\left\{ \frac{N(\alpha+1)}{(\alpha+\beta+2)}\right\} ^{l+1}\\
 & = & \frac{\left[N(N-1).....(N-l)\right]\left[(\alpha+1)(\alpha+2).....(\alpha+l+1)\right]}{(\alpha+\beta+2)(\alpha+\beta+3).....(\alpha+\beta+l+2)}-\left\{ \frac{N(\alpha+1)}{(\alpha+\beta+2)}\right\} ^{l+1}\end{array}\label{eq:d(l)broy}\end{equation}
The physical condition $N\geq l+1$ ensures that all the terms in
$d(l)$ are positive. The expression of $d(l)$ is quite complex and
it depends on various parameters (e.g. $\alpha,\,\beta$ and $N$).
Fig 1 shows that for particular values of these parameters HOA is
possible. As it is expected from the earlier works on the properties
of HOA {[}\ref{the:martin}], the depth of nonclassicality is more in
case of $d_{GBS}(9)$ than in $d_{GBS}(8)$. This is consistent with
earlier observation. A systematic study reveals that the probability
of observing HOA increases with the increase of $\alpha$ but it
decreases (i.e the probability of higher order bunching increases)
with the increase of $\beta$. This can be seen clearly in Fig2 and
Fig3. Further it is observed (from Fig1 and Fig2) that for lower
values of $\alpha$ the probability of bunching increases with the
increase of $N$ but for a comparatively large values of $\alpha$
(larger compare to $\beta$) the probability of HOA increases with
the increase in $N$ {[}see Fig 2] but the situation is just opposite
in the case of $\beta$ (as it is seen from Fig3).

\begin{figure}[h]
\centering \scalebox{0.6}{\includegraphics{gbs3_a2_b1.jpg}}
\caption{Higher order antibunching can be seen in Generalized
binomial state. Existence of 8th and 9th order antibunching (for
$\alpha=2$ and $\beta=1$)and variation of depth of nonclassicality
with $N$ has been shown.}\label{Fig1}
\end{figure}
\begin{figure}[h]
\centering \scalebox{0.6}{\includegraphics{GBS_1.jpg}} \label{Fig2}
\caption{Variation of $d_{GBS}(2)$ with $\alpha$ and $N$ for
$\beta$=1.} \label{fig2}
\end{figure}
\begin{figure}[h]
\centering \scalebox{0.6}{\includegraphics{GBS_2.jpg}}
\caption{Variation of $d_{GBS}(2)$ with $\beta$ and $N$ for
$\alpha$=10.} \label{Fig3}
\end{figure}
While studying different limiting cases of Roy and Roy generalized
binomial state, we have observed that binomial state and number
state always show HOA and $d(l)=0$ for coherent state. This is
consistent with the physical expectation and the conclusions of the
last section.

\section{Other Intermediate States}

As it is mentioned in the earlier sections, there exist several
different intermediate states. For the systematic study of
possibility of observing HOA in intermediate states, we have studied
all the well known intermediate states. Since the procedure followed
for the study of different states are similar, mathematical detail
has not been shown in the subsections below. But from the expression
of $d(l)$ and the corresponding plots it would be easy to see that
the HOA can be observed in all the intermediate states studied
below.

\subsection{Reciprocal binomial state}

Reciprocal binomial state (RBS) can be defined as
{[}\ref{the:reciprocalBS}]
\begin{equation}
|\phi\rangle=\frac{1}{N}\begin{array}{c}
{\scriptstyle N}\\
\sum\\
{\scriptscriptstyle {\scriptstyle k=0}}\end{array}\left(\begin{array}{c}
N\\
K\end{array}\right)^{-1/2}e^{ik(\theta-\pi/2)}|k\rangle\label{eq:rbs1}\end{equation}
where \emph{N} is a normalization constant. Using antinormal
ordering and procedure adapted in section 3, we can
obtain;\begin{equation}
\begin{array}{lcl}
d_{RBS}(l) & = & \begin{array}{c}
{\scriptstyle l+1}\\
\sum\\
i{\scriptstyle =0}\end{array}(-1)^{i}\frac{(l+1)!^{2}}{(l+1-i)!^{2}i!}\frac{(N+(l+1-i\,)!}{N!}-N^{l+1}\\
 & = & \frac{\pi\csc\left(\pi\left(l+N\right)\right)(l+1)!^{2}Gamma(l+1-N)}{M!\left(Gamma(2+l)Gamma(-N)\right)^{2}}-N^{l+1}\end{array}\label{eq:rbs2}\end{equation}
Where $Gamma$ denotes the Gamma function. The possibility of
observing HOA in reciprocal binomial state can be clearly seen from
the Figure 4. But it is interesting to note that the nature of
singularity and zeroes present in the simplified expression of
$d_{RBS}(l)$ as expressed in the last line of (\ref{eq:rbs2}) can
provide us some important information. For example, the underlying
mathematical structure of the criterion of HOA and that of
reciprocal binomial state demands that $l$ and $N$ be integers but
if both of them are integer then $d_{RBS}(l)$ has a singularity as
the $\csc[\pi(l+N)]$ term present in the numerator blows up. But
this local singularity can be circumvented by assuming $l\rightarrow
integer$ and $N\rightarrow integer$. In this situation (i.e. when
$l$ and $N$ tends to integer value) $\langle N^{(l)}\rangle$ is
finite and consequently $d(l)$ is also finite. This is the reason
that the singular nature of the simplified expression of
$d_{RBS}(l)$ is not reflected in the Figure 4. In the analysis of
the $d_{RBS}(l)$ it is also interesting to observe that
$Gamma(-N)=\infty$ for $N=integer$, and in an approximated situation
when $l\rightarrow integer$ and $M=integer$, the $\csc[\pi(l+N)]$
term in the numerator is no more singular and as a result$\langle
N^{(l)}\rangle=0$ and $d_{RBS}=-N^{(l+1)}$. In this situation one
can observe HOA for arbitrarily large values of $l$ and $N$. Thus
physically, it is expected that in reciprocal binomial state higher
order antibunching of any arbitrary order will be seen and HOA will
not be destroyed with the increase of $N$, as it happens (for some
particular values of $\alpha$ and $\beta$) in the case of
generalized binomial state (see Figure 1).
\begin{figure}[h]
\centering
 \scalebox{0.6}{\includegraphics{RBS-1.jpg}} \caption{Variation of $d_{RBS}(8)$ and $d_{RBS}(9)$ with
photon number $N$ } \label{Fig5}
\end{figure}

\subsection{Negative Binomial State}

Following Barnett {[}\ref{the:Barnett}] we can define Negative
Binomial state (NBS) as

\begin{equation}
|\eta,M\rangle=\sum_{n=M}^{\infty}C_{n}(\eta,M)|n\rangle\label{eq:NBS1}\end{equation} where
$C_{n}(\eta,M)=\left[\left(\begin{array}{c}
n\\
M\end{array}\right)\eta^{M+1}(1-\eta)^{n-M}\right]^{1/2}$, $0\leq\eta\leq1$ and M is a non-negative integer.
This intermediate state interpolates between number state and geometric state. Following the mathematical
techniques adopted in the earlier sections we obtain
\begin{equation}
d_{NBS}(l)=\eta^{-l}\left(\frac{\left(l+M+1\right)!\,_{2}F_{1}\left(-l-1,-l-1;-l-M-1;\eta\right)}{M!}-\frac{(M+1)^{l+1}}{\eta}\right)\label{eq:NBS2}\end{equation}
where $\,_{2}F_{1}(a,b;c;z)$ is a conventional hypergeometric function. Variation of $d(l)$ with various
parameters such as $\eta,\, l$ and $M$ have been studied and are shown in Figure 5- Figure6. From these figures
one can observe that the state is not always antibunched, rather the plot of $d_{NBS}(8)$ has a very sharp rise
near $\eta\approx.15$ and $M\approx10$. From Figure 5 we can observe that the possibility of higher order
bunching get destroyed with the increase of $M$ or $\eta$. The second feature has clearly been shown in Figure
6, where we can see that $d_{NBS}(8)$ become positive for a small region for $M=10$ but for a bigger region
($\eta\approx.3$ to $\eta=1$) it remains negative and thus shows the existence of HOA. The broad features remain
same for the other orders (other values of $l$) of antibunching.
\begin{figure}[h]
\centering
 \scalebox{0.5}{\includegraphics{NBS_3D.jpg}} \caption{Variation of $d_{NBS}(8)$  with $\eta$ and $M$
} \label{Fig6}
\end{figure}

\begin{figure*}[h]
\centering
\begin{tabular}{cc}
 \scalebox{0.5}{\includegraphics{NBS_M10_1.jpg}}& \scalebox{0.5}{\includegraphics{NBS_M10_2.jpg}}
\end{tabular}\caption{Variation of
$d_{NBS}(8)$  with $\eta$ for $M=10$ }\label{Fig7}
\end{figure*}

In the limit $N\rightarrow0$ the negative binomial state reduces to
geometric state (GS). In this limit $d_{NBS}(l)$ reduces to
\begin{equation}
d_{GS}(l)=\frac{1}{\eta^{l+1}}\left((1-\eta)^{l+1}\eta(l+1)!-1\right).\label{eq:geometric
1}\end{equation} It is interesting to observe that the above
expression has a singularity at $\eta=0$ and
$d_{GS}(l)\rightarrow-\infty$ in the limit $\eta\rightarrow 0$.
Consequently one negative values of $d_{GS}$ can be seen at very
very small values $\eta$ but this is not the signatures of HOA,
rather this is the signature of the existence of a strong
singularity in the neighborhood. This can further justified by the
fact that for any finite value of $l$ there does not exist any real
root, (whose value is close to zero or which is negligibly small
compared to 1) of $d_{GS}(l)=0$. Thus there is no oscillation
between bunching and antibunching. We further observe that
$d_{GS}(l)\rightarrow-1$ in the limit $\eta\rightarrow1$ and for
$l\geq3$ there exists only one physically acceptable real root of
$d_{GS}(l)=0$. By physically acceptable real root we mean that it
lies in $[0,1]$. Before this value of $\eta$ (or before the
physically acceptable real root) the state shows higher order super
poissonian photon statistics but immediately after the root it
becomes negative and thus shows HOA or higher order subpoissonian
photon statistics. As we increase $l$ the real root shifts in right
side of the real axis (i.e towards $\eta=1$). It can be clearly seen
in the Figure 7. From this figure it can be easily seen that it
satisfies all the physical properties of HOA derived in
{[}\ref{the:martin}].
\begin{figure*}[h]
\centering
 \scalebox{0.7}{\includegraphics{NBS-GS-2.jpg}} \caption{Variation of $d_{GS}(8)$, $d_{GS}(9)$ and
$d_{GS}(10)$ with respect to $\eta$. \label{Fig4}}
\end{figure*}

\subsection{Photon added coherent state}

Photon added coherent state (\ref{eq:photonadded1}) or PACS, which
was introduced by Agarwal and Tara {[}\ref{the:agarwal-photonadded}]
can be defined as
\begin{equation}
|\alpha,m\rangle=\frac{\exp(-\frac{|\alpha|^{2}}{2})}{\sqrt{L_{m}(-|\alpha|^{2})m!}}\sum_{n=0}^{\infty}\alpha^{n}\,\frac{\sqrt{(m+n)!}}{n!}|n+m\rangle\label{eq:photonadded1n}\end{equation}
where $L_{m}(x)$ is Lauguere polynomial of $mth$ order. Rigorous operator algebra yields

\begin{equation}
\begin{array}{lcl}
d_{PACS}(l) & = & \frac{\exp(-\alpha^{2})\alpha^{2l+2}((l+m+1)!)^{2}\,_{P}F_{Q}\left(\left\{ 1,2+l+m,2+l+m\right\} ;\left\{ 2+l,2+l,m+1\right\} ;\alpha^{2}\right)}{(m!(l+1)!)^{2}\,_{1}F_{1}\left(-m;1;-\alpha^{2}\right)}\\
 & - & \left(\frac{\exp(-\alpha^{2})\left(-m+m\,_{1}F_{1}\left(1+m;1;\alpha^{2}\right)+(1+m)\alpha^{2}\,_{1}F_{1}\left(2+m;2;\alpha^{2}\right)\right)}{\,_{1}F_{1}\left(-m;1;-\alpha^{2}\right)}\right)^{l+1},\end{array}\label{eq:photonadded2}\end{equation}
where, $\,_{P}F_{Q}$ is the generalized Hypergeometric function. The
analytic expression for $d_{PACS}(l)$ is quite complicated and it is
difficult to conclude anything regarding its photon statistics
directly from (\ref{eq:photonadded2}) but we have investigated the
variation of $d_{PACS}$ with $\alpha,\, l$ and $m$ and could not
find any region which does not show HOA. Therefore, HOA can be seen
in this particular intermediate state. This fact is manifested in
Fig 8 and Fig 9. From these figures it is easy to observe that depth
of nonclassicality increases monotonically with the increase of $m$
and $l$. The variation of depth of nonclassicality with $\alpha$ has
a deep for a small value of $\alpha$ (see Fig 8. and Fig. 9).
Although $d_{PACS}$ is always negative, initially its magnitude
($d_{PACS}(l)$ without the negative sign) increases, then decreases
and then become a monotonically increasing function. Actually for
the smaller values of $\alpha$, an effective contribution from the
combination of all the hypergeometric functions appears and
dominates but as soon as $\alpha$ increases a bit, the
$\exp(-\alpha^{2})$ term starts dominating and as a consequence
depth of nonclassicality increases monotonically. Here we would also
like to note that in contrast to the photon added coherent state,
$d(l)$ is always positive for the analogous state
$|\alpha,-m\rangle$, introduced by Sivakumar
{[}\ref{the:Shiavakumar}]. Thus $|\alpha,-m\rangle$ always shows
higher order superpoisonian photon statistics. Further, we would
like to note that photon added coherent state which, is intermediate
between coherent and fock state has already been experimentally
generated in 2004 {[}\ref{the:photonadded-experiment}]. Therefore,
it is technically feasible to observe higher order antibunching for
an intermediate state.

\begin{figure}[h]
\centering
 \scalebox{0.6}{\includegraphics{PACS-3D.jpg}} \caption{Variation of $d_{PACS}(4)$  with
 $\alpha$ and $m$} \label{Fig8}
\end{figure}
\begin{figure}[h]
\centering
 \scalebox{0.6}{\includegraphics{PACS.jpg}} \caption{Variation of $ 10 d_{PACS}(3)$ and $d_{PACS}(4)$ with
 $\alpha$ for $m=15$, the solid line denotes $d_{PACS}(3)$ and the dashed line denotes $d_{PACS}(4)$ to keep
 show the variation in the same scale $d_{PACS}(3)$ is multiplied by 10. The plot shows that depth of
 nonclassicality of $d_{PACS}(4)$ is always greater than that of $d_{PACS}(3)$ which is consistent with the
 properties of HOA} \label{Fig9}
\end{figure}
\subsection{Hypergeometric state}

Following {[}\ref{the:hypergeomtric}] hypergeometric state (HS) can
be defined as \begin{equation}
|L,M,\eta\rangle=\sum_{n=0}^{M}H_{n}^{M}(\eta,L)|n\rangle\label{eq:hyper1}\end{equation}
where the probability $\eta$ is a real parameter satisfying
$0<\eta<1.$ $L$ is a real number satisfying \begin{equation} L\geq
max\left\{ M\eta^{-1},\, M(1-\eta)^{-1}\right\}
,\label{eq:hyper2}\end{equation}
\begin{equation}
H_{n}^{M}(\eta,L)=\left[\left(\begin{array}{c}
L\eta\\
n\end{array}\right)\left(\begin{array}{c}
L\left(1-\eta\right)\\
M-n\end{array}\right)\right]^{\frac{1}{2}}\left(\begin{array}{c}
L\\
M\end{array}\right)^{-\frac{1}{2}},\label{eq:hyper3}\end{equation} and \begin{equation} \left(\begin{array}{c}
\alpha\\
n\end{array}\right)=\frac{\alpha(\alpha-1)...(\alpha-n+1)}{n!},\,\,\,\,\,\,\,\left(\begin{array}{c}
\alpha\\
0\end{array}\right)\equiv1.\label{eq:hyper4}\end{equation} Here
$\alpha$ is not necessarily an integer. Using the techniques adopted
in the earlier sections and a bit of operator algebra we can obtain
a closed form analytic expression for $d(l)$ as \begin{equation}
d_{HS}(l)=-\left(M\eta\right)^{l+1}+\frac{(L-l-1)!M!(L\eta)!}{L!(M-l-1)!(L\eta-1-l)!}\label{eq:hyper5}\end{equation}
From Fig. 10 it is clear that HOA can be observed in hypergeometric
state. It is also observed that the depth of nonclassicality
increases with the increase in $\eta$ and $M$. Hypergeometric state
reduces to binomial state, coherent state, number state and vacuum
state in different limits of $M,$ $L$ and $\eta$. It has been
verified that if we impose those limits on $d_{HS}$ then we obtain
corresponding photon statistics.

\begin{figure}[h]
\centering
 \scalebox{0.6}{\includegraphics{Hyper_3D.jpg}} \caption{Variation of $d_{HS}(8)$  with
 $\eta$ and $M$, when lowest allowed values of $L$ have been chosen at every point. } \label{Fig10}
\end{figure}

\section{Conclusions}

In essence all the intermediate states studied in the present paper
\textcolor{black}{show HOA i.e. higher order subpoissonian photon
statistics. But it neither mean that all the intermediate states are
higher order antibunched (for example $|\alpha,-m\rangle$ is always
higher order bunched) nor an intermediate state which shows HOA will
show it for all possible values of the control parameters (for
example negative binomial state and generalized binomial state shows
both higher order bunching and higher order antibunching for various
parametric values). Thus we can conclude that, as far as HOA is
concerned there does not exist any common characteristics among the
different intermediate states but most of them show HOA. Further, we
have seen from (Figure1-Figure10) that the depth of nonclassicality
of a higher order antibunched state varies with different control
parameters (e.g. $\alpha$, $N$, $m$ etc.). These parameters
represent some physical quantity and their value may be controlled
and consequently by controlling these parameters we can control the
depth of nonclassicality.
 Photon statistics (factorial moment) of an intermediate state
can be obtained experimentally by using homodyne detection (photon
counting) technique. These facts along with the recent success in
experimental production of intermediate state open up the
possibility of experimental observation of HOA in intermediate
state. Thus the present work strongly establishes the fact that HOA
is not a rare phenomenon. }

\textcolor{black}{The prescription followed in the present work is easy and straight forward and it can be used
to study the possibilities of observing higher order antibunching in other intermediate states (such as negative
hyper geometric state, excited binomial state and odd excited binomial state) and other physical systems. Thus
it opens up the possibility of studying higher order nonclassical effects from a new perspective. This is also
important from the application point of view because any probabilistic single photon source used for quantum
cryptography has to satisfy the condition for higher order antibunching. }

\textbf{Acknowledgement}: AP thanks to DST, India for partial financial
support through the project grant SR\textbackslash{}FTP\textbackslash{}PS-13\textbackslash{}2004.

redpajama_set_name: RedPajamaArXiv