The reggeometric pomeron and exclusive production

of $J/\psi$ and $\psi(2S)$ in ultraperipheral collisions at the LHC

László Jenkovszky a⋆, Vladyslav Libova⋄, and Magno V. T. Machadob♠

a Bogolyubov ITP, National Academy of Sciences of Ukraine, Kiev 03143 Ukraine

b HEP Phenomenology Group, CEP 91501-970, Porto Alegre, RS, Brazil

###### Abstract

By using a Regge-pole model for vector meson production (VMP), successfully
describing the HERA data, we analyse the correlation between VMP cross
sections in photon-induced reactions at HERA and those in ultra-peripheral
collisions at the Large Hadron Collider (LHC). The rapidity distributions of
proton-proton collisions at 13 TeV and lead-lead collisions at 2.76 and 5.02
TeV are investigated. The transverse momentum distribution in nuclear coherent
vector meson production is also addressed. Predictions for future experiments
on production of $J/\psi$ and $\psi(2S)$ are presented.

$\begin{array}[]{ll}{}^{\star}\mbox{{\it e-mail
<EMAIL_ADDRESS>{}^{\diamond}\mbox{{\it e-mail
<EMAIL_ADDRESS>{}^{\spadesuit}\mbox{{\it
e-mail<EMAIL_ADDRESS>\end{array}$

##  1 Introduction

Recently, new data on $J/\psi$ production in ultra-peripheral lead-lead (PbPb)
collisions became public [1, 2, 3, 4, 5]. Interestingly, ALICE collaboration
has measured for the first time the $|t|$-distribution in the coherent nuclear
process $\gamma+A\rightarrow J/\psi+A$ [4]. The photo-production at low
transverse momentum has been studied in peripheral PbPb collisions by LHCb
collaboration at 5 TeV [5]. Vector meson production in ultra-peripheral
collisions (UPC) were studied in a number of papers (see comprehensive reviews
in Refs. [6, 7, 8, 9]). Ultra-peripheral vector meson production was studied,
both experimentally and theoretically in proton-proton ($pp$), proton-nucleus
($pA$) and nucleus-nucleus ($AA$) collisions at the LHC energies and below.
Production of various vector mesons – from light to heavy ones, such as
$\psi(nS)$ and $\Upsilon(nS)$ states were studied. In these studies the
incoming energy, rapidity and momentum transfer of the produced system were
scrutinized. These experimental studies have triggered a large number of
theoretical investigations. For quarkonium production basically three classes
of models exist: those based on (generalized) vector dominance model
(GVDM/VDM) [10, 11, 12], perturbative QCD/BFKL calculations [13, 14, 15, 16]
and color dipole framework [17, 18, 19, 20, 21, 22, 23, 24]. The first one is
quite appealing due to its simplicity and deeply connection with reggeon
phenomenology and reggeon field theory [25].

In this study the VMP at the CERN LHC will be scrutinized, specializing for
the $J/\psi$ and $\psi(2S)$ photo-production in $pp$ and PbPb collisions. It
updates our previous analysis [26, 27, 28] which was restricted to $pp$
collisions at 7 TeV. In Ref. [26] a single-component Reggeometric Pomeron
model has been considered and the focus was on the energy dependence of cross
section. The inclusion of a two-pomeron (soft + hard) model and the soft
dipole pomeron model was performed in Refs. [27, 28]. Here, the analysis will
be extended to include $pp$ collisions at 13 TeV. By using VDM approximation
and Glauber multiple scattering formalism, we analyse the production in PbPb
collisions at 2.76 and 5.04 TeV (predictions for XeXe collisions are also
provided). The transverse momentum distribution will be investigated in the
PbPb mode.

This paper is organized as follows. In Sec. 2 we shortly review the VMP in the
context of the single-component Reggeometric Pomeron model. In Subsection 2.1
the formalism to obtain the rapidity distribution in $pp$ collisions is shown.
The numerical results both for $pp$ and PbPb reactions are compared to
experimental measurements in Sec. 3 and discussion on the theoretical
uncertainties is done. In the last section we summarize the main achievements
and prospects for near future are presented.

##  2 Theoretical formalism

At DESY-HERA experiment, vector meson production (VMP) was investigated in
details both by the H1 and ZEUS collaborations. The same is true for the
deeply virtual Compton scattering (DVCS) process. In the kinematical regime
corresponding to diffractive scattering, those processes can be described by a
Pomeron exchange (see Fig. 1). Pomeron dominance is especially clean in
$J/\Psi$ production because, due to the Zweig (OZI) rule any exchange of
secondary trajectories is forbidden. In this case, one has an uncontaminated
Pomeron exchange which gives us the opportunity to study the nature of the
Pomeron. An important feature is the two-fold nature of the Pomeron, which
seems to be soft or hard depending on the virtuality, $Q^{2}$, of the incident
photon and/or the mass of the produced vector meson, $M_{V}$.

Figure 1: Diagrams of DVCS (a) and VMP (b) in $e^{\pm}p$ scattering; (c) DVCS
(VMP) amplitude in a Regge-factorized form (including representation for
particle vertices, $V_{1,2}$).

It is common in literature to assume that two Pomerons exist: a hard (or pQCD)
Pomeron, resulting from perturbative Quantum Chromodynamics (QCD)
calculations, and a soft one, somewhat misleadingly called non-perturbative.
We instead believe that there is only one Pomeron in Nature, but it has two
components, whose relative weight is regulated by relevant
$\widetilde{Q}^{2}$-dependent factors in front of the two, where the measure
of the hardness, $\widetilde{Q}^{2}=Q^{2}+M_{V}^{2}$, is the sum of the
squared photon virtuality $Q^{2}$ and the squared mass $M^{2}_{V}$ of the
produced vector meson.

A specific model realizing this idea was constructed and tested against the
experimental data few years ago (see Refs. [30, 31] and earlier references
therein). The relevant VMP scattering amplitude reads

$\displaystyle A(s,t,Q^{2},{M_{V}}^{2})$ $\displaystyle=$
$\displaystyle\widetilde{A_{s}}e^{-i\frac{\pi}{2}\alpha_{s}(t)}\left(\frac{s}{s_{0}}\right)^{\alpha_{s}(t)}e^{b_{s}t-n_{s}\ln{\left(1+\frac{\widetilde{Q^{2}}}{\widetilde{Q_{s}^{2}}}\right)}}$
(1) $\displaystyle+$
$\displaystyle\widetilde{A_{h}}e^{-i\frac{\pi}{2}\alpha_{h}(t)}\left(\frac{s}{s_{0}}\right)^{\alpha_{h}(t)}e^{b_{h}t-(n_{h}+1)\ln{\left(1+\frac{\widetilde{Q^{2}}}{\widetilde{Q_{h}^{2}}}\right)}+\ln{\left(\frac{\widetilde{Q^{2}}}{\widetilde{Q_{h}^{2}}}\right)}},$

where $\alpha_{s}(t)$ and $\alpha_{h}(t)$ are the soft and hard Pomeron
trajectories, $s=W^{2}$, $W$ being the energy of the VMP (the center-of-mass
energy of photon-nucleon system). Let us stress that the Pomeron is unique in
all reactions, but its components (and parameters) vary. Examples with
detailed fits can be found in the papers of Ref. [30, 31].

The integrated VMP cross section is given by

$\displaystyle\sigma_{\gamma p\rightarrow
Vp}(\widetilde{Q}^{2},W)=\int_{t_{\mathrm{m}}}^{t_{\mathrm{thr}}}\frac{d\sigma(\widetilde{Q}^{2},W,t)}{dt},$
(2)

where the upper limit is $t_{\mathrm{thr}}=0$ GeV2 and the lower limit is the
kinematical one, $t_{\mathrm{m}}=-s/2$. The total cross section can be simply
calculated, without integration for an exponential diffraction cone, according
to the formula

$\sigma_{el}(s)=\frac{1}{B(s)}\frac{d\sigma}{dt}\biggr{\rvert}_{t=0},$ (3)

where $B$ is the forward slope.

Figure 2: Feynman diagram of vector meson production in a hadronic collision.

For practical reasons, an effective Pomeron which can be soft or hard,
depending on kinematical domain or process defining its hardness. In the
single-component Reggeometric Pomeron the amplitude takes the simple form:

$\displaystyle
A(Q^{2},s,t)=\widetilde{H}(\widetilde{Q^{2}})e^{-\frac{i\pi\alpha(t)}{2}}\left(\frac{s}{s_{0}}\right)^{\alpha(t)}e^{2\left(\frac{a}{\widetilde{Q^{2}}}+\frac{b}{2m_{N}^{2}}\right)t},\quad\widetilde{H}(\widetilde{Q^{2}})\equiv\frac{\widetilde{A_{0}}}{\left(1+\frac{\widetilde{Q^{2}}}{{Q_{0}^{2}}}\right)^{n}},$
(4)

where the exponent in the exponential factor in Eq. (4) reflects the
geometrical nature of the model. The quantities $a/\widetilde{Q}^{2}$ and
$b/2m_{N}^{2}$ correspond to the effective sizes of upper and lower vertices
in Fig. 1-c. The corresponding elastic differential cross section is given by:

$\displaystyle\frac{d\sigma_{el}}{dt}=\frac{\pi}{s^{2}}|A(Q^{2},s,t)|^{2},\quad\frac{d\sigma_{el}}{dt}=\frac{A_{0}^{2}}{\left(1+\frac{\widetilde{Q^{2}}}{{Q_{0}^{2}}}\right)^{2n}}\left(\frac{s}{s_{0}}\right)^{2(\alpha(t)-1)}e^{4\left(\frac{a}{\widetilde{Q^{2}}}+\frac{b}{2m_{N}^{2}}\right)t}.$
(5)

Since our primary goal is the comparison between the energy (rapidity)
dependence of VMP at HERA and the LHC, we start with very simple ansätz for
the $\gamma p\rightarrow Vp$ cross section, postponing the use of the more
elaborated model given by Eq. (1) to a future study. Here, we follow the
studies on the single-component Reggeometric Pomeron model which has been
applied to deeply virtual Compton scattering (DVCS) in Ref. [29]. Apart from
$W$ and $t$, it contains also the dependence on virtuality, $Q^{2}$. The model
was fitted to the HERA data on DVCS, but it can be applied also to the VMP by
refitting its parameters. In this paper we use the first (simpler) version of
the Reggeometric model of VMP and DVCS, proposed in Refs. [30, 31].

The Reggeometric Pomeron model applies to photoproduction ($Q^{2}=0$) limit,
where $\tilde{Q}^{2}=M_{V}^{2}$. The total photoproduction cross section for
$J/\psi$ meson is given by [30, 31]:

$\displaystyle\sigma_{\gamma p\to J/\Psi p}(W,\widetilde{Q}^{2})$
$\displaystyle=$
$\displaystyle\frac{A_{0}^{2}}{\left(1+\widetilde{Q}^{2}/\widetilde{Q}^{2}_{0}\right)^{2n}}\frac{\left(\frac{W}{W_{0}}\right)^{4(\alpha_{0}-1)}}{B_{V}\left(W,\widetilde{Q}^{2}\right)},$
(6) $\displaystyle B_{V}\left(W,\widetilde{Q}^{2}\right)$ $\displaystyle=$
$\displaystyle
4\left[\alpha^{\prime}\ln(W/W_{0})+\left(\frac{a}{\widetilde{Q}^{2}}+\frac{b}{2m_{N}^{2}}\right)\right].$
(7)

The parameters, fitted [30, 31] to the $J/\psi$ photoproduction data, are
presented in Table 1. Note that compared to the original formula, $s$ was
replaced by $W^{2}$, since $W$ is used in this paper to denote the photon-
proton centre-of-mass energy (in contrast, $\sqrt{s}$ in this paper is the
proton-proton centre-of-mass energy). Our baseline model is defined by the set
of parameters fitted in Ref. [30]. In next subsection we summarize the main
expressions used to describe the rapidity distribution of exclusive meson
production in proton-proton collisions. The input is the production cross
section presented above, Eq. (6).

Table 1: Values of the parameters presented in Refs. [30, 31] fitted to data on VMP at HERA. Our baseline prediction is based in parameter from Ref. [30]. Ref. |  $A_{0}$ $\left[\frac{\sqrt{\text{nb}}}{\text{GeV}}\right]$ |  $\widetilde{Q^{2}_{0}}$ $\left[\text{GeV}^{2}\right]$ | $n$ | $\alpha_{0}$ |  $\alpha^{\prime}$ $\left[\text{GeV}^{-2}\right]$ | $a$ | $b$
---|---|---|---|---|---|---|---
[30] |  29.8 $\pm$ 2.8 |  2.1 $\pm$ 0.4 |  1.37 $\pm$ 0.14 |  1.20 $\pm$ 0.02 |  0.17 $\pm$ 0.05 |  1.01 $\pm$ 0.11 |  0.44 $\pm$ 0.08
[31] |  30$\pm$31 |  2.3$\pm$2.2 |  1.45$\pm$0.32 |  1.21$\pm$0.09 |  0.077$\pm$0.072 | 1.72 | 1.16

###  2.1 Rapidity distribution in proton-proton collisions

Vector meson production cross section in hadronic collisions can be written in
a factorized form in the equivalent photon approximation (EPA) [6]. The
distribution in rapidity $y$ of the production of a vector meson $V$ in the
reaction $h_{1}+h_{2}\rightarrow h_{1}+V+h_{2}$ (where $h_{i}$ stands for a
hadron) is calculated according to a standard prescription based on the
factorization of the photon flux and photon-proton cross section (see Fig. 2).

Generally speaking, the $\gamma p$ differential cross section depends on three
variables: the total energy of the $\gamma p$ system, $W$, the squared
momentum transfer at the proton vertex, $t$, and
$\widetilde{Q}^{2}=Q^{2}+M_{V}^{2}$, where, as previously said, $Q^{2}$ is the
photon virtuality and $M_{V}$ is the mass of the produced vector meson. Since,
by definition, in ultra-peripheral collisions we have $b>>R_{1}+R_{2}$, where
$b$ is the impact parameter, i.e. the closest distance between the centres of
the colliding particles/nuclei with radii $R_{i}$ (i = 1,2), photons are
nearly real, $Q^{2}\simeq 0$, and $M_{V}^{2}$ remains the only measure of
hardness. Notice that this might not be true for peripheral collisions, where
$b\sim R_{1}+R_{2}$, and in the Pomeron or Odderon exchange instead of the
photon.

Figure 3: The cross section for exclusive $J/\Psi$ photoproduction, $\gamma
p\to J/\Psi p$, as a function of the photon-proton c.m.s. energy compared with
different data sets. In addition to data considered in Refs. [26, 27, 28], we
added the latest ALICE [32, 33] and LHCb [34] data. Dotted line is the
prediction from the Reggeometric Pomeron model.

The differential cross section, $d\sigma/dy$, as function of the meson
rapidity $y$ reads

$\displaystyle\frac{d\sigma(h_{1}+h_{2}\rightarrow h_{1}+V+h_{2})}{dy}$
$\displaystyle=$ $\displaystyle
S_{\mathrm{gap}}(\omega_{+})\,\omega_{+}\frac{dN_{\gamma/h_{1}}(\omega_{+})}{d\omega}\sigma_{\gamma
h_{2}\rightarrow Vh_{2}}(\omega_{+})$ (8) $\displaystyle+$ $\displaystyle
S_{\mathrm{gap}}(\omega_{-})\,\omega_{-}\frac{dN_{\gamma/h_{2}}(\omega_{-})}{d\omega}\sigma_{\gamma
h_{1}\rightarrow Vh_{1}}(\omega_{-}),$

where $dN_{\gamma/h}/d\omega$ is the equivalent photon flux [6] for protons.
The quantity $S_{\mathrm{gap}}(\omega_{\pm})$ is the gap survival factor.
Here, we will use the following expression for the photon flux:

$\frac{dN_{\gamma/h}(\omega)}{d\omega}=\frac{\alpha_{em}}{2\pi\omega}[1+(1-\frac{2\omega}{\sqrt{s}})^{2}](\ln\Omega-\frac{11}{6}+\frac{3}{\Omega}-\frac{3}{2\Omega^{2}}+\frac{1}{3\Omega^{3}}),$
(9)

and $\sigma_{\gamma p\rightarrow Vp}(\omega)$ is the total cross section of
the vector meson photoproduction sub-process. Here $\omega$ is the photon
energy, $\omega=W^{2}_{\gamma p}/2\sqrt{s}$, where $\sqrt{s}$ denotes now the
c.m.s. energy of the proton-proton system;
$\omega_{min}=M_{V}^{2}/(4\gamma_{L}m_{p}),$ where
$\gamma_{L}=\sqrt{s}/(2m_{p})$ is the Lorentz factor (Lorentz boost of a
single beam). Furthermore, $\Omega=1+Q_{0}^{2}/Q_{min}^{2},$ with
$Q_{min}^{2}=(\omega/\gamma_{L})^{2}$ and $Q_{0}^{2}=0.71GeV^{2}$. The photon
energy can be written in terms of the meson rapidity as
$\omega=\frac{M_{v}}{2}e^{y}$.

In Refs. [26, 27, 28], we have addressed the rapidity distribution for
$J/\psi$ and $\psi(2S)$ in $pp$ collisions at the LHC for the energy of 7 TeV.
Comparison of different models for the scattering amplitude has been done,
where an overall gap survival factor $S_{\mathrm{gap}}=0.8$ was utilized. In
particular, analysis was done using the Reggeometric Pomeron model, a
soft+hard Pomeron and the soft dipole Pomeron model as well. Moreover, the
energy dependence of the extracted cross section using hadronic $pp$
scattering was scrutinized. For simplicity, we restrict ourselves to the
Reggeometric Pomeron model in present analysis. In order to update the studies
done in Refs. [26, 27, 28] in Fig. 3 is presented the cross section
$\sigma_{\gamma p\to J/\Psi p}(W,\widetilde{Q}^{2})$ as a function of $W$.
From experimental side we added recent data from ALICE [32, 33] and LHCb [34].
The result is compared to our prediction in Fig. 3 (dotted line). Also the
ZEUS and H1 data are shown.

In what follows we investigate the rapidity distribution for $J/\psi$ and
$\psi(2S)$ in proton-proton collisions at 13 TeV. Moreover, we address the
rapidity and transverse momentum distributions in PbPb ultraperipheral
collisions at 2.76 and 5.02 TeV. Nuclear effects are described by vector
dominance model and Glauber approach for multiple scattering.

##  3 Results and discussions

Figure 4: Rapidity distribution of exclusive $J/\psi$ production in $pp$
collisions at the LHC ($\sqrt{s_{pp}}=13$ TeV). Predictions of Reggeometric
Pomeron are represented by solid and dashed curves, respectively. Comparison
is done with the H1 power-like fit (dotted line) and STARLIGHT fit (dot-
dashed). Data from LHCb collaboration [34].

In this section the numerical results are presented and compared to recent
experimental measurements. Let us start with the rapidity distribution of
$J/\psi$ and $\psi(2S)$ production in proton-proton collisions at the LHC (see
Eq. (8)). In Fig. 4 the predictions are shown for $J/\psi$ production at 13
TeV against the LHCb data [34]. The results for the Reggeometric Pomeron model
are represented by the solid and dashed curves, respectively. The upper curve
(solid) is the prediction using the values for the parameters determined in
Ref. [31]. On the other hand, the lower curve (dashed) corresponds to input
presented in Ref. [30] (see Table 1). The $y$-dependent gap survival
probabilities, $S_{\mathrm{gap}}(\omega_{\pm})$, are taken from Ref. [35]. In
order to compare our prediction to other approaches in literature, it is shown
also the result using the H1 power-law fit [36] to exclusive photoproduction,
$\sigma(\gamma p\rightarrow J/\psi)=\sigma_{0}(W/W_{0})^{\delta}$ (with
$\sigma_{0}=81$ nb, $W_{0}=90$ GeV and $\delta=0.67$). This is labeled by the
dotted curve in figure. Moreover, a similar fit including threshold
corrections to the production cross section provided in the STARLIGHT Monte
Carlo [37] is also presented (dot-dashed curve). Now, the fit takes the form
$\sigma(\gamma p\rightarrow J/\psi)=\sigma_{0}(W/W_{0})^{\delta}(1-x_{V})^{2}$
(with $\sigma_{0}=4.06$ nb, $W_{0}=1$ GeV and $\delta=0.65$). Here,
$x_{V}=(m_{p}+m_{J/\psi})^{2}/W^{2}$. The predictions using the Reggeometric
Pomeron are consistent with experimental results. Of course, there is room for
future adjustment of the parameters using the high precision HERA data and
threshold correction, which are important at very large $y$ in the
$\omega_{-}$ contribution. In Fig. 5, the rapidity distribution for $\psi(2S)$
production is presented. We have multiplied the $J/\psi$ production cross
section by the measured value for the ratio
$R_{\psi}=\sigma(\psi(2S))/\sigma(J/\psi)$ at very low $Q^{2}$. Namely, we
have used $R_{\psi}=0.17$ [38] and the upper limit for the model. The
corresponding data description is reasonable given the simplification done.

Figure 5: Rapidity distribution of $\psi(2S)$ production in $pp$ collisions at
the LHC ($\sqrt{s_{pp}}=13$ TeV). The curve stands for the Reggeometric
Pomeron model. Data from LHCb collaboration [34].

Let us move now to the exclusive $J/\psi$ and $\psi(2S)$ production in PbPb
collisions. In the calculations we employ traditional vector dominance model
(VDM) and the classical mechanics Glauber formula for multiple scattering of
vector meson in the nuclear medium [39]. In this approach, the differential
cross section for the process, $\gamma+A\rightarrow J/\psi+A$ (nuclear
coherent scattering), is given by:

$\displaystyle\left.\frac{\mathrm{d}\sigma\left(\gamma+A\to
J/\psi+A\right)}{\mathrm{d}t}\right|_{t=0}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{em}\sigma_{tot}^{2}(J/\psi A)}{4f^{2}_{J/\psi}},$
(10) $\displaystyle\sigma_{tot}\left(J/\psi A\right)$ $\displaystyle=$
$\displaystyle\int\mathrm{d}^{2}\textbf{b}\left[1-\exp\left(-\sigma_{tot}\left(J/\psi
p\right)T_{A}\left(\textbf{b}\right)\right)\right],$ (11)

where $T_{A}(b)$ is the nuclear thickness function and $f_{J/\Psi}$ is the
vector-meson coupling (where $f^{2}_{J/\Psi}/4\pi=10.4$). The approach above
has been extended in the generalized vector dominance model (GVDM) (see for
instance the comprehensive review of Ref. [40]). The input for the Glauber
model calculation in Eq. (11) is the cross section for the process
$J/\psi+p\rightarrow J/\psi+p$ is given by:

$\displaystyle\sigma_{tot}\left(J/\psi
p\right)=\left[\frac{4f^{2}_{J/\psi}}{\alpha_{em}}\left.\frac{\mathrm{d}\sigma\left(\gamma+p\to
J/\psi+p\right)}{\mathrm{d}t}\right|_{t=0}\right]^{1/2}.$ (12)

Figure 6: Rapidity distribution of $J/\psi$ and $\psi(2S)$ production in PbPb
collisions at the LHC ($\sqrt{s_{\mathrm{NN}}}=2.76$ TeV). Curves represent
prediction of the Reggeometric Pomeron model. Data from ALICE and LHCb
collaborations [42, 43, 44, 45].

The total cross section for $\gamma+A\rightarrow J/\psi+A$ process can be
written as:

$\displaystyle\sigma_{\gamma+A\rightarrow
J/\psi+A}=\left.\frac{d\sigma(\gamma+A\rightarrow
J/\psi+A)}{dt}\right|_{t=0}\int\limits_{t_{min}}^{\infty}\mathrm{d}|t|\,\left|F_{A}\left(t\right)\right|^{2},$
(13)

where $F_{A}$ is the nuclear form factor which is the Fourier transform of the
nuclear density profile, $F_{A}(k^{2})=\int
d^{3}\vec{r}\,e^{i\vec{k}\cdot\vec{r}}\rho_{A}(\vec{r})$. In our calculation
we consider an analytic form factor given by a hard sphere of radius,
$R_{A}=1.2A^{1/3}$ fm, convoluted with a Yukawa potential with range $a$
($a=0.7$ fm) [41]:

$\displaystyle
F_{A}(|k|)=\frac{4\pi\rho_{0}}{A|k^{3}|}\left(\frac{1}{1+a^{2}k^{2}}\right)\left[\sin{(|k|R_{A})}-|k|R_{A}\cos{(|k|R_{A})}\right],$
(14)

where $A$ is the mass number, $k$ is the momentum transfer ($k^{2}=|t|)$ and
$\rho_{0}=A/\frac{4}{3}\pi R_{A}^{3}$.

The coherent production of $J/\psi$ in $AA$ collisions is easily computed in
the equivalent photon approximation (EPA). The rapidity distribution is given
by:

$\displaystyle\frac{d\sigma(A+A\rightarrow
A+V+A)}{dy}=\omega_{+}\frac{dN_{\gamma/A}(\omega_{+})}{d\omega}\sigma_{\gamma
A\rightarrow
VA}(\omega_{+})+\omega_{-}\frac{dN_{\gamma/A}(\omega_{-})}{d\omega}\sigma_{\gamma
A\rightarrow VA}(\omega_{-}),$ (15)

where the rapidity $y$ of the vector meson $V$ is related to the center-of-
mass energy of the photon-nucleus system, $W_{\gamma
A}^{2}=\sqrt{s_{\mathrm{NN}}}M_{V}e^{Y}$. The quantity
$\sqrt{s_{\mathrm{NN}}}$ is the center-of-mass energy per nucleon pair in the
AA system and $\omega_{\pm}=\frac{M_{V}}{2}e^{\pm Y}$. For simplicity, we use
the analytical photon flux corresponding to a point-like charge [10, 9]:

$\displaystyle\frac{\mathrm{d}N_{\gamma/A}(\omega)}{\mathrm{d}\omega}=\frac{2Z^{2}\alpha_{em}}{\pi\omega}\left(XK_{0}(X)K_{1}(X)-\frac{X^{2}}{2}\left[K_{1}^{2}(X)-K_{0}^{2}(X)\right]\right),$
(16)

where $Z$ is the nucleus atomic number, $K_{0,1}$ are the modified Bessel
function and $X=2R_{A}\omega/\gamma_{L}$.

Figure 7: Rapidity distribution of $J/\psi$ and $\psi(2S)$ production in PbPb
collisions at the LHC ($\sqrt{s_{\mathrm{NN}}}=5.02$ TeV). Curves represent
predictions of the Reggeometric Pomeron model. Data from ALICE and LHCb
collaborations [1, 2, 3].

In Fig. 6 results are shown for exclusive $J/\psi$ and $\psi(2S)$ production
for PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV. The experimental
measurements from ALICE [42, 43, 44] and CMS [45] collaborations are
presented. The prediction for the Reggeometric Pomeron are labeled by the
solid ($J/\psi$) and dashed ($\psi(2S)$) curves, respectively. The input for
the calculation is Eq. (6) using parameters of Ref. [30] (see Table 1). For
$J/\psi$ the result overestimates the experimental measurements. This is
expected as the value of $\sigma(J/\psi p)$ is small which leads to the rather
small nuclear suppression. This feature is discussed in details in Refs. [46,
47]. In the case of $\psi(2S)$ production the agreement with the data is
reasonable.

In Fig. 7 results are shown for the energy of $\sqrt{s_{\mathrm{NN}}}=5.02$
TeV. The experimental measurements from ALICE [1, 2] and LHCb [3]
collaborations are presented. As observed, there is some indication that the
ALICE and LHCb data at forward rapidity are mutually inconsistent. The
Reggeometric Pomeron model is in agreement with ALICE data at forward
rapidities and overestimates data at mid-rapidity. The same occurs for
$\psi(2S)$ at central rapidity. The predictions follow the trend observed at
2.76 TeV.

Motivated by the recent ALICE measurement for coherent $\rho^{0}$ production
in XeXe ultraperipheral collisions (UPC) at 5.44 TeV [52], we estimate the
rapidity distribution for the $J/\psi$ production at 5.44 and 5.86 TeV. This
is presented in Fig. 9. In this case the shadowing effects are expected to
depend on the $A$ and UPC measurements for different atomic number extend the
study of the shadowing. Along these lines was reported performance projections
for measurements that can be carried out with an Oxigen-Oxigen run during the
LHC Run 3 [53, 54]. In the impulse approximation, $d\sigma_{\gamma
A}(t=0)/dt\approx A^{2}d\sigma_{\gamma p}(t=0)/dt$, and in case of similar
energies it is expected that the XeXe cross section be a factor 0.4 smaller
than for PbPb collisions.

Finally, we present our predictions for the transverse momentum distribution
for $J/\psi$ production which has been recently measured by ALICE
collaboration [4] in PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV. The
$p_{T}$-distribution is obtained theoretically in the following way:

$\displaystyle\frac{d^{2}\sigma(A+A\rightarrow
A+V+A)}{dydp_{T}^{2}}=\omega_{+}\frac{dN_{\gamma/A}(\omega_{+})}{d\omega}\frac{d\sigma_{\gamma
A\rightarrow
VA}(\omega_{+})}{d|t|}+\omega_{-}\frac{dN_{\gamma/A}(\omega_{-})}{d\omega}\frac{d\sigma_{\gamma
A\rightarrow VA}(\omega_{-})}{d|t|},$ (17)

where the $t$-dependence is obtained directly from Eq. (13) by removing the
integration in $|t|$. In Fig. 8 the result is presented and compared to ALICE
data at $y=0$. As it can be seen, the Reggeometric Pomeron model overestimates
the data in consistency with the rapidity distribution at mid-rapidity in Fig.
7. Similar overestimation is also observed for the prediction from STARLIGHT
Monte Carlo, where the shape of the $|t|$-distribution is wider than that of
experimental measurement. The normalization and shape of data seem to be
better described with the leading-twist approximation (LTA) of nuclear
shadowing considering a low shadowing prediction (see Ref. [47]).

Figure 8: Transverse momentum of $J/\psi$ production in PbPb collisions at the
LHC ($\sqrt{s_{\mathrm{NN}}}=5.02$ TeV). Prediction of the Reggeometric
Pomeron model is compared to data from ALICE collaboration [4]. Figure 9:
Rapidity distribution of $J/\psi$ production in XeXe collisions at the LHC
($\sqrt{s_{\mathrm{NN}}}=5.44$ and $5.86$ TeV).

##  4 Conclusions

In this paper predictions for exclusive $J/\Psi$ and $\psi(2S)$ meson
production in ultra-peripheral collisions at the LHC are presented and
compared with the recent experimental data collected by different
collaborations. We update our previous studies in Refs. [26, 27, 28] by
analysing the measurements in proton-proton collisions at 13 TeV. The cross
section $\sigma(\gamma+p\rightarrow J/\psi+p)$ is described in terms of the
single-component Reggeometric Pomeron model. It describes nicely the recent
indirect measurements (see Fig. 3) and the HERA experimental data. A better
description can be obtained by further tuning the power. The LHCb data at 13
TeV are described reasonably given the band of uncertainty for the model. The
description of $\psi(2S)$ data is good. By using the Glauber model and VDM
approach, the predictions for PbPb collisions at 2.76 and 5.02 TeV are
presented. Concerning the rapidity distributions, the model tends to describe
better the forward rapidity region and overestimates the data at central
rapidity. This feature is similar to the approaches that not consider leading-
twist nuclear shadowing. Prediction for XeXe UPC collision has been provided.
We provide prediction for the transverse momentum distribution and compared it
to the first measurement of $|t|$-distribution performed by ALICE
collaboration [4].

The present study can be extended by:

  * •

the inclusion of the $t$ dependence of the differential cross section with
deviations from the usual linear parameterization (non-linear Regge
trajectories) [48];

  * •

the inclusion of other models beyond the single-component Reggeometric Pomeron
like the two-pomeron models and the soft dipole pomeron [31];

  * •

studies on incoherent vector meson production, $A+A\rightarrow A+V+A^{*}$, in
ultraperipheral heavy ion collisions [37, 20, 49, 50, 22, 51];

  * •

studies on production at low transverse momentum and its dependence on
centrality in PbPb collisions [5].

## Acknowledgements

L.J. was supported by the grant 1230 Fundamental properties of matter in
relativistic nuclear collisions and early Universe of the National Academy of
Sciences of Ukraine, 0120U100935. MVTM was supported by funding agencies CAPES
and CNPq (Brazil).

## References

  * [1] S. Acharya et al. [ALICE], Phys. Lett. B 798, 134926 (2019) doi:10.1016/j.physletb.2019.134926 [arXiv:1904.06272 [nucl-ex]].
  * [2] S. Acharya et al. [ALICE], Eur. Phys. J. C 81, no.8, 712 (2021) doi:10.1140/epjc/s10052-021-09437-6 [arXiv:2101.04577 [nucl-ex]].
  * [3] R. Aaij et al. [LHCb], [arXiv:2107.03223 [hep-ex]].
  * [4] S. Acharya et al. [ALICE], Phys. Lett. B 817, 136280 (2021) doi:10.1016/j.physletb.2021.136280 [arXiv:2101.04623 [nucl-ex]].
  * [5] R. Aaij et al. [LHCb], [arXiv:2108.02681 [hep-ex]].
  * [6] G. Baur et al., Phys. Rept. 364 (2002) 359; K. Hencken et al., Phys. Rept. 458 (2008) 1.
  * [7] J. G. Contreras and J. D. Tapia Takaki, Int. J. Mod. Phys. A 30, 1542012 (2015) doi:10.1142/S0217751X15420129
  * [8] W. Schäfer, Eur. Phys. J. A 56, 231 (2020).
  * [9] S. Klein and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 70, 323-354 (2020) doi:10.1146/annurev-nucl-030320-033923 [arXiv:2005.01872 [nucl-ex]].
  * [10] S. Klein and J. Nystrand, Phys. Rev. C 60, 014903 (1999) doi:10.1103/PhysRevC.60.014903 [arXiv:hep-ph/9902259 [hep-ph]].
  * [11] V. Guzey, E. Kryshen and M. Zhalov, Phys. Rev. C 93, no.5, 055206 (2016) doi:10.1103/PhysRevC.93.055206 [arXiv:1602.01456 [nucl-th]].
  * [12] M. Kłusek-Gawenda and A. Szczurek, Phys. Rev. C 93, no.4, 044912 (2016) doi:10.1103/PhysRevC.93.044912 [arXiv:1509.03173 [nucl-th]].
  * [13] I. Bautista, A. Fernandez Tellez and M. Hentschinski, Phys. Rev. D 94, no.5, 054002 (2016) doi:10.1103/PhysRevD.94.054002 [arXiv:1607.05203 [hep-ph]].
  * [14] M. Hentschinski and E. Padrón Molina, Phys. Rev. D 103, no.7, 074008 (2021) doi:10.1103/PhysRevD.103.074008 [arXiv:2011.02640 [hep-ph]].
  * [15] S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, JHEP 11, 085 (2013) doi:10.1007/JHEP11(2013)085 [arXiv:1307.7099 [hep-ph]].
  * [16] S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, J. Phys. G 41, 055009 (2014) doi:10.1088/0954-3899/41/5/055009 [arXiv:1312.6795 [hep-ph]].
  * [17] V. P. Gonçalves, M. V. T. Machado, B. D. Moreira, F. S. Navarra and G. S. dos Santos, Phys. Rev. D 96, no.9, 094027 (2017) doi:10.1103/PhysRevD.96.094027 [arXiv:1710.10070 [hep-ph]].
  * [18] G. Sampaio dos Santos and M. V. T. Machado, J. Phys. G 42, no.10, 105001 (2015) doi:10.1088/0954-3899/42/10/105001 [arXiv:1411.7918 [hep-ph]].
  * [19] T. Lappi and H. Mantysaari, Phys. Rev. C 83, 065202 (2011) doi:10.1103/PhysRevC.83.065202 [arXiv:1011.1988 [hep-ph]].
  * [20] T. Lappi and H. Mantysaari, Phys. Rev. C 87, no.3, 032201 (2013) doi:10.1103/PhysRevC.87.032201 [arXiv:1301.4095 [hep-ph]].
  * [21] J. Cepila, J. G. Contreras and J. D. Tapia Takaki, Phys. Lett. B 766, 186-191 (2017) doi:10.1016/j.physletb.2016.12.063 [arXiv:1608.07559 [hep-ph]].
  * [22] J. Cepila, J. G. Contreras and M. Krelina, Phys. Rev. C 97, no.2, 024901 (2018) doi:10.1103/PhysRevC.97.024901 [arXiv:1711.01855 [hep-ph]].
  * [23] A. Łuszczak and W. Schäfer, Phys. Rev. C 99, no.4, 044905 (2019) doi:10.1103/PhysRevC.99.044905 [arXiv:1901.07989 [hep-ph]].
  * [24] M. B. Gay Ducati and S. Martins, Phys. Rev. D 97, no.11, 116013 (2018) doi:10.1103/PhysRevD.97.116013 [arXiv:1804.09836 [hep-ph]].
  * [25] L. Jenkovszky, R. Schicker and I. Szanyi, Int. J. Mod. Phys. E 27, no.08, 1830005 (2018) doi:10.1142/S0218301318300059 [arXiv:1902.05614 [hep-ph]].
  * [26] R. Fiore, L. Jenkovszky, V. Libov and M. Machado, Teor. Mat. Fiz. 182, no.1, 171-181 (2014) doi:10.1007/s11232-015-0252-8 [arXiv:1408.0530 [hep-ph]].
  * [27] R. Fiore, L. Jenkovszky, V. Libov, M. V. T. Machado and A. Salii, [arXiv:1506.01990 [hep-ph]]. Contribution to: Diffraction 2014.
  * [28] R. Fiore, L. Jenkovszky, V. Libov, M. V. T. Machado and A. Salii, AIP Conf. Proc. 1654, no.1, 090002 (2015) doi:10.1063/1.4916009
  * [29] M. Capua et al., Phys. Lett. B645 (1997) 161, hep-ph/0605319.
  * [30] S. Fazio, R. Fiore, A. Lavorini, L. Jenkovszky and A. Salii, Acta Phys. Polon. B 44, 1333-1353 (2013), [arXiv:1304.1891 [hep-ph]].
  * [31] S. Fazio, R. Fiore, L. Jenkovszky, and A. Salii, Phys. Rev. D 90, no.1, 016007 (2014), [arXiv:1312.5683 [hep-ph]].
  * [32] B. B. Abelev et al. [ALICE], Phys. Rev. Lett. 113, no.23, 232504 (2014) doi:10.1103/PhysRevLett.113.232504 [arXiv:1406.7819 [nucl-ex]].
  * [33] S. Acharya et al. [ALICE], Eur. Phys. J. C 79, no.5, 402 (2019) doi:10.1140/epjc/s10052-019-6816-2 [arXiv:1809.03235 [nucl-ex]].
  * [34] R. Aaij et al. [LHCb], JHEP 10, 167 (2018) doi:10.1007/JHEP10(2018)167 [arXiv:1806.04079 [hep-ex]].
  * [35] S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, J. Phys. G 44, no.3, 03LT01 (2017) doi:10.1088/1361-6471/aa56ea [arXiv:1611.03711 [hep-ph]].
  * [36] C. Alexa et al. [H1], Eur. Phys. J. C 73, no.6, 2466 (2013) doi:10.1140/epjc/s10052-013-2466-y [arXiv:1304.5162 [hep-ex]].
  * [37] S. R. Klein, J. Nystrand, J. Seger, Y. Gorbunov and J. Butterworth, Comput. Phys. Commun. 212, 258-268 (2017) doi:10.1016/j.cpc.2016.10.016 [arXiv:1607.03838 [hep-ph]].
  * [38] H. Abramowicz et al. [ZEUS], Nucl. Phys. B 909, 934-953 (2016) doi:10.1016/j.nuclphysb.2016.06.010 [arXiv:1606.08652 [hep-ex]].
  * [39] T. H. Bauer, R. D. Spital, D. R. Yennie and F. M. Pipkin, Rev. Mod. Phys. 50, 261 (1978) [erratum: Rev. Mod. Phys. 51, 407 (1979)] doi:10.1103/RevModPhys.50.261
  * [40] L. Frankfurt, M. Strikman and M. Zhalov, Acta Phys. Polon. B 34, 3215-3254 (2003) [arXiv:hep-ph/0304301 [hep-ph]].
  * [41] K. T. R. Davies and J. R. Nix, Phys. Rev. C 14, 1977-1994 (1976) doi:10.1103/PhysRevC.14.1977
  * [42] E. Abbas et al. [ALICE], Eur. Phys. J. C 73, no.11, 2617 (2013) doi:10.1140/epjc/s10052-013-2617-1 [arXiv:1305.1467 [nucl-ex]].
  * [43] B. Abelev et al. [ALICE], Phys. Lett. B 718, 1273-1283 (2013) doi:10.1016/j.physletb.2012.11.059 [arXiv:1209.3715 [nucl-ex]].
  * [44] J. Adam et al. [ALICE], Phys. Lett. B 751, 358-370 (2015) doi:10.1016/j.physletb.2015.10.040 [arXiv:1508.05076 [nucl-ex]].
  * [45] V. Khachatryan et al. [CMS], Phys. Lett. B 772, 489-511 (2017) doi:10.1016/j.physletb.2017.07.001 [arXiv:1605.06966 [nucl-ex]].
  * [46] V. Guzey, E. Kryshen, M. Strikman and M. Zhalov, Phys. Lett. B 726, 290-295 (2013) doi:10.1016/j.physletb.2013.08.043 [arXiv:1305.1724 [hep-ph]].
  * [47] V. Guzey, E. Kryshen, M. Strikman and M. Zhalov, Phys. Lett. B 816, 136202 (2021) doi:10.1016/j.physletb.2021.136202 [arXiv:2008.10891 [hep-ph]].
  * [48] S. Fazio, R. Fiore, L. Jenkovszky and A. Lavorini, Phys. Rev. D 85, 054009 (2012) doi:10.1103/PhysRevD.85.054009 [arXiv:1109.6374 [hep-ph]].
  * [49] T. Toll and T. Ullrich, Comput. Phys. Commun. 185, 1835-1853 (2014) doi:10.1016/j.cpc.2014.03.010 [arXiv:1307.8059 [hep-ph]].
  * [50] M. Strikman, M. Tverskoy and M. Zhalov, Phys. Lett. B 626, 72-79 (2005) doi:10.1016/j.physletb.2005.08.083 [arXiv:hep-ph/0505023 [hep-ph]].
  * [51] M. C. Traini and J. P. Blaizot, Eur. Phys. J. C 79, no.4, 327 (2019) doi:10.1140/epjc/s10052-019-6826-0 [arXiv:1804.06110 [hep-ph]].
  * [52] S. Acharya et al. [ALICE], Phys. Lett. B 820, 136481 (2021) doi:10.1016/j.physletb.2021.136481 [arXiv:2101.02581 [nucl-ex]].
  * [53] ALICE Collaboration, ALICE physics projections for a short oxygen-beam run at the LHC, ALICE-PUBLIC-2021-004.
  * [54] J. Brewer, A. Mazeliauskas and W. van der Schee, [arXiv:2103.01939 [hep-ph]].


# The estimation of bias and variance in clustering coefficient streaming
algorithms

Roohollah Etemadi, Jianguo Lu  School of Computer Science, University of
Windsor
Windsor, ON, Canada
etemadir<EMAIL_ADDRESS>
###### Abstract

Clustering coefficient is one of the most important metrics to understand the
complex structure of networks. This paper addresses the estimation of
clustering coefficient in network streams. There have been substantial work in
this area, most of conducting empirical comparisons of various algorithms. The
variance and the bias of the estimators have not been quantified. Starting
with a simple yet powerful streaming algorithm, we derived the variance and
bias for the estimator, and the estimators for the variances and bias. More
importantly, we simplify the estimators so that it can be used in practice.
The variance and bias estimators are verified extensively on 49 real networks.

###### Index Terms:

Estimation; Clustering Coefficient; Streaming algorithms; Bias; Variance.

##  I Introduction

The amount of data in digital world is growing faster than before in the age
of the Internet. This deluge of data results in massive networks with the size
of billions. Recently analyzing such real-world networks has captivated a
great attention among practitioners and scholars. Clustering coefficient
(hereafter $\mathcal{C}$) is one of the most important metrics to analyze such
networks. It has been used in many applications including graph clustering,
community detection, spam detection, link prediction, wireless and ad-hoc
networks analysis, microarray data and DNA sequence analysis, word-learning,
risk analysis, etc. Computing $\mathcal{C}$ on large networks is an intensive
task, i.e. the time complexity of the state-of-the-art method is
$\Theta(N^{1.41})$, here $N$ is the number of nodes in the network [1]. Thus,
sampling-based algorithms are indispensable.

This paper focuses on the streaming algorithms where data arrive sequentially
in an arbitrary order. A number of techniques have been proposed to estimate
$\mathcal{C}$ in such a streaming model [2, 3, 4, 5]. The aim is to provide an
accurate estimation using limited memory budgets over single or multi passes
of the stream of edges of networks.

Despite extensive work in this area, there is a lack of formal analyses of the
algorithms, in particular the lack of estimators for the variance and bias of
the estimation algorithms. Most algorithms are compared empirically, thus
their performances are often data dependent. Some algorithms do discuss the
bias and variance problems, but they all fail to give the exact formulas.
E.g., [2] [3] gives a variance that is generic for every problem; [2, 3, 4, 5]
mentioned about the existence of the bias.

In addition to the bias and variance formulas, a more important issue is their
estimation without the global data. In order to have the confidence interval
for an estimate, we need to know the variance and bias during the sampling and
estimation process, not the variance and bias that is derived from the
knowledge of the global data. I.e., we need to estimate the variance and bias.

This paper derived the variance and bias so that algorithms can be compared
analytically. The formulas are along and tedious, hard to be applied in
practice. A more interesting result of the paper is the simplification of the
formulas so that it can be applied easily. Based on the assumption that the
data is big, we show that the RSE (Relative standard error, the normalized
square root of variance) can be estimated by the reciprocal of the square-root
of the triangles observed in samples. In other words, the accuracy grows at
the rate of $\sqrt{\Delta_{g}}$, where $\Delta_{g}$ is the observed sample
closed wedges in the streaming process.

The simplified result is especially important in the era of big data. When
networks are small, we can hardly predict the behaviour of the sampling
algorithms. However, when the data is big, the performance can be
characterized simply by the closed wedges observed in the process. This result
is also supported by our extensive experiments on 49 real network of different
size and structure.

The second contribution is the correction of the bias of the traditional
estimator. This is an extension of the work from direct sampling [6] to
streaming algorithms. The existence of the bias in streaming algorithms was
also observed in [3], but it is not quantified, thus not corrected. We proved
that the bias depends on the structure of the sample.

##  II Background and related work

###  II-A Notations

Suppose $G(\mathcal{V},\mathcal{E})$ be a simple undirected graph, where
$\mathcal{V}$ is the set of nodes, and $\mathcal{E}$ the set of edges. For
simplicity, we assume that the graph is not a multi-graph, and does not have
self-loops. Let $N=|\mathcal{V}|$, $M=|\mathcal{E}|$. Let $1,2,\dots,M$ be the
labels of the edges in $\mathcal{E}$ according to their occurrence order in
the stream. E.g. edge $M$ is the last edge in the stream. A wedge
$\mathcal{W}$ is a path $(u,v,w)$ of length two, where nodes
$u,v,w\in\mathcal{V}$, and edges $(u,v)\in\mathcal{E}$ and
$(v,w)\in\mathcal{E}$. A wedge $\mathcal{W}$ is closed if
$(u,w)\in\mathcal{E}$. Otherwise it is open. A closed wedge $\mathcal{W}$ is
also called a triangle. Note that each triangle has three (closed) wedges. Let
$\Delta$ denote the number of closed wedges and $\Lambda$ the count of wedges
in $G$.

###  II-B Related work

A number of sampling-based methods have been proposed in recent years to
estimated clustering coefficient in networks. A straightforward method is
wedge sampling [7]. It selects wedges uniformly at random and the fraction of
closed ones is used as an estimation for $\mathcal{C}$. To sample a random
wedge, two passes over an edge stream are required in a streaming model.
Furthermore, one additional pass over edge stream is needed to check the
closeness of the random wedge [8] [9]. Therefore, sampling a random wedge from
a large graph is an intensive task. A scalable technique to estimate
$\mathcal{C}$ is edge sampling [2, 3, 4]. It samples some edges uniformly at
random and create a subgraph $g$. Then, the number of wedges and closed ones
in $g$ are used to estimate $\mathcal{C}$.

Estimating $\mathcal{C}$ using edge-based sampling methods is biased. It is
also noticed in [2][3][5]. Recently, authors in [6] quantify the bias in none-
streaming model when edges of original graph are randomly accessible. However,
to the best of our knowledge, there is no study to quantify the bias in the
streaming model.

Input: $p$

Output: $\widehat{\Delta}$, RSE($\widehat{\Delta}$)

1 begin

2 $\Delta_{g}=0$, $\Lambda_{g}=0$, $g=\\{\phi\\}$.

3 while  _new edge e_ do

4 Add $e$ into $g$ with probability $p$.

5 foreach  _wedge $w$ formed using $e$ and edges in $g$_ do

6 $\Lambda_{g}+=1.$

7 end foreach

8 foreach  _wedge $w\in g$ closed by $e$_ do

9 $\Delta_{g}+=1.$

10 end foreach

11

12 end while

13 $\widehat{\mathcal{C}}=3\Delta_{g}(p\Lambda_{g})^{-1}$.

14 $\widehat{RSE}(\widehat{\mathcal{C}}$)$\approx\Delta_{g}^{-1/2}$.

15 end

16

Algorithm 1 Naive edge sampling (NES)

###  II-C The Algorithm

In order to discuss the variance and bias problem rigorously, we give a simple
edge-based streaming algorithm (hereafter we call it NES–Naive Edge Streaming)
as described in Alg. 1. It is adapted from a direct sampling algorithm [6] to
the stream model. It is also a special case of [3] where all edges are sampled
with the same priority.

In NES, edges arrive in an arbitrary sequence over an edge stream of the
original graph $G$. Over the stream, NES samples some edges uniformly at
random with an equal probability $p$, and adds them into the subgraph $g$.
Once a new edge $e$ arrives, the number of wedges closed by $e$ in $g$ is
counted. Let $\Delta_{g}$ denote the number of such closed wedges observed in
the stream. At the same time, the number of wedges formed using $e$ and edges
in $g$ is enumerated. Let $\Lambda_{g}$ be the total number of such wedges.
Then, $\mathcal{C}$ is estimated using $\Delta_{g}$ and $\Lambda_{g}$. A key
difference with [6] is that $\Delta_{g}$ and $\Lambda_{g}$ are counted along
the streaming process, not after the sampling.

Given a subgraph $g$ of $G$, for every wedge $(u,v,w)$ in $g$, we check
whether or not $(u,w)$ is in the rest of the stream. For every closed wedge in
$G$, if its two edges are sampled in $g$, the probability of observing the
third edge in the rest of the stream is $\frac{1}{3}$. Thus, the probability
of identifying a closed wedges by NES is $\frac{1}{3}p^{2}$. Suppose
$\delta_{i}$ be an indicator for the $i^{th}$ closed wedge in $G$. Variable
$\delta_{i}$ is 1 when two edges of the $i^{th}$ closed wedge are sampled and
the third edge is observed in the rest of the stream; otherwise it is 0.
Recall that $\Delta_{g}$ is the number of closed wedges identified by NES over
a single pass on an edge stream of $G$. The expectation of $\Delta_{g}$ is
$\mathbb{E}(\Delta_{g})=\mathbb{E}(\sum_{i=1}^{\Delta}\delta_{i})=\sum_{i=1}^{\Delta}\mathbb{E}(\delta_{i})=\sum_{i=1}^{\Delta}\frac{1}{3}p^{2}=\frac{1}{3}p^{2}\Delta$.
Thus, an unbiased estimator for $\Delta$ is given by
$\widehat{\Delta}=\frac{3\Delta_{g}}{p^{2}}$. Similar estimators for $\Delta$
have also been proposed in [10, 11, 12].

Next, we give an unbiased estimator for the number of wedges in $G$. To sample
a wedge by NES, one of its two edges needs to be added into $g$ with
probability $p$, and its second edge is required to be observed in the rest of
the stream. Thus, the probability of identifying a wedge based on $g$ is $p$.
Suppose $\lambda_{i}$ be an indicator for the $i^{th}$ wedge in the input
graph. Clearly, $\lambda_{i}$ is 1 when its two edges are observed; otherwise
it is 0. Recall that $\Lambda_{g}$ is the number of wedges identified based on
$g$ by NES. Its expectation is
$\mathbb{E}(\Lambda_{g})=\mathbb{E}(\sum_{i=1}^{\Lambda}\lambda_{i})=\sum_{i=1}^{\Lambda}\mathbb{E}(\lambda_{i})=\sum_{i=1}^{\Lambda}p=p\Lambda$.
Thus, an unbiased estimation for $\Lambda$ is given by
$\widehat{\Lambda}=\frac{\Lambda_{g}}{p}.$

Now we can use the unbiased estimators for $\Delta$ and $\Lambda$ to estimate
$\mathcal{C}$. Although both $\widehat{\Delta}$ and $\widehat{\Lambda}$ are
unbiased, the following estimator is biased as we will correct it later in
this paper.

$\widehat{\mathcal{C}}=\frac{\widehat{\Delta}}{\widehat{\Lambda}}=\frac{3p\Delta_{g}}{p^{2}\Lambda_{g}}=\frac{3\Delta_{g}}{p\Lambda_{g}}.$
(1)

##  III Estimator of the variance

We derive the variance of the estimator using the Delta method. Applying $var$
on Eq. 1, we get

$\displaystyle var(\widehat{\mathcal{C}})$
$\displaystyle=var\bigg{(}\frac{3\;\Delta_{g}}{p\;\Lambda_{g}}\bigg{)}=\frac{9}{p^{2}}var\bigg{(}\frac{\Delta_{g}}{\Lambda_{g}}\bigg{)}.$
(2)

Applying Taylor expansion in the neighbourhood of $(a,b)$, we have:

$\displaystyle var\biggl{(}\frac{\Delta_{g}}{\Lambda_{g}}\biggr{)}\approx$
$\displaystyle\frac{1}{b^{2}}var(\Delta_{g})+\frac{a^{2}}{b^{4}}var(\Lambda_{g})-\frac{2a}{b^{3}}cov(\Delta_{g},\Lambda_{g}).$

Let $a=\mathbb{E}(\Delta_{g})$ and $b=\mathbb{E}(\Lambda_{g})$, and using the
fact that $\mathbb{E}(\Delta_{g})=\frac{1}{3}\Delta p^{2}$, and
$\mathbb{E}(\Lambda_{g})=\Lambda p$, we obtain the variance and present it in
the form of relative standard error (RSE=$\sqrt{var}/\mathcal{C}$) as follows:

$\displaystyle RSE(\widehat{\mathcal{C}})$
$\displaystyle\approx\biggl{[}\frac{9\;var(\Delta_{g})}{\Delta^{2}p^{4}}+\frac{var(\Lambda_{g})}{\Lambda^{2}p^{2}}-\frac{6\;cov(\Delta_{g},\Lambda_{g})}{\Delta\Lambda
p^{3}}\biggr{]}^{-1}.$

The RSE depends on the variance of $\Delta_{g}$ and $\Lambda_{g}$, and the
covariance between them. Note that this is where [2, 3] stops. We continue the
derivation of the variances and covariances in Appendix. When the networks are
large, $p$ is a very small value. Hence, we can assume that $1-p\approx
1-p^{2}\approx 1$, and the results of the Lemmas 1, 2, and 3 are simplified as
follows:

$\displaystyle
var(\Delta_{g})\approx\frac{1}{3}\Delta\;p^{2}+\frac{16}{15}\Phi\;p^{3},$ (3)
$\displaystyle var(\Lambda_{g})\approx\frac{2}{3}\Psi\;p,$ (4) $\displaystyle
cov(\Delta_{g},\Lambda_{g})\approx\frac{5}{12}\Omega^{\prime}\;p^{2}.$ (5)

Here $\Phi$ is the number of pairs of dependent triangles, $\Psi$ is the
number of pairs of shared wedges, and $\Omega^{\prime}$ is the number of pairs
of wedges and triangles with one common edge in the original graph $G$.
Replace Eq.s 3, Eq. 4, and Eq. 5 in the RSE above and after some math
simplifications we obtain the RSE as:

$\displaystyle RSE(\widehat{\mathcal{C}})$
$\displaystyle\approx\biggl{[}\frac{3}{\Delta p^{2}}+\frac{48\Phi
p^{3}}{5\Delta^{2}p^{4}}+\frac{2\Psi\;p}{3\Lambda^{2}p^{2}}-\frac{5\Omega^{\prime}p^{2}}{2\Delta\Lambda
p^{3}}\biggr{]}^{-1}.$ (6)

The RSE in Eq. 6 depends on the properties of the original graph $G$, i.e.
$\Delta$, $\Lambda$, $\Phi$, $\Psi$, and $\Omega^{\prime}$. Note that those
properties are unknown for the third party. However, practitioners need to
know the error bound of the estimator using the properties in the sample. To
do so, we give the estimation of the variance as follows.

Based on NES sampling scheme we have $\mathbb{E}(\Delta_{g})=\frac{1}{3}\Delta
p^{2}$, $\mathbb{E}(\Lambda_{g})=\Lambda p$,
$\mathbb{E}(\Phi_{g})=\frac{8}{15}\Phi p^{3}$,
$\mathbb{E}(\Psi_{g})=\frac{2}{6}\Psi p$, and
$\mathbb{E}(\Omega^{\prime}_{g})=\frac{5}{12}\Omega^{\prime}p^{2}$. Thus,
substitute $\Delta$, $\Lambda$, $\Phi$, $\Psi$, and $\Omega^{\prime}$ in Eq. 6
by their estimations, and after some math work we obtain the estimator of the
RSE as:

$\displaystyle\widehat{RSE}(\widehat{\mathcal{C}})$
$\displaystyle\approx\biggl{[}\frac{1}{\Delta_{g}}+\frac{2\Phi_{g}}{\Delta^{2}_{g}}+\frac{2\Psi_{g}}{\Lambda^{2}_{g}}-\frac{2\Omega^{\prime}_{g}}{\Delta_{g}\Lambda_{g}}\biggr{]}^{-1}.$
(7)

where $\Phi_{g}$ is the number of pairs of dependent triangles, $\Psi_{g}$ is
the number of pairs of shared wedges, and $\Omega^{\prime}_{g}$ is the number
of pairs of wedges and triangles with one common edge observed based on
subgraph $g$. The RSE in Eq. 7 hangs on several variables in the sample, i.e.
$\Delta_{g}$, $\Lambda_{g}$, $\Phi_{g}$, $\Psi_{g}$, and
$\Omega^{\prime}_{g}$. To have better understating the RSE of the estimator,
we need to simplify it further more. We claim that when sampling probability
$p$ is small, the first term in Eq. 7, i.e. $\Delta_{g}^{-1}$, is dominant.
Therefore, we give the following Theorem as a simplified estimator for the RSE
of $\mathcal{C}$.

######  Theorem 1.

When $p$ is small, the RSE of $\widehat{\mathcal{C}}$ is approximated by

$\widehat{RSE}(\widehat{\mathcal{C}})\approx\Delta_{g}^{-1/2}.$ (8)

##  IV The bias-corrected estimator

To quantify the bias, we apply the expectation on $\widehat{\mathcal{C}}$:

$\displaystyle\mathbb{E}\big{(}\widehat{\mathcal{C}}\;\big{)}=\mathbb{E}\biggl{(}\frac{3\Delta_{g}}{p\Lambda_{g}}\biggr{)}=\frac{3}{p}\mathbb{E}\biggl{(}\frac{\Delta_{g}}{\Lambda_{g}}\biggr{)}.$
(9)

The approximation of $\mathbb{E}\big{(}\Delta_{g}/\Lambda_{g}\big{)}$ using
the quadratic Taylor expansion of $\Delta_{g}/\Lambda_{g}$ in the neighborhood
of $(a,b)$ is:

$\displaystyle\mathbb{E}\biggl{(}\frac{\Delta_{g}}{\Lambda_{g}}\biggr{)}\approx$
$\displaystyle\frac{a}{b}+\frac{a}{b^{3}}var(\Lambda_{g})-\frac{1}{b^{2}}cov(\Delta_{g},\Lambda_{g}).$
(10)

Replace $a=\mathbb{E}\big{(}\Delta_{g}\big{)}$ and
$b=\mathbb{E}\big{(}\Lambda_{g}\big{)}$ and take
$\mathbb{E}\big{(}\Delta_{g}\big{)}=\Delta p^{2}/3$ and
$\mathbb{E}\big{(}\Lambda_{g}\big{)}=\Lambda p$; substitute Eq. 10 in Eq. 9,
we obtain:

$\displaystyle\mathbb{E}\big{(}\widehat{\mathcal{C}}\big{)}\approx$
$\displaystyle\frac{\Delta}{\Lambda}\biggl{(}1+\frac{var(\Lambda_{g})}{\Lambda^{2}p^{2}}-\frac{3cov(\Delta_{g},\Lambda_{g})}{\Delta\Lambda
p^{3}}\biggr{)}.$ (11)

Note that $\mathcal{C}=\Delta/\Lambda$. Thus, we get:

$\displaystyle\mathbb{E}\big{(}\widehat{\mathcal{C}}\big{)}\approx$
$\displaystyle\mathcal{C}\biggl{(}1+\frac{var(\Lambda_{g})}{\Lambda^{2}p^{2}}-\frac{3cov(\Delta_{g},\Lambda_{g})}{\Delta\Lambda
p^{3}}\biggr{)}.$ (12)

To quantify the bias of $\widehat{\mathcal{C}}$, the relative bias, i.e. RB=
$(\mathbb{E}\big{(}\widehat{\mathcal{C}}\big{)}-\mathcal{C})/\mathcal{C}$, is
used. By substituting Eq. 12 in RB, and simplifying the result we obtained:

$\displaystyle
RB\approx\frac{var(\Lambda_{g})}{\Lambda^{2}p^{2}}-\frac{3cov(\Delta_{g},\Lambda_{g})}{\Delta\Lambda
p^{3}}.$ (13)

The bias depends on the variance of $\Lambda_{g}$ and the covariance between
$\Delta_{g}$ and $\Lambda_{g}$. Using the same treatment in the variance of
$\widehat{\mathcal{C}}$ for the $var(\Lambda_{g})$ and
$cov(\Delta_{g},\Lambda_{g})$, i.e. Eq. 4 and Eq. 5, and replacing them in Eq.
13 we have

$\displaystyle RB\approx$
$\displaystyle\frac{2\Psi\;p}{3\Lambda^{2}p^{2}}-\frac{5\Omega^{\prime}\;p^{2}}{4\Delta\Lambda
p^{3}}.$ (14)

Thus, we give our result in the following theorem.

######  Theorem 2.

The RB of the estimator is approximated by

$\displaystyle
RB\approx\frac{1}{p}\biggl{[}\frac{2\Psi}{3\Lambda^{2}}-\frac{5\Omega^{\prime}}{4\Delta\Lambda}\biggr{]}.$
(15)

and its estimation is

$\displaystyle\widehat{RB}\approx\frac{2\Psi_{g}}{\Lambda_{g}^{2}}-\frac{\Omega^{\prime}_{g}}{\Delta_{g}\Lambda_{g}}.$
(16)

Figure 1:  The observed RSEs of $\widehat{\mathcal{C}}$ fit perfectly the
estimated ones in representative graphs. Estimated RSEs are obtained based on
Eq. 7.

Using $\widehat{RB}$, a bias-corrected estimator of $\widehat{\mathcal{C}}$ is

$\displaystyle\widehat{\mathcal{C}}^{+}=\frac{\widehat{\mathcal{C}}}{1+\widehat{RB}}.$
(17)

TABLE I: Properties of the networks in our experiments, sorted by graph size
$N$.

Network | N$({\times 10^{6}})$ | M | $\mathcal{C}$ | $\Delta({\scriptsize\times 10^{6}})$ | $\Lambda(\times 10^{9})$ | $\Phi(\times 10^{9})$ | $\Psi(\times 10^{12})$ | $\Omega^{\prime}(\times 10^{10})$ | Description
---|---|---|---|---|---|---|---|---|---
Ego-facebook1 | 0.004 | 88,234 | 0.519 | 4.8 | 0.009 | 0.2 | 0.003 | 0.1 | OSN
CA-GrQc 1 | 0.005 | 14,484 | 0.629 | 0.1 | 0.0002 | 0.002 | 0.00001 | 0.0009 | Collaboration
Wiki-vote1 | 0.007 | 100,762 | 0.125 | 1.8 | 0.014 | 0.04 | 0.006 | 0.08 | OSN
AstroPh 2 | 0.01 | 198,050 | 0.31 | 4.0 | 0.012 | 0.07 | 0.002 | 0.07 | Citation
CA-CondMat1 | 0.02 | 93,439 | 0.264 | 0.5 | 2 | 0.002 | 0.0001 | 0.002 | Coauthorship
HepPh 2 | 0.02 | 3,148,447 | 0.279 | 587 | 2 | 92 | 4 | 106 | Coauthorship
Enron-email2 | 0.03 | 183,831 | 0.085 | 2 | 0.025 | 0.03 | 0.01 | 0.09 | E-communication
Brightkite1 | 0.05 | 214,078 | 0.110 | 1.4 | 0.013 | 0.02 | 0.005 | 0.03 | OSN
Facebook 2 | 0.06 | 817,035 | 0.147 | 10.5 | 0.07 | 0.1 | 0.02 | 0.2 | OSN
Epinions 2 | 0.07 | 405,740 | 0.065 | 4.8 | 0.07 | 0.1 | 0.06 | 0.3 | OSN
Slashdot-Zoo 2 | 0.07 | 467,731 | 0.023 | 1.6 | 0.06 | 0.02 | 0.05 | 0.09 | OSN
Prosper 2 | 0.08 | 3,330,022 | 0.003 | 3.4 | 1.1 | 0.06 | 1 | 0.5 | Interaction
Livemocha 2 | 0.1 | 2,193,083 | 0.014 | 10.0 | 0.716 | 0.1 | 0.8 | 1 | OSN
Douban 2 | 0.1 | 327,162 | 0.01 | 0.1 | 0.011 | 0.0001 | 0.001 | 0.001 | OSN
Gowalla 1 | 0.1 | 950,327 | 0.023 | 6.8 | 0.290 | 0.1 | 2 | 0.7 | OSN
Libimseti 2 | 0.2 | 17,233,142 | 0.007 | 207 | 28 | 19 | 262 | 151 | OSN
Digg 2 | 0.2 | 1,548,126 | 0.061 | 42 | 0.69 | 3 | 2 | 6 | OSN
Web-Stanford 2 | 0.2 | 1,992,636 | 0.008 | 33 | 3.94 | 9 | 75 | 17 | Web graph
Dblp-Coau1 | 0.3 | 1,049,866 | 0.306 | 6 | 0.021 | 0.1 | 0.001 | 0.06 | Coauthorship
Web-NotreDame1 | 0.3 | 1,090,108 | 0.087 | 26 | 0.304 | 1.5 | 1 | 1 | Web graph
Amazon1 | 0.3 | 925,872 | 0.205 | 2 | 0.009 | 0.003 | 0.0005 | 0.004 | Co-purchasing
Actor 2 | 0.3 | 15,038,083 | 0.166 | 1,040 | 6.26 | 91 | 10 | 157 | Collaboration
Citeseer2 | 0.3 | 1,736,145 | 0.049 | 4 | 0.081 | 0.01 | 0.02 | 0.06 | Citation
Dogster 2 | 0.4 | 8,543,549 | 0.014 | 250 | 17 | 42 | 378 | 191 | OSN
Catster 2 | 0.6 | 15,695,166 | 0.028 | 1,969 | 69 | 1,017 | 3,651 | 1,528 | OSN
Web-Berkeley 2 | 0.6 | 6,649,470 | 0.0069 | 194 | 27.9 | 105 | 1,148 | 176 | Web graph
Web-Google 2 | 0.8 | 4,322,051 | 0.055 | 40 | 0.727 | 0.6 | 1 | 1 | Web graph
Youtube1 | 1.1 | 2,987,624 | 0.006 | 9 | 1 | 0.2 | 17 | 2 | OSN
Dblp2 | 1.3 | 5,362,414 | 0.170 | 36 | 0.214 | 0.4 | 0.05 | 0.4 | Coauthorship
Hyves 2 | 1.4 | 2,777,419 | 0.001 | 2 | 1.4 | 0.02 | 32 | 0.1 | OSN
Wiki-Polish 2 | 1.5 | 42,188,631 | 0.01 | 3,402 | 308 | 697 | 62,290 | 3,466 | Web graph
Trec-wt10g 2 | 1.6 | 6,679,248 | 0.014 | 63 | 4.3 | 11 | 50 | 13 | Web graph
Wiki-Japanese 2 | 1.6 | 56,231,610 | 0.021 | 3,863 | 180 | 685 | 12,652 | 2,784 | Web graph
Pokec 2 | 1.6 | 22,301,964 | 0.046 | 97 | 2.08 | 0.7 | 3 | 2 | OSN
As-skitter1 | 1.6 | 11,095,298 | 0.005 | 86 | 16 | 20 | 291 | 48 | Internet topology
Hudong 2 | 1.9 | 14,428,382 | 0.003 | 64 | 18.7 | 5 | 357 | 20 | Web graph
Hollywood 3 | 1.9 | 24,337,642 | 0.152 | 614 | 4 | 27 | 2 | 37 | OSN
Baidu 2 | 2.1 | 17,014,946 | 0.002 | 75 | 30.8 | 4 | 1,717 | 69 | Web graph
Flicker2 | 2.3 | 22,838,276 | 0.107 | 2,512 | 23 | 613 | 120 | 1,079 | OSN
Flixster 2 | 2.5 | 7,918,801 | 0.013 | 23 | 1.7 | 0.3 | 1 | 1 | OSN
Wiki-Russian2 | 2.8 | 63,058,425 | 0.015 | 5,697 | 370 | 1,180 | 56,397 | 5,420 | Web graph
Wiki-Franch 2 | 3.0 | 83,455,052 | 0.015 | 6,843 | 455 | 4,237 | 54,242 | 7,677 | Web graph
Orkut2 | 3.0 | 117,185,083 | 0.041 | 1,882 | 45 | 67 | 320 | 347 | OSN
USpatent 2 | 3.7 | 16,518,947 | 0.067 | 22 | 0.33 | 0.08 | 0.02 | 0.1 | Citation
LiveJournal1 | 3.9 | 34,681,189 | 0.125 | 533 | 4 | 39 | 7 | 30 | OSN
Web-Arabic3 | 22 | 553,903,073 | 0.031 | 110,686 | 3,531 | 112,260 | 986,071 | 91,940 | Web graph
Twitter2 | 41 | 1,202,513,046 | 0.0008 | 104,474 | 123,435 | 176,266 | 154,818,391 | 956,719 | OSN
MicrosoftAc.G.4 | 46 | 528,463,861 | 0.015 | 1,734 | 115 | 19 | 4,776 | 107 | Citation
Friendster2 | 65 | 1,806,067,135 | 0.017 | 12,521 | 720 | 185 | 1,138 | 1,284 | OSN

  * 1

SNAP: [13]

  * 2

Konect: [14]

  * 3

[15, 16]

  * 4

[17]

The bias-corrected estimator depends on the structure of sampled graph $g$.
Therefore, using a single pass over an edge stream of the input graph we need
to count $\Delta_{g}$, $\Lambda_{g}$, $\Psi_{g}$, and $\Omega^{\prime}_{g}$.
Recall that $\Psi_{g}$ is the number of pairs of shared wedges and
$\Omega^{\prime}_{g}$ is the number of pairs of wedges and triangles with one
common edge observed based on subgraph $g$. The authors in [6] have shown that
the cost of counting $\Psi_{g}$ and $\Omega^{\prime}_{g}$ is the same as
counting $\Lambda_{g}$ and $\Delta_{g}$. Note that the meaning of
$\Omega^{\prime}$ is different from the one in [6], i.e. $\Omega$ (the number
of pairs of wedges and closed ones with a common edge) and here it is
$\Omega^{\prime}=(\Omega-6\Delta)/2$.

Figure 2:  The estimated RSEs obtained based on Eq. 8 are apt estimations for
the observed RSEs of $\widehat{\mathcal{C}}$. Figure 3:  The observed RBs of
$\widehat{\mathcal{C}}$ (biased estimator) support our estimations of RB based
on Eq. 16.  Figure 4:  Our biased-corrected $\widehat{\mathcal{C}}^{+}$
removes the bias perfectly.

##  V Experiments

The code along with all the data are publicly available at
myweb.cs.uwindsor.ca/~etemadir/cbiasstream. The observed RSEs and RBs are
obtained over $k$ independent runs of the algorithm, where $k=1000$ except for
the three largest networks with $k=100$ for RSE, and k=50,000 for RB. Carrying
out the experiments, including obtaining the ground truth of the network
properties, is computationally expensive. The experiments were accomplished on
two big servers each with 24 cores and 256 GB RAM.

Figure 5:  Comparison of sample sizes, i.e. number of sampled edges, of the
methods when RSE=0.2. The sample size of NES is comparable with GPS-In [3].

###  V-A Datasets

Our analytical results show that the variance and bias hangs on the structure
of networks. Due to the fact that different networks have different
structures, the bias and variance vary greatly from graph to graph.
Furthermore, several approximations were made in the derivation of the
theorems. To find out the patterns behind, we need to experiment extensively
with graphs from variety of domains with different sizes and structures. We
use 49 real network graphs in total. The networks are from several domains
including online social networks (OSNs), web graph, co-authorship, citation,
etc. The network sizes range from 14 thousand to 2 billion edges. For each
network graph, self-loops were removed and directionally was ignored in
digraphs. The statistics of the graphs are listed in Table I.

###  V-B Verification of Theorem 1

To verify the approximations made in the derivation of the RSE of the
estimator in Theorem 1, the parameters of the estimator were set up to obtain
the RSEs between 0.1 and 0.4. First, we report the observed RSEs along with
its estimations based on Eq. 7 in Fig. 1. Only 4 representative graphs are
reported for lack of space. Similar pattern are observed for other graphs. We
make several observations as follows.

  * •

The estimated RSEs based on Eq. 7 fit perfectly the observed RSEs not only for
large graphs when $p$ is small but also for small ones. Take Ego-facbook, for
example, the smallest graph in our datasets, the observed RSEs (red circles)
fit very well our estimations based on Eq. 7.

  * •

When $p$ is small the first term in Eq. 7, i.e. $1/\Delta_{g}$, is a dominant
term compared to the remaining three. The reason is as follows. Firstly, when
$p$ is small, the probability of identifying shared triangles, i.e.
$\Phi_{g}$, and pairs of triangle and wedge with one common edge
($\Omega^{\prime}_{g}$) based on the sampled $g$ is very small. Note that
$\mathbb{E}(\Phi_{g})=\frac{8}{15}\Phi p^{3}$, and
$\mathbb{E}(\Omega^{\prime}_{g})=\frac{5}{12}\Omega^{\prime}p^{2}$. Moreover,
the term $2\Psi_{g}/\Lambda^{2}_{g}$ is ignorable due to the fact that their
expectations depends on the sampling probability $p$ which it is in the same
order in the both expectations. Therefore, in the sample
$\Psi_{g}<<\Lambda_{g}^{2}$.

  * •

By increasing sampling probability $p$, terms $+2\Phi_{g}/\Delta^{2}_{g}$ and
$-2\Omega^{\prime}_{g}/\Delta_{g}\Lambda_{g}$ are increasing in the same
order. Surprising that the two terms neutralize each other. Furthermore, term
$2\Psi_{g}/\Lambda^{2}_{g}$ remains ignorable compared to the other terms.

The observations above support our claim in Theorem 1 to simplify the RSE of
the estimator. To verify the result in Theorem 1, we report the observed RSEs
and their estimations based on Eq. 8 in Fig. 2. It can be seen that the
observed RSEs support our estimations not only for large graphs but also for
small ones. In a few small graphs, i.e. Web-Standford, Web-Berkeley, and As-
skitter, there are small gaps between the observed and estimated RSEs.
However, by increasing the size of networks (see the last row in the figure)
the estimated RSEs match perfectly the observed ones. Thus, we believe that
our result in Theorem 1 can be used in practice to determine the size of
samples to achieve a given accuracy level of the estimation.

###  V-C The bias

To understand the bias, we set the parameters of the estimator to achieve the
RSEs between 0.1 and 0.4. The estimator run on the graphs and the observed RBs
(relative bias) and the estimated ones based on Eq. 16 were computed. The
results were reported in Fig. 3 for the graphs with the RB more than 0.7%. We
make several observations as follows: 1) The observed RBs support the
estimated RBs based on Eq. 16 for all graphs; 2) Both negative and positive
biases are observed; 3) The largest RBs were observed on small graphs and it
can be as high as 2%; 4) In most of the graphs in our datasets, the bias is
very small and it is ignorable.

We also report the observed RBs of our biased estimator and the bias-corrected
one in Fig. 4. It can be seen that the bias of $\widehat{\mathcal{C}}^{+}$ was
removed in all the plots.

##  VI Discussions and Conclusions

This paper addresses the estimation of the bias and variance in a streaming
algorithm for estimating clustering coefficient. Essentially it is about
estimating the properties of an estimator. It is important since it is the
only way to know how good an estimate is during the estimation process.

Our result is obtained based on two simplifications. One is the Taylor
expansion–we take only the first two terms. The other is $1-p\approx
1-p^{2}\approx 1$, assuming that sampling probability is very small if the
data is very large. With such simplifications, we can characterize the
variance with a single variable, i.e., $\Delta_{g}$. Although in theory
variances depends on graph structures characterized by $\Phi$, $\Psi$, and
$\Omega^{\prime}$, all these variables can be neglected when estimations are
performed on very large graph. This simple yet powerful result is very useful
in practice– we can give a confidence interval when an estimate is given.

For the bias part, we conclude that it is small overall, and it is more
observable for smaller graphs. Interestingly, the bias can be either positive
or negative, depending on the graph structure.

Although our result is developed on NES algorithm, the same method can be
extended to numerous other streaming algorithms. Besides, NES itself is a very
powerful algorithm. Despite its simplicity, its performance is comparable to
the state-of-the-art algorithm GPS-In as illustrated in Fig.5. Hence, the
variance estimator and NES are a good combination to be used in practice.

##  VII Acknowledgments

The research is supported by NSERC and CZI (Chan Zuckerberg Initiative).

## References

  * [1] M. Latapy, “Main-memory triangle computations for very large (sparse (power-law)) graphs,” _Theoretical Computer Science_ , vol. 407, no. 1, pp. 458–473, 2008.
  * [2] N. K. Ahmed, N. Duffield, J. Neville, and R. Kompella, “Graph sample and hold: A framework for big-graph analytics,” in _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, $1446$-$1455$_. ACM, 2014.
  * [3] N. K. Ahmed, N. Duffield, T. L. Willke, and R. A. Rossi, “On sampling from massive graph streams,” _Proceedings of the VLDB Endowment_ , vol. 10, no. 11, pp. 1430–1441, 2017.
  * [4] M. Jha, C. Seshadhri, and A. Pinar, “A space efficient streaming algorithm for triangle counting using the birthday paradox,” in _Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , ser. KDD ’13. New York, USA: ACM, 2013, pp. 589–597.
  * [5] ——, “A space-efficient streaming algorithm for estimating transitivity and triangle counts using the birthday paradox,” _ACM Trans. Knowl. Discov. Data_ , vol. 9, no. 3, pp. 15:1–15:21, Feb. 2015.
  * [6] R. Etemadi and J. Lu, “Bias correction in clustering coefficient estimation,” in _2017 IEEE International Conference on Big Data (Big Data)_ , Dec 2017, pp. 606–615.
  * [7] T. Schank and D. Wagner, “Approximating clustering coefficient and transitivity,” _Journal of Graph Algorithms and Applications_ , vol. 9, no. 2, pp. 265–275, 2005.
  * [8] H. Jowhari and M. Ghodsi, “New streaming algorithms for counting triangles in graphs,” in _International Computing and Combinatorics Conference, $710$-$716$_. Springer, 2005.
  * [9] L. S. Buriol, G. Frahling, S. Leonardi, A. Marchetti-Spaccamela, and C. Sohler, “Counting triangles in data streams,” in _Proceedings of the twenty-fifth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database systems, $253$-$262$_. ACM, 2006.
  * [10] Y. Lim and U. Kang, “Mascot: Memory-efficient and accurate sampling for counting local triangles in graph streams,” in _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, $685$-$694$_. ACM, 2015.
  * [11] R. Etemadi, J. Lu, and Y. H. Tsin, “Efficient estimation of triangles in very large graphs,” in _Proceedings of the 25th ACM International on Conference on Information and Knowledge Management_ , ser. CIKM ’16. ACM, 2016, pp. 1251–1260.
  * [12] L. De Stefani, A. Epasto, M. Riondato, and E. Upfal, “$triest$: Counting local and global triangles in fully-dynamic streams with fixed memory size,” in _Proceedings of the 22th ACM KDD international conference on Knowledge discovery and data mining_. ACM, 2016.
  * [13] J. Leskovec and A. Krevl, “SNAP Datasets: Stanford large network dataset collection,” http://snap.stanford.edu/data, Jun. 2014.
  * [14] J. Kunegis, “Konect - the koblenz network collection,” http://konect.uni-koblenz.de/networks, May 2016.
  * [15] P. Boldi, M. Rosa, M. Santini, and S. Vigna, “Layered label propagation: A multiresolution coordinate-free ordering for compressing social networks,” in _Proceedings of the 20th international conference on WWW, $587$-$596$_. ACM, 2011.
  * [16] P. Boldi and S. Vigna, “The WebGraph framework I: Compression techniques,” in _Proceeding of the Thirteenth International World Wide Web Conference, $595$-$601$_. ACM, 2004\.
  * [17] A. Sinha, Z. Shen, Y. Song, H. Ma, D. Eide, B.-j. P. Hsu, and K. Wang, “An overview of microsoft academic service (mas) and applications,” in _Proceedings of the 24th International Conference on World Wide Web_. ACM, 2015, pp. 243–246.

## Appendix

######  Lemma 1.

Let $\Delta_{g}$ be the number of closed wedges identified based on $g$ using
NES. The variance of $\Delta_{g}$ is

$var(\Delta_{g})=\frac{1}{3}\left(\Delta(p^{2}-\frac{1}{3}p^{4})+8\Phi(\frac{2}{5}p^{3}-\frac{1}{3}p^{4})\right).$
(18)

###### Proof.

Applying the variance we have

$\displaystyle
var(\Delta_{g})=var(\sum\limits_{i=1}^{\Delta}\delta_{i})=\sum\limits_{i=1}^{\Delta}\sum\limits_{j=1}^{\Delta}cov(\delta_{i},\delta_{j})$
$\displaystyle=\sum\limits_{i=1}^{\Delta}var(\delta_{i})+\sum\limits_{i\neq
j}cov(\delta_{i},\delta_{j}).$

We remind the reader that indicator $\delta_{i}$ is 1 if $i^{th}$ closed wedge
is identified based on $g$; otherwise it is 0. The probability of identifying
a closed wedge by NES is $\frac{1}{3}p^{2}$. Hence, the variance of
$\delta_{i}$ is $\frac{1}{3}p^{2}-\frac{1}{9}p^{4}$. Therefore the cost of the
variance term is $\Delta(\frac{1}{3}p^{2}-\frac{1}{9}p^{4})$. The covariance
between two independent closed wedges is zero. Thus, we need to find the
covariance between dependent closed wedges. The probability of identifying
such a dependent case is $\frac{2}{15}p^{3}$. Hence, the covariance between
two shared closed wedges in such cases is
$\frac{2}{15}p^{3}-\frac{1}{9}p^{4}$. Recall that the total number of pairs of
shared triangles is denoted by $\Phi$. For each pair of shared triangles there
are four dependent closed wedges. Thus, in summation above, there are $8\Phi$
pairs of shared closed wedges. Therefore, the cost of covariance term is
$8\Phi(\frac{2}{15}p^{3}-\frac{1}{9}p^{4})$. Add the costs of the two terms,
we obtain the Lemma.

######  Lemma 2.

Suppose $\Lambda_{g}$ be the number of wedges identified based on $g$ using
NES. The variance of $\Lambda_{g}$ is obtained by

$\displaystyle var(\Lambda_{g})$
$\displaystyle=\Lambda(p-p^{2})+\frac{2}{3}\Psi(p-p^{2}).$ (19)

where $\Lambda$ and $\Psi$ are the number of wedges and the count of pairs of
shared wedges in $G$ respectively.

###### Proof.

Recall that $\lambda_{i}$ is an indicator for $i^{th}$ wedge in $G$. Indicator
$\lambda_{i}$ is 1 when $i^{th}$ wedge is identified based on $g$, and 0 else.
The variance of $\Lambda_{g}$ is

$\displaystyle
var(\Lambda_{g})=var(\sum\limits_{i=1}^{\Lambda}\lambda_{i})=\sum\limits_{i=1}^{\Lambda}\sum\limits_{j=1}^{\Lambda}cov(\lambda_{i},\lambda_{j})$
$\displaystyle=\sum\limits_{i=1}^{\Lambda}var(\lambda_{i})+\sum\limits_{i\neq
j}cov(\lambda_{i},\lambda_{j}).$

The probability of identifying a wedge by NES is $p$. When $i=j$ holds, the
covariance term is equal to the variance of $\lambda_{i}$, which is $p-p^{2}$.
Thus, we get $var(\widehat{\Lambda})=\Lambda(p-p^{2})+\sum\limits_{i\neq
j}cov(w_{i},w_{j}).$ Next, we need to understand the covariance between two
wedges. When the two wedges $w_{i}$ and $w_{j}$ are independent the covariance
between them is zero. Hence, we need to consider the covariance between
dependent wedges. The probability of identifying a pair of dependent wedges is
$p$. Thus, the covariance between two shared closed wedges is $(p-p^{2})$.
Suppose $\Psi$ be the number of pairs of dependent wedges in $G$. The chance
to identify such a dependent case is $1/3$. Because
$cov(\lambda_{i},\lambda_{j})=cov(\lambda_{j},\lambda_{i})$ , we need to
multiply the covariance term by 2. Thus, the lemma is proved.

######  Lemma 3.

Covariance between $\Delta_{g}$ and $\Lambda_{g}$ is given by

$cov(\Delta_{g},\Lambda_{g})=2\Delta(p^{2}-p^{3})+\frac{5}{12}\Omega^{\prime}(p^{2}-p^{3}).$
(20)

where $\Omega^{\prime}$ is the number of pairs of wedges and triangles with
one common edge.

###### Proof.

Let indicators $\delta_{i}$ and $\lambda_{j}$ be the same as defined before.
The covariance is given by

$cov(\Delta_{g},\Lambda_{g})=\sum\limits_{i=1}^{\Delta}\sum\limits_{j=1}^{\Lambda}cov(\delta_{i},\lambda_{j}).$

Suppose $\Omega^{\prime}$ be the exact number of pairs of wedges and triangles
with one common edge in $G$. The sampling probabilities of such a pair is
$p^{2}$. Moreover, the chance to identify such a pair in the streaming model
by NES is $5/12$. Therefore, the total cost is
$\frac{5}{12}\Omega^{\prime}(p^{2}-p^{3})$. In addition, for each closed wedge
of a triangle, there are two cases of pairs of wedges and closed ones with a
common edge. The total cost for such cases is
$6\Delta(\frac{1}{3}p^{2}+\frac{1}{3}p^{3})$. Therefore, we get the lemma.