#  Self- similar behavior of the Neutron Fracture Functions

Samira Shoeibi Mohsenabadi1 ID<EMAIL_ADDRESS>F.
Taghavi-Shahri1 ID<EMAIL_ADDRESS>(Corresponing author) (1)Department
of Physics, Ferdowsi University of Mashhad, P.O.Box 1436, Mashhad, Iran

###### Abstract

In this article, we have employed fractal formalism to calculate the Fracture
Functions of the Leading neutron produced in ep collisions. The fractal
concept describes the self-similar behavior of the proton structure at Leading
neutron production of semi-inclusive Deep Inelastic Scattering at the low
values of the fractional momentum variable $\beta$. The Fracture Functions
(FFs) parameterized the non-perturbative part of the fragmentation process at
the initial scale of Q${}_{0}^{2}$. In this analysis, we benefit from the
Leading neutron (Ln) experimental data published by H1 Collaboration and in
order to estimate the uncertainty of neutron FFs and corresponding
observables, we used the Hessian method. As a consequence of this Next-to-
Leading order QCD analysis, we achieve a nice agreement between the prediction
of our model for neutron FFs and experimental data. It seems that the fractal
approach based on the self-similar behavior of these conditional parton
distribution functions at low x can nicely describe the experimental data.

###### pacs:

12.38.Bx, 12.39.-x, 14.65.Bt

###### Contents

  1. I Introduction
  2. II Theory setup
    1. II.1 Fracture Functions and the Leading neutron Production
    2. II.2 Fractals
    3. II.3 The Validity test of using the Fractal Model In Leading neutron Production Mechanism
  3. III Essential parts of the Fractal Neutron FFs and their uncertainties
    1. III.1 Neutron FFs in the Fractal approach
    2. III.2 Experimental data
    3. III.3 $\chi^{2}$ minimization and nFFs Uncertainties
  4. IV Results and discussion
  5. V Summary and Conclusions

##  I Introduction

Proton is known as a strongly bound state of quarks and gluons in Quantum
Chromodynamics (QCD). Considering its collisions with electron not only
clarify its structure, but also reveal the dynamics of strong interaction at
high energy physics Bourrely:2018qit ; AbdulKhalek:2019bux ; Zaidi:2019mfd ;
Harland-Lang:2014zoa ; Ball:2014uwa ; Martin:2009iq ; Alekhin:2018pai ;
Salajegheh:2019zos . Investigation of the Parton Distribution Functions (PDFs)
and the proton structure functions at small values of x (the Bjorken scaling
variable, which represents the proton momentum fraction carried by the struck
quark) is interesting and the subject of some debates during the last decades
Heidari:2019fio ; Abt:2016vjh ; Abdolmaleki:2018jln . One of the fascinating
phenomenological models that well describes experimental data at low x is
based on the scaling behavior or fractal behavior of PDFs Dremin:1988mu ;
Lipa:1989yh ; Florkowski:1990ba ; Bjorken:1994uf ; Bhattacharya:2000uh .
Fractal geometry provides a new perspective to describe the natural phenomena
and they are based on self-similarity symmetry with respect to scale or size:
A small part of the system is approximately similar to the entire system.
Benoit Mandelbrot, as a pioneer scientist, defined the fractal as a set of
self-similar objects that characterized by the fractional dimension (Hausdorff
dimension) which is different than the topological one Mandelbrot . The low x
behavior of the proton structure is dominated by the gluon-gluon interactions
which increase the gluon and the sea quark densities, leading to the self-
similarity in PDFs. Recently, relevant work has been done in the investigation
of the Deep Inelastic Scattering (DIS) processes Dremin:1992zc . The idea of
considering the self-similar features in the multi partons interactions was
proposed by Lastovcka Lastovicka:2002hw ; Lastovicka:2004mq  and others have
used it to investigate the proton structure functions Choudhury:2003yy ;
Choudhury:2005vy ; Jahan:2011ig ; Choudhury:2013ita ; Jahan:2014ova ;
Jahan:2014sqa ; Choudhury:2016fjy ; Deppman:2016prl . In these references, the
fractal formalism are applied to proton structure functions up to Leading
Order (LO) and they do not consider the gluon densities on their analysis.

Here in this paper, we intend to use the fractal model to describe the Leading
neutron production processes and calculate the neutron Fracture Functions
(nFFs) in the positron-proton collisions. These Leading neutrons are high-
energy particles that carry a remarkable portion of the incoming proton
energy. They are also produced at a small polar angle with respect to the
positron-proton collision axis. The FFs are conditional probabilities of
finding partons with flavor i and momentum fraction x while simultaneously
hadron h is detected in the final states of the interaction between lepton and
proton. These FFs obey evolution equations with an additional inhomogeneous
term as compared to the standard Altarilli-Parisi equations deFlorian:1997wi .
They also describe the semi-inclusive DIS processes in the target
fragmentation region deFlorian:2017lwf ; Bertone:2017tyb ; Anderle:2015lqa ;
Ethier:2017zbq ; Altarelli:1979kv ; deFlorian:1997wi ; Daleo:2003jf ;
Trentadue:1993ka ; Trentadue:1994iw . Although they play a fundamental role in
describing the semi-inclusive DIS processes and the global analysis requires
their shapes at a given initial scale, QCD is unable to provide or predict
them yet. Therefore, it is not only crucial to parametrize them
phenomenologically, but also it is necessary to use experimental data in
parallel, to make a connection with perturbative QCD.
From the experimental point of view, the ZEUS and H1 collaborations published
the experimental Leading baryons data Adloff:1998yg ; Aaron:2010ab ;
Chekanov:2008tn ; Chekanov:2002yh ; Chekanov:2002pf ; Rinaldi:2006mf . The
semi-inclusive cross-sections for the production of Leading baryons can be
studied by using FFs concept Shoeibi:2017lrl ; Shoeibi:2017zha ;
Ceccopieri:2014rpa .
In the following, we will investigate the cross-section of the Leading neutron
production by adopting the fractal model and defining the neutron Fracture
Functions, ”nFFs”. In Sec. II, fractals are introduced and the FFs in Leading
neutron production are discussed. Furthermore, the emergence of the fractal
features in Leading neutron production is examined. In Sec. III, the
parametrization of the neutron FFs are provided within the fractal model and
we aso compared them with experimental data. The uncertainties of nFF are
given, using the Hessian method. In Sec. IV, our findings are checked against
H1 experimental data. In Sec. V, we finish with the summary and our
conclusions.

##  II Theory setup

In this section, we will review the theoretical aspects of the QCD analysis
performed in this paper. To begin, let us have a short review on the nFFs
concept in QCD. They are also reviewed in Refs. Ceccopieri:2014rpa ;
Shoeibi:2017lrl ; Shoeibi:2017zha . We then move to the discussion of the
fractal concept and its footprint in the Leading neutron production mechanism.

###  II.1 Fracture Functions and the Leading neutron Production

Inclusive DIS processes are described by three kinematic variables:
$x=\frac{Q^{2}}{2p.q}$, $Q^{2}=-q^{2}$ and
$y=\frac{p.q}{p.k}\simeq\frac{Q^{2}}{sx}$ where $p$, $k$, and $q$ are the
four-momenta of the incident proton, the incident positron and the exchanged
virtual photons, respectively. Here $\sqrt{s}$ is the center of mass energy
for positron and proton. However, semi-inclusive DIS describes the mechanism
of the Leading neutrons produced in the final state by using two more
kinematic variables Altarelli:1979kv :

$\displaystyle x_{L}$ $\displaystyle=$ $\displaystyle
1-\frac{q.(p-p_{n})}{q.p}\simeq\frac{E_{n}}{E_{p}}$ $\displaystyle t$
$\displaystyle=$
$\displaystyle(p-p_{n})^{2}\simeq-\frac{p_{T}^{2}}{x_{L}}-(1-x_{L})(\frac{m_{n}^{2}}{x_{L}}-m_{p}^{2})$
(1)

where $x_{L}$ is the momentum fraction carried by the outgoing neutron and $t$
is squared four-momentum transfer between the incoming proton and the final
state neutron. In Eq. (II.1), $m_{p}$ is the proton mass and $p_{n}$ is the
four-momentum of the final state neutron. The outgoing neutron energy is
$E_{n}$ and its transverse momentum for the incoming proton is denoted by
$p_{T}$.

The four-fold differential cross-section for the Leading neutron production,
$e^{+}p\rightarrow e^{+}nX$ processes, can be written as a function of semi-
inclusive Leading neutron transverse and longitudinal structure functions,
$F_{2}^{\rm Ln(4)}$ and $F_{L}^{\rm Ln(4)}$, which is defined as Aaron:2010ab
; Chekanov:2002pf ; deFlorian:1998rj ; deFlorian:1997wi ; Chekanov:2002yh :

$\displaystyle\frac{d^{4}\sigma(ep\to
e^{\prime}nX)}{d\beta\,dQ^{2}\;dx_{L}\,dt}$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}}{\beta
Q^{4}}(1-y+\frac{y^{2}}{2})F_{2}^{\rm Ln(4)}(\beta,Q^{2};x_{L},t)$ (2)
$\displaystyle+$ $\displaystyle F_{L}^{\rm Ln(4)}(\beta,Q^{2};x_{L},t)\,.$

The scaled fractional momentum variable $\beta$ is defined by:

$\beta=\frac{x}{1-x_{L}}$ (3)

In Eq. (2), $F_{2,L}^{Ln(4)}(\beta,Q^{2};x_{L},p^{2}_{T})$ can also be
interpreted as the sub-structure function of the proton structure functions,
$F_{2,L}(x,Q^{2})$ in the inclusive DIS. The $t$ integrated differential cross
section can be obtained by:

$\displaystyle\frac{d^{3}\sigma(ep\to e^{\prime}nX)}{d\beta\,dQ^{2}\,dx_{L}}$
$\displaystyle=$ $\displaystyle\int_{t_{0}}^{t_{min}}\frac{d^{4}\sigma(ep\to
e^{\prime}nX)}{d\beta\,dQ^{2}\ dx_{L}\,dt}dt$ (4) $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}}{\beta
Q^{4}}(1-y+\frac{y^{2}}{2})F_{2}^{\rm Ln(3)}(\beta,Q^{2};x_{L})$
$\displaystyle+$ $\displaystyle F_{L}^{\rm Ln(3)}(\beta,Q^{2};x_{L})\,,$

where the integration limits are:

$\displaystyle t_{min}$ $\displaystyle=$
$\displaystyle-(1-x_{L})(\frac{m_{N}^{2}}{x_{L}}-m_{p}^{2})\,,$ $\displaystyle
t_{0}$ $\displaystyle=$ $\displaystyle
t_{min}-\frac{(p_{T}^{max})^{2}}{x_{L}}\,.$ (5)

$p_{T}^{max}$ is the upper limit of the neutron transverse momentum used for
the $F_{2,L}^{\rm Ln(3)}$ measurement. In this paper we define the reduced
$e^{+}p$ cross section $\sigma_{r}^{\rm Ln(3)}$ in term of leading neutron
transverse $F_{2}^{\rm Ln(3)}$ and the longitudinal structure functions
$F_{L}^{\rm Ln(3)}$ as Aaron:2010ab ; Chekanov:2002pf :

$\displaystyle\sigma_{r}^{Ln(3)}$ $\displaystyle=$ $\displaystyle
F_{2}^{Ln(3)}(\beta,Q^{2};x_{L})$ (6) $\displaystyle-$
$\displaystyle\frac{y^{2}}{1+(1-y)^{2}}F_{L}^{Ln(3)}(\beta,Q^{2};x_{L})\,.$

By applying the factorization theorem to the semi-inclusive DIS, we can write
the semi-inclusive Leading neutron transverse and longitudinal structure
functions, $F_{2}^{\rm Ln(4)}$ and $F_{L}^{\rm Ln(4)}$,
$F_{2}^{Ln(4)}(\beta,Q^{2};x_{L},p^{2}_{T})$ as the convolution of conditional
parton distribution functions, ${\cal
M}^{n}_{i/p}(\beta,\mu_{F}^{2};x_{L},p_{T}^{2})$: ”proton to neutron Fracture
Functions” and Wilson coefficient functions, $C_{ki}$ (k = 2, L)
Ceccopieri:2014rpa :

$\displaystyle
F_{k}^{Ln(4)}(\beta,Q^{2};x_{L},p_{T}^{2})=\sum_{i}\int_{\beta}^{1}\frac{d\xi}{\xi}{\cal
M}^{n}_{i/p}(\beta,\mu_{F}^{2};x_{L},p_{T}^{2})\times
C_{ki}(\frac{\beta}{\xi},\frac{Q^{2}}{\mu_{F}^{2}},\alpha_{s}(\mu^{2}_{R}))+{\cal
O}(\frac{1}{Q^{2}})\,.$ (7)

here $\mu_{F}^{2}$ and $\mu^{2}_{R}$ are factorization and renormalization
scales, respectively and ${\cal
M}^{n}_{i/p}(\beta,\mu_{F}^{2};x_{L},p_{T}^{2})$ describes the probability of
detecting Leading neutron with momentum fraction $x_{L}$ and transverse
momentum $p_{T}^{2}$ in the final state, while simultaneously the specific
parton i with fraction momentum of $\beta$ interacts with the virtual photon.
Similar to the usual parton densities, these functions are independent of the
hard scattering process and represent the objective structure of the target.
They also contain rich information about the proton structure functions and
non-perturbative dynamics. Integrating
$F_{2}^{Ln(4)}(\beta,Q^{2};x_{L},p_{T}^{2})$ with respect to $p_{T}^{2}$, up
to the small value of $p_{T}^{max}\ll Q^{2}$, enhances the relative
contribution of the exchanged object. Therefore, one may interpret nFFs,
${\cal M}^{n}_{i/p}(\beta,\mu_{F}^{2};x_{L})$, as the parton distribution of
the exchanged object (normally pion) between the target and the final state
Leading neutron. Using limiting fragmentation hypothesis Benecke:1969sh ,
$F_{2}^{Ln(4)}(\beta,Q^{2};x_{L},p_{T}^{2})$ can be written as the pion flux
factor, $f$, (which depends on Leading neutron variable ($x_{L},p_{T}^{2}$))
and the pion structure function, $F_{2}^{\pi}$, (which depends on the lepton
variable ($\beta,Q^{2}$)) Chekanov:2002pf :

$F_{2}^{Ln(4)}(\beta,Q^{2};x_{L},p_{T}^{2})=f(x_{L},p_{T}^{2})F_{2}^{\pi}(\beta,Q^{2})$
(8)

Recently different groups Arash:2003js ; McKenney:2015xis ; Barry:2018ort
have used this equation in order to extract the pionic distribution functions
by using the experimantal data related to the Leading neutron production
Chekanov:2002pf ; Aaron:2010ab . It is worth mentioning that for the Leading
neutron production process, the pion exchange mechanism is usually considered
tobe dominated at large neutron longitudinal momentum fractions $x_{L}$ and
nonpionic contributions become large at small $x_{L}$. The sensitivity to the
one-pion exchange contribution is considered in detail at Ref McKenney:2015xis
. The scale dependence of the neutron FFs can be extracted from Leading
neutron production processes. Integrating the proton-to-neutron FFs over
$p^{2}_{T}$ up to the small value of $p^{2}_{T_{max}}=\epsilon Q^{2}$
($\epsilon<1$), in fixed value of $x_{L}$, obey the standard DGLAP evolution
equations Ceccopieri:2014rpa ; Ceccopieri:2007th ; Shoeibi:2017lrl :

$Q^{2}\frac{\partial M^{B}_{i/p}(\beta,Q^{2};x_{L})}{\partial
Q^{2}}=\frac{\alpha_{s}(Q^{2})}{2\pi}\int_{\beta}^{1}\frac{du}{u}P_{i}^{j}(u)M^{B}_{j/p}(\frac{\beta}{u},Q^{2};x_{L})$
(9)

where $P_{i}^{j}(u)$ are regularized Altarelli-Parisi splitting functions
Altarelli:1977zs . Although the Wilson coefficient functions are the same as
those obtained in inclusive DIS and they are calculable within pQCD
Vermaseren:2005qc , the nFFs have non-perturbative nature and their
parameterized form should be inspired by and obtained from experimental data.

###  II.2 Fractals

Fractal geometry is a new language used to describe the models and analyze the
complex forms found in the nature. Plants, clouds, weather, turbulence,
chaotic phenomena in fluids, biological phenomena such as biological time
series or the growth pattern of bacteria and many other phenomena in nature
can be modeled by the fractals. Fractals have three general properties: a)
Self-similarity. b) Iterative formation. c) Fractional or fractal dimension.
In mathematical language, similarity means that when the size is varied, the
form would not be changed. Fig. (1) shows two samples of exact self-similar
fractals. By iteration, they are made of smaller and smaller shapes. In this
figure, M is the number of each side segment named as ”magnification factor”
and N is the number of self-similar shapes created by the segmentation. The
self-similarity dimension is introduced in Ref. Mandelbrot  and used to
consider the proton structure function at low values of x in Refs.
Lastovicka:2002hw ; Lastovicka:2004mq . It defines as:

$D=\frac{log(number~{}of~{}self-
similar~{}objects)}{log(magnification~{}factor)}=\frac{log~{}M^{D}}{log~{}M}$
(10)

Figure 1:  (Color online) Equal segmentation of square and Sierpinski triangle
fractal. M is the number of each side segment named it as magnification factor
and N is the number of self-similar shapes created by the segmentation.

Using Eq. (10) for square and Sierpinski triangle fractal in Fig. (1), we
achieve the dimension values of these figures as 2 and 1.5849, respectively.
The fractional or fractal dimension for the Sierpinski triangle fractal
indicates the degree of detail in this object and measures its complexity. The
fractal dimension also demonstrates how much space it occupies between the
euclidean dimensions.

Let’s look at fractal dimension in more details. Eq. (10) is the definition of
a fractal dimension for discrete fractals. In general, $\zeta$ is a real
number and the number of self-similar objects is represented by f($\zeta$).
Since the dimension may change with scaling, a local dimension is defined as:

$D_{f}(\zeta)=\frac{\partial log(f(\zeta))}{\partial log(\zeta)}$ (11)

For ideal mathematical fractals, mentioned above, $D_{f}(\zeta)$ is constant
for the whole fractal ($D_{f}(\zeta)$$\equiv$D). However, there are many
fractals in nature that they are not mathematically ideal and most of the time
they have mono fractal structure only in a certain ragion of magnification.as
a result, introducing a scale dependent dimension is not weird. If we consider
a region on which the the dimension is approximately constant, we can rewrite
the Eq.(11) as:

$D=\frac{\partial log~{}f(\zeta)}{\partial log~{}\zeta}$ (12)

Therefore, the density function can be written as:

$log~{}f(\zeta)=D~{}.~{}log~{}\zeta+D_{0}$ (13)

In this relation, D0 defines the normalization factor, and $f(\zeta)$ is a
power-law function ($f(\zeta)\propto\zeta^{D}$) where $D$ is the fractal
dimension of this density function. The linear behavior of $logf(\zeta)$ as
$log\zeta$ has a key role to find the magnification factor of the density
function. For two independent magnification factors, $\zeta$ and $\eta$, this
relation can be extended as:

$log~{}f(\zeta,\eta)=D_{\zeta}~{}.~{}log~{}\zeta+D_{\eta}~{}.~{}log~{}\eta+D_{\zeta\eta}~{}.~{}log~{}\eta~{}log~{}\zeta+D_{0}$
(14)

$D_{\zeta\eta}$ demonstrates the dimensional correlation related to two
magnification factors, $\zeta$ and $\eta$.

###  II.3 The Validity test of using the Fractal Model In Leading neutron
Production Mechanism

It is known that at low x region, scale violation occurs and the parton
densities rise. If we assume this region obey the fractal features, the
unintegrated parton densities related to the Leading neutron production can be
described by the fractal model mentioned in section II.2. In order to check
whether the fractal model can be useful in describing the neutron production
mechanism or not, we employ the Ref. Shoeibi:2017lrl  parameterization of
nFFs. This enables us not only to check the validity of the model but also
provides us with the magnification factor(s) of the unintegrated nFFs. The
nFFs, ${\cal M}^{n}_{i/p}(\beta,Q^{2};x_{L})$ are defined as follows:

${\cal M}^{n}_{i/p}(\beta,Q^{2};x_{L})=\int_{0}^{Q^{2}}{\cal
M}^{n}_{i/p}(\beta,k_{t}^{2};x_{L})~{}dk_{t}^{2}$ (15)

where $k_{t}^{2}$ is the parton transverse momentum. Since unintegrated nFFs,
${\cal M}^{n}_{i/p}(\beta,k_{t}^{2};x_{L})$, depending on
$\beta=\dfrac{x}{1-x_{L}}$, $x_{L}$ and $k_{t}^{2}$ therefore, it seems that
there are three candidates for magnification factors of unintegrated nFFs:
$\beta$ (or x), $x_{L}$ and $k_{t}^{2}$. However, accourding to factorization
hypotheses mentioned in Refs. H1:1998vul , the leading neutron transverse
structure functions follow the general forms of

$F_{k}^{Ln(3)}(\beta,Q^{2};x_{L})=f(x_{L})~{}F_{k}^{Ln(2)}(\beta,Q^{2})$ (16)

where the discrete-function $f(x_{L})$ is known as chiral splitting function
McKenney:2015xis  and it describes the splitting of a proton into a $\pi n$
system. In this analysis this function is expressed as five free parameters.
The parametrized form for F${}_{2}^{Ln(2)}(\beta,Q^{2})$ is then based on the
proton structure functions. Therefore, the two variabe of $\beta$ (or x) and
$k_{t}^{2}$ are the final candidates for magnification factors of unintegrated
nFFs. The behaviour of logarithmic unintegrated parton density in order of
$log(x)$ in fixed value of $Q^{2}=15~{}GeV^{2}$ , $x_{L}=0.5$ and its
logarithmic behaviour in order of $log(Q^{2})$ in fixed value of $x=0.001$ ,
$x_{L}=0.5$ are shown in Fig. (2). The linear behavior is detected for
$x<0.001$ and $Q^{2}>7GeV^{2}$. If we consider the unintegrated densities as
fractal, in these two regions, we can find the constant value of fractal
dimension by using Eq. (12). This region known as mono fractal region, while
in the other regions we expect the multi fractal behavior for those
unintegrated PDFs. It means that in Eq. (11), the fractal dimension should be
local one.

As mentioned in Refs. Lastovicka:2002hw ; Lastovicka:2004mq , magnification
factors have some properties, namely, they should be positive, nonzero and
dimensionless. As a result, in order to satisfy the two latter features, we
choose $1+\frac{Q^{2}}{q_{0}^{2}}$ instead of $Q^{2}$. In addition, when the
structure is probed deeper, x goes to zero while the magnification factor
should increase. Therefore, we choose $\frac{1}{x}$ instead of x.
Consequently, we choose $1+\frac{Q^{2}}{q_{0}^{2}}$ and $\frac{1}{x}$ as
magnification factors of unintegrated parton densities in Leading neutron
production mechanism.

Figure 2:  (Color online) The validity test of the fractal model for Leading
neutron production mechanism by using nFFs proposed in Ref Shoeibi:2017lrl

##  III Essential parts of the Fractal Neutron FFs and their uncertainties

The concept of the Fracture Functions, FFs, is useful to describe and
understand the physics in the forward region. Since these universal functions
cannot be calculated directly in pQCD, one needs to assume a functional form
for them with a number of free parameters and then, attempt to fix the
parameters by using the experimental data. Here, we introduce the
parameterized form of the nFFs at an initial $Q_{0}^{2}$ based on the fractal
model. We also briefly review the experimental data and estimation of the
uncertainties.

###  III.1 Neutron FFs in the Fractal approach

Following the sections II.2 and II.3 and Refs. Lastovicka:2002hw ;
Lastovicka:2004mq , we introduce the unintegrated nFFs as follows:

$log({\cal
M}^{n}_{i/p}(\beta,k_{t}^{2};x_{L}))=\theta(x_{L})[D_{1}^{i}log(\frac{1}{\beta})log(1+\frac{{k_{t}^{2}}}{q_{0}^{2}})+D_{2}^{i}log(\frac{1}{\beta})+D_{3}log(1+\frac{{k_{t}^{2}}}{q_{0}^{2}})+D_{0}^{i}-log(M^{2})]$
(17)

where q0 is a free parameter to be determined in the analysis and
$\theta(x_{L})$ is a $x_{L}$-dependent function. The flavour dependent
parameters in Eq. (17) are D${}_{0}^{i}$, D${}_{1}^{i}$ and D${}^{i}_{2}$.
D${}_{0}^{i}$ fixes the normalization of ${\cal
M}^{n}_{i/p}(\beta,k_{t}^{2};x_{L})$. D${}_{1}^{i}$ is the dimensional
correlation relating to the $\beta$ and Q2 factors. D${}_{2}^{i}$ and D3 are
the dimensions related to the magnification factors $\beta$ and Q2,
respectively. The superscript i denotes the gluon and light quarks (q=u,
$\bar{u}$, d, $\bar{d}$, s, and $\bar{s}$). Using the parameterized form given
in Eq. (17) and Eq. (15), the initial integrated nFFs (up to Q${}_{0}^{2}$
value) are written as:

$\displaystyle\beta{\cal M}^{n}_{q/p}(\beta,Q_{0}^{2};x_{L})$ $\displaystyle=$
$\displaystyle
g(x_{L}).\frac{e^{D_{0}^{q}}q_{0}^{2}\beta^{-D_{2}^{q}+1}}{(M^{2})(1+D_{3}-D_{1}^{q}log(\beta))}(\beta^{-D_{1}^{q}log(1+\dfrac{Q_{0}^{2}}{q_{0}^{2}})}(1+\dfrac{Q_{0}^{2}}{q_{0}^{2}})^{D_{3}+1}-1)$
$\displaystyle\beta{\cal M}^{n}_{g/p}(\beta,Q_{0}^{2};x_{L})$ $\displaystyle=$
$\displaystyle
g(x_{L}).\frac{e^{D_{0}^{g}}q_{0}^{2}\beta^{-D_{2}^{g}+1}}{(M^{2})(1+D_{3}-D_{1}^{g}log(\beta))}(\beta^{-D_{1}^{g}log(1+\dfrac{Q_{0}^{2}}{q_{0}^{2}})}(1+\dfrac{Q_{0}^{2}}{q_{0}^{2}})^{D_{3}+1}-1)$
(18)

here $g(x_{L})$ is defined as:

$g(x_{L})=N_{i}x_{L}^{A_{i}}(1-x_{L})^{B_{i}}(1+C_{i}x_{L}^{D_{i}})$ (19)

Note that the additional term, -logM2, in Eq. (17), having dimension of energy
squared, so the integrated parton distribution function defined in Eq. ( 15)
is dimensionless. We set M2=1 GeV2 Jahan:2014sqa .

The power-law behavior of the nFFs can be considered in different point of
view: Accourding to Eq. (15), at low value of $\beta$, the unintegrated gluon
distribution g($\beta$, $k_{t}^{2}$; $x_{L}$) is defined as:

$xg(\beta,Q^{2};x_{L})=\int_{0}^{Q^{2}}g(\beta,k_{t}^{2};x_{L})~{}dk_{t}^{2}$
(20)

g($\beta$, $k_{t}^{2}$; $x_{L}$) satisfy the BFKL evolution equation:

$\dfrac{dg(\beta,k_{t}^{2};x_{L})}{dln(1/x)}=K_{L}\otimes g=\lambda g$ (21)

where $\otimes$ stands for convolution and $K_{L}$ is the Lipatov kernel which
represent the sum over powers of $\alpha_{s}ln(1/x)$ terms Lastovicka:2004mq .
From Eq. (21) it is obvious that the function g follows a power-law behaviour
in the variable $\beta$. Therefore, the gluon distribution can be expressed
as:

$xg(\beta,Q^{2};x_{L})\propto f(Q^{2})\beta^{-\lambda}$ (22)

The special case for which DGLAP equations can be solved analytically is the
double leading log approximation (DLL) Ninomiya:1976qt , which is based on the
small x limit. The DLL solution represents as:

$g(x,Q^{2})~{}=~{}e^{2\sqrt{\alpha_{s}^{{}^{\prime}}.ln\dfrac{Q^{2}}{Q_{0}^{2}}.ln\dfrac{1}{x}}}.g(x,Q_{0}^{2})$
(23)

here $\alpha_{s}^{{}^{\prime}}=\dfrac{3\alpha_{s}}{\pi}$ and $g(x,Q_{0}^{2})$
is the gluon distribution at factorisation scale Q0 Ninomiya:1976qt .

In addition, The rise of $F_{2}(x,Q^{2})$ at low x can be studied using the
following structure function derivative H1:2001ert :

$\lambda(x,Q^{2})=-(\dfrac{\partial lnF_{2}(x,Q^{2})}{\partial lnx})_{Q^{2}}$
(24)

Accourding to the H1 publication H1:2000muc , $\lambda(x,Q^{2})$ does not
depend on x, for x ¡ 0.01. If we consider the fitting procedure of the form of
$F_{2}(x,Q^{2})=c(Q^{2})~{}.~{}x^{-\lambda(Q^{2})}$ to the H1 structure
function data H1:2000muc , we see that in the region of x¡0.01, there is a
nice agreement with the form of
$\lambda(Q^{2})=a~{}ln\dfrac{Q^{2}}{\Lambda^{2}}$ and H1 data. Now Let’s come
back to our input in Eq. (III.1). There are two terms in our proposed model,
that they explain nFFs behavior at small value of $\beta$ :
$\beta^{-D_{2}^{i}+1}$ and
$\beta^{-D_{1}^{i}log(1+\dfrac{Q_{0}^{2}}{q_{0}^{2}})}$. When ones compare the
latter term with the above information, they would find out that Eq. (III.1)
is not a weird form for gluon and sea quark densities. It is just like a
window that gives us a new perspective on these distribution functions.

As mentioned in Ref. Lastovicka:2002hw , the structure functions of proton is
the subsequence of quark density which depend on x and $Q_{0}^{2}$ (in this
paper $Q_{0}^{2}$=$Q^{2}$). In this reference the Quark-Parton model is used
to extract the proton structure function (which depend on x and $Q^{2}$) and
the gluon density as an input is not considered. Since nFFs are conditional
PDFs, we propose our fractal model by inspiring from Ref. Lastovicka:2002hw .
As outlined earlier, nFFs obey the DGLAP evolution equations at a fixed value
of $x_{L}$. Therefore, we set the value of $Q_{0}^{2}$=1 GeV2 in Eq. (III.1)
to be able to compare the result of this analysis with the previous one (which
are aslo in initial value of $Q_{0}^{2}$=1 GeV2) Ceccopieri:2014rpa ;
Shoeibi:2017lrl . Convolution of nFFs with the Wilson coefficients, $C_{i}$,
gives the Leading neutron transverse structure functions Ceccopieri:2014rpa ;
Shoeibi:2017lrl :

$\displaystyle
F_{k}^{Ln(3)}(\beta,Q^{2};x_{L})=\sum_{i}\int_{\beta}^{1}\frac{d\xi}{\xi}{\cal
M}^{n}_{i/p}(\beta,\mu_{F}^{2};x_{L})\times
C_{ki}(\frac{\beta}{\xi},\frac{Q^{2}}{\mu_{F}^{2}},\alpha_{s}(\mu^{2}_{R}))+{\cal
O}(\frac{1}{Q^{2}})\,.$ (25)

Since there are a few experimental data related to the Leading neutron
production, in this analysis we assume the symmetric light quark
distributions: ${\cal M}^{n}_{u/p}={\cal M}^{n}_{\bar{u}/p}={\cal
M}^{n}_{d/p}={\cal M}^{n}_{\bar{d}/p}={\cal M}^{n}_{s/p}={\cal
M}^{n}_{\bar{s}/p}$. So, at the initial scale of $Q^{2}_{0}=1~{}GeV^{2}$, we
have the singlet and gluon distributions as inputs. As fully describe in Ref.
Ceccopieri:2014rpa , $g(x_{L})$ can be described by common parameters of
Bg=B${}^{q}\equiv$B, Cg=C${}^{q}\equiv$C and Dg=D${}^{q}\equiv$D. Therefore
the Eq. 19 can be rewritten as:
$g(x_{L})=N_{i}x_{L}^{A_{i}}(1-x_{L})^{B}(1+Cx_{L}^{D})$. There are two
parameters that control the normalization of inputs: Ni and D${}^{i}_{0}$. We
considered these parameters as only one free parameter and we name it as:
$\Gamma_{i}$ = N${}_{i}e^{D_{0}^{i}}$. As a results, It turns out that the
model contains 13 free parameters to be determined using the experimental data
($\Gamma_{i}$, Ai, B, C, D, q0, D${}_{1}^{i}$, D${}_{2}^{i}$, and D3 ). To
perform the analysis, we used QCDNUM package Botje:2010ay  to not only solve
the DGLAP evolution equations for the integrated nFFs in Eq. (III.1), but also
we use it to extract the Leading neutron transverse structure functions,
$F_{2}^{Ln(3)}(\beta,Q^{2};x_{L})$ up to NLO QCD approximation. We did our
analysis in Variable Flavour Number Scheme (VFNS), which means that when the
$Q^{2}$ value is equal to the pole mass of the heavy quarks, c, b and t
(mc=1.4 GeV, mb=4.5 GeV and mt=175 GeV), the number of flavours change from
$n_{f}$ to $n_{f}$+1. In addition, we set the strong coupling constant on it’s
recently updated average, $\alpha_{s}(M_{z})=0.1177$ dEnterria:2015kmd .

###  III.2 Experimental data

The production of Leading neutrons in the semi-inclusive process $e^{+}$p
$\longrightarrow$ $e^{+}$nX are considered by the ZEUS Chekanov:2002pf  and H1
Aaron:2010ab  collaborations at HERA. The H1 experiment measured
F${}_{2}^{LN(3)}$ in the phase space defined as 6 $<Q^{2}<$100 GeV2, 1.5
*10${}^{-4}<x<3*10^{-2}$, the longitudinal momentum fraction 0.32$<x_{L}<$0.95
and neutron transverse momentum $p_{T}$¡200 MeV. In this interaction, the
proton, and positron beam energies are $E_{p}$=920 GeV and $E_{e}$=27.6 GeV,
respectively. Notice that the small value of the transverse momentum increases
the relative contribution of pion exchange Holtmann:1994rs ;
Kopeliovich:1996iw  and the maximum transverse momentum, $p_{T}^{max}$, is set
to 200 MeV Adloff:1998yg ; Aaron:2010ab . ZEUS data also cover the large range
of kinematics of 7$<Q^{2}<$1000 GeV2, 1.1 *10${}^{-4}<x<3.2*10^{-2}$ and
0.2$<x_{L}<$1\. In this process, positron and proton energies are $E_{e}$=27.5
GeV and $E_{p}$=820 GeV, respectively and the value of the neutron scattering
angle,$\theta_{n}$ , is smaller than 0.8 mrad. The transverse momentum of the
neutron is given by $p_{T}\sim x_{L}E_{p}\theta_{n}$. In this globl analysis,
only H1-2010 data are included in fitting procedure. Fig. (2) confirms the
validity of the fractal model in the kinematic region of x$<$0.01 and
Q${}^{2}>$7GeV2. Therefore, we have limited ourselves to the 189 data points
of H1 collaboration, which are compatible with this kinematic region.

###  III.3 $\chi^{2}$ minimization and nFFs Uncertainties

The optimum parameter values of the nFFs at the initial scale of
Q${}_{0}^{2}$=1 GeV2, we used the $\chi^{2}$ function:

$\chi^{2}=\sum_{i=1}^{N}\frac{(D_{i}-T_{i})^{2}}{\sigma_{i}^{2}}$ (26)

where i labels the number of data points used in our analysis, $D_{i}$ is the
data points and $T_{i}$ is the theoretical prediction, which depends on the
input nFFs parameters. $\sigma_{i}^{2}$ are the experimental and statistical
errors on data points, added in quadrature. In our analysis, the $\chi^{2}$
function is optimized by the CERN program MINUIT James:1994vla . In table 1,
we give the best fit values of nFFs parameters at initial scale of
Q${}_{0}^{2}$=1 GeV2. Those parameters marked by (∗) are set as fixed values
in the final analysis. In Sec. IV we will discuss these obtained parameters in
more detail. In high-energy processes with initial state hadrons, PDFs and
reliable knowledge of their uncertainties are needed. PDF uncertainties stem
from two sources: a) Theoretical uncertainties which are the consequences of
several assumptions made, such as the truncation of the DGLAP perturbation
expansion Martin:2003sk  b) Experimental uncertainties which are due to the
statistical and systematical errors of the experimental data used in the
global fit Martin:2002aw . Besides various techniques such as Lagrange
Multiplier method Stump:2001gu  and Neural Network, the Hessian method is also
a reliable method to calculate the experimental uncertainties of PDFs
Pumplin:2001ct ; Martin:2009iq . We have used it to extract the nFF
uncertainties. We refer the reader to Refs Pumplin:2001ct ; Martin:2009iq  for
technical detail. In our NLO QCD analysis, the covariance matrix elements for
9 free parameters are given in table 2.

Table 1:  The best values, obtained for free parameters in Eq. (III.1) at the initial scale of Q${}_{0}^{2}$ = 1 GeV2 using the H1 experimental data.  Parameters | $\beta M^{n}_{\frac{i=q}{P}}(\beta,Q_{0}^{2},x_{L})$ | $\beta M^{n}_{\frac{i=g}{P}}(\beta,Q_{0}^{2},x_{L})$
---|---|---
$\Gamma_{i}$ | $0.14\pm 0.02$ | $0.61\pm 0.04$
$A_{i}$ | $0^{*}$ | $0^{*}$
$B$ | $1.75\pm 0.02$ | $1.75\pm 0.02$
$C$ | $26.3\pm 1.1$ | $26.3\pm 1.1$
$D$ | $6.67\pm 0.10$ | $6.67\pm 0.10$
$q_{0}$ | $0.32\pm 0.02$ | $0.32\pm 0.02$
$D_{1}^{i}$ | $0.033\pm 0.009$ | $0^{*}$
$D_{2}^{i}$ | $1.015\pm 0.008$ | $1^{*}$
$D_{3}$ | $-0.62\pm 0.02$ | $-0.62\pm 0.02$
Table 2:  The correlation matrix elements for the 9 free parameters at NLO fit  |  | $\Gamma_{q}$ | $\Gamma_{g}$ | $B$ | $C$ | $D$ | $q_{0}$ | $D_{1}^{q}$ | $D_{2}^{q}$ | $D_{3}$
---|---|---|---|---|---|---|---|---|---|---
$\Gamma_{q}$ |  | 1 |  |  |  |  |  |  |  |
$\Gamma_{g}$ |  | 0.274 | 1 |  |  |  |  |  |  |
$B$ |  | 0.070 | 0.003 | 1 |  |  |  |  |  |
$C$ |  | -0.060 | 0.026 | 0.795 | 1 |  |  |  |  |
$D$ |  | -0.022 | 0.024 | 0.254 | 0.534 | 1 |  |  |  |
$q_{0}$ |  | -0.503 | -0.806 | 0.127 | 0.097 | 0.138 | 1 |  |  |
$D_{1}^{q}$ |  | 0.497 | 0.152 | -0.008 | -0.055 | -0.049 | -0.271 | 1 |  |
$D_{2}^{q}$ |  | -0.206 | 0.063 | -0.019 | 0.014 | 0.003 | -0.041 | 0.145 | 1 |
$D_{3}$ |  | -0.340 | -0.583 | 0.093 | 0.068 | 0.099 | 0.243 | -0.012 | -0.021 | 1

##  IV Results and discussion

Let us now give the details of the fractal analysis. The optimum values
obtained in our analysis are summarized in Table 1.

In our analysis, the data set does not include the small values of $x_{L}$, so
we choose the zero value for Aq and Ag parameters. We also obtained a small
value with a large error for D${}_{1}^{g}$. Therefore, we set its value equal
to zero. Now the significant parameter is D${}_{2}^{g}$ which controls the
small-$\beta$ behavior of gluon FFs. This parameter has so much effect on
$\chi^{2}$, so we set its value equal to 1. In this analysis, we obtained the
value of normalized $\chi^{2}$ or the $\chi^{2}$ per degree of freedom as
$\frac{\chi^{2}}{n.d.f}=\frac{217.91}{180}=1.21$ which is almost reasonable
value. Our inputs in Eq (III.1) consist of 13 parameters. In final global
analysis we decided to fix 3 parameters. In comparison with the previous
analyses Ceccopieri:2014rpa ; Shoeibi:2017lrl , the dependence of input to the
fixed parameters reduced (in Ref. Ceccopieri:2014rpa , input describe by 11
parameters and 4 parameters are considered to be fixed and in Ref.
Shoeibi:2017lrl  there is 16 parameters and 8 parameters are chosen to be
fixed). The singlet and gluon densities, $\beta{\cal
M}^{n}_{i/p}(x,Q_{0}^{2};x_{L})$, resulting from our QCD analysis are shown in
Fig. (3) at the initial scale of Q${}_{0}^{2}$=1 GeV2 and at two different
values of $x_{L}$=0.455 and 0.725. The shaded bands correspond to the
estimated uncertainties of at $\Delta\chi^{2}$=1 and $\Delta\chi^{2}$=10.41
for confidence region of 68% (see appendix A for more detail).

Figure 3: (Color online). The singlet and gluon FFs as a function of $\beta$
at the initial scale of Q${}_{0}^{2}$=1 GeV2 for two different values of
$x_{L}$=0.455, 0.725.

As shown in Fig. (4), our results, based on the fractal model, for $\beta{\cal
M}^{n}_{i/p}(x,Q_{0}^{2};x_{L})$ is compared with those of SKTJ17 and F.
Ceccopieri Models Ceccopieri:2014rpa ; Shoeibi:2017lrl . The fractal nFFs seem
valid up to the vertical line is shown in the figures, which confirm the
validity of this model at low values of $\beta$. It appears that the upward
behavior of the singlet fractal densities at low values of $\beta$ (or x)
consistent with those we expected from the behavior of the usual PDFs at low
x.

Figure 4:  (Color online) Comparing the singlet and gluon momentum
distribution resulting from the fractal model with SKTJ17 model
Shoeibi:2017lrl  and prediction proposed by F.Ceccopieri Ceccopieri:2014rpa

Fig. (5) nicely shows the self-similarity behavior of the nFFs in the fixed
value of $x_{L}$. The left panel represents the log-log plot of the nFFs for
singlet and gluon densities as a function of x for different values of Q2. The
same log-log plot as a function of Q2 for different values of x is shown in
the right panel too.

Figure 5:  (Color online) The self-similarity behavior of nFFs in the fixed
value of $x_{L}$=0.5 for different values of x in some fixed values of Q2
(left panel) and also in different values of Q2 for some fixed values of x
(right panel).

In Fig. (6), we present the comparison of our theoretical predictions for
Leading neutron transverse structure functions,
$F_{2}^{Ln}(\beta,x_{L},Q^{2})$, with H1-2010 experimental data Aaron:2010ab .
In general, there is agreement between our analysis based on the fractal
approach with these experimental data in different values of $\beta$, Q2 and
$x_{L}$.

Figure 6:  (Color online) The comparison of H1-2010 experimental data with
Leading neutron transverse structure function, $F_{2}(\beta,Q^{2};x_{L})$, as
a function of $\beta$ for some selected values of Q2 at fixed values of xL =
0.455, 0.725. To facilitate the graphical presentation, we have plotted
$F_{2}(\beta,Q^{2};x_{L})$*3i.

Comparison between our result with those from experimental data and SKTJ17
model are shown in Fig. (7). The shaded bands of both models correspond to the
estimated uncertainties of $\Delta\chi^{2}$=1. It seems that the uncertainties
in the fractal approach reduce about 50%.

Figure 7:  (Color online) The comparison of H1-2010 experimental data with the
experimental observables resulting from the fractal model and SKTJ17 model
Shoeibi:2017lrl , as a function of $\beta$ for some selected values of Q2=11,
37, 82 GeV2 at fixed value of xL = 0.725.

Finally, as a result of this analysis, we show that the nFFs can be described
with a fractal or self-similar model at the low values of the fractional
momentum $\beta$. Recently, we used the fractal approach to explore the
behavior of un-polarized fragmentation functions (FFs) of the pion in the
regions of small momentum fractions z Mohamaditabar:2020btr . We believe that
it is possible to use the fractal framework to explore the low x behavior of
the PDFs too. Work along this line is in progress. In summary, we show that
the upward behavior of the nFFs for $x<0.001$ and $Q^{2}>7GeV^{2}$ may
described by the self-similar behavior of these conditional PDFs. Such
behavior has not obtained in the last studies.

##  V Summary and Conclusions

Fractals provide a different way to observe and model complex phenomena. Their
applications in strong interaction open detail of high energy collisions such
as ep scattering. In 2010 and the years before, several experiments are
dedicated to the positron-proton collider at HERA in order to collect
precision data for events with the final state related to the Leading neutron
production. In this paper, we present the NLO QCD analysis of neutron FFs
based on the fractal approach. We have performed an analysis using the H1-2010
Leading neutron production data. As a result, we have shown that the
parameterized forms of the neutron FFs based on the fractal concept can give
an appropriate description of the Leading neutron production mechanism. We
used the Hessian method in order to explain the uncertainties of the nFFs and
the corresponding observables.

## Acknowledgements

We gratefully acknowledge to Professor Firooz Arash for carefully reading the
manuscript, many helpful discussions, and comments. F. Taghavi Shahri also is
thankful to Ferdowsi University of Mashhad for the financial support provided
for this research. This work is supported by the Ferdowsi University of
Mashhad under grant number 73209(25/12/1397).

## Appendix A

The experimental uncertainties of fractal nFFs are estimated by a general
statistical method named as Hessian method. If the optimized parameters of
each distribution obtained in this analysis present by $p^{0}_{i}$ (i=1,2,3,
..,N), the quadratic expansion of the $\chi^{2}$ around the minimum point
$p^{0}$ is written as:

$\Delta\chi^{2}=\chi^{2}(p^{0}+\delta
p)-\chi^{2}(p^{0})=\sum_{i,j}H_{ij}\delta p_{i}\delta p_{j}\\\ $

In this Taylor series expansion, the first derivative terms vanish at the
minimum point and $H_{ij}$ is the Hessian matrix defined as

$H_{ij}=\frac{1}{2}\frac{\partial^{2}\chi^{2}}{\partial p_{i}\,\partial
p_{j}}\Bigg{|}_{min}\\\ $

The uncertainty on a quantity $F(\textit{x},p^{0})$ is then calculated by
linear error propagation:

$\Delta F=\left[\Delta\chi^{2}_{\rm global}\,\sum_{i,j=1}^{k}\frac{\partial
F}{\partial p_{i}}\,H^{-1}_{ij}\,\frac{\partial F}{\partial
p_{j}}\right]^{\frac{1}{2}}\,.\\\ $

in order to determine the nFFs uncertainties, the confidence region of
$\chi^{2}$ distribution should be estimated. This region can be identified by
an ellipsoid on this normal distribution and it is defined by
$\Delta\chi^{2}$. The confidence level P for the $\chi^{2}$ distribution with
N degrees of freedom is written as:

$P=\int_{0}^{\Delta\chi^{2}}\frac{(\chi^{2})^{\frac{N}{2}-1}}{2^{\frac{N}{2}}\Gamma(\frac{N}{2})}e^{\frac{-\chi^{2}}{2}}d\chi^{2}=P(\approx
0.68)~{},\\\ $

In this analysis, we consider the one-$\sigma$-error range, therefore the
confidence level value is $P=0.68$. Since there are 9 free parameters in the
final fitting procedure, by using the above equation, we obtain the value of
$\Delta\chi^{2}=10.41$.

## References

  * (1) C. Bourrely, Phys. Rev. C 98, no. 5, 055202 (2018) doi:10.1103/PhysRevC.98.055202 [arXiv:1810.08445 [hep-ph]].
  * (2) R. Abdul Khalek et al., arXiv:1905.04311 [hep-ph].
  * (3) F. Zaidi, H. Haider, M. Sajjad Athar, S. K. Singh and I. Ruiz Simo, arXiv:1903.09000 [hep-ph].
  * (4) L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, “Parton distributions in the LHC era: MMHT 2014 PDFs,” Eur. Phys. J. C 75, no. 5, 204 (2015) doi:10.1140/epjc/s10052-015-3397-6 [arXiv:1412.3989 [hep-ph]].
  * (5) R. D. Ball et al. [NNPDF Collaboration], “Parton distributions for the LHC Run II,” JHEP 1504, 040 (2015) doi:10.1007/JHEP04(2015)040 [arXiv:1410.8849 [hep-ph]].
  * (6) A. D. Martin, W. J. Stirling, R. S. Thorne and G. Watt, “Parton distributions for the LHC,” Eur. Phys. J. C 63, 189 (2009) doi:10.1140/epjc/s10052-009-1072-5 [arXiv:0901.0002 [hep-ph]].
  * (7) S. Alekhin, J. Blümlein and S. Moch, Eur. Phys. J. C 78, no. 6, 477 (2018) doi:10.1140/epjc/s10052-018-5947-1 [arXiv:1803.07537 [hep-ph]].
  * (8) M. Salajegheh, M. Nejad, A. Mirjalili and S. Atashbar Tehrani, arXiv:1904.03439 [hep-ph].
  * (9) I. Abt, A. M. Cooper-Sarkar, B. Foster, V. Myronenko, K. Wichmann and M. Wing, Phys. Rev. D 94, no. 3, 034032 (2016) doi:10.1103/PhysRevD.94.034032 [arXiv:1604.02299 [hep-ph]].
  * (10) H. Abdolmaleki et al. [xFitter Developers’ Team], Eur. Phys. J. C 78, no. 8, 621 (2018) doi:10.1140/epjc/s10052-018-6090-8 [arXiv:1802.00064 [hep-ph]].
  * (11) S. Heidari, B. Rezaei, J. K. Sarma and G. R. Boroun, Nucl. Phys. A 986, 195 (2019). doi:10.1016/j.nuclphysa.2019.03.005
  * (12) I. M. Dremin, Mod. Phys. Lett. A 3, 1333 (1988). doi:10.1142/S0217732388001604
  * (13) P. Lipa and B. Buschbeck, Phys. Lett. B 223, 465 (1989). doi:10.1016/0370-2693(89)91634-1
  * (14) W. Florkowski and R. C. Hwa, Phys. Rev. D 43, 1548 (1991). doi:10.1103/PhysRevD.43.1548
  * (15) J. D. Bjorken, In *Eilat 1994, Proceedings, Deep inelastic scattering and related subjects* 151-176, and SLAC Stanford - SLAC-PUB-6477 (94/04,rec.Jul.) 33 p
  * (16) A. Bhattacharya, B. Chakrabarti, T. Sarkhel and S. N. Banerjee, Int. J. Mod. Phys. A 15, 2053 (2000). doi:10.1016/S0217-751X(00)00083-5
  * (17) B. B. Mandelbrot, The Fractal Geometry of Nature., W. H. Freeman, New York 1977.
  * (18) I. M. Dremin and B. B. Levchenko, Phys. Lett. B 292, 155 (1992). doi:10.1016/0370-2693(92)90623-C
  * (19) M. Ninomiya and K. Watanabe, Nuovo Cim. A 37, 343 (1977) doi:10.1007/BF02909267
  * (20) T. Lastovicka, Eur. Phys. J. C 24, 529 (2002) doi:10.1007/s10052-002-8893-9 [hep-ph/0203260].
  * (21) T. Lastovicka, doi:10.3204/DESY-THESIS-2004-016
  * (22) D. K. Choudhury and R. Gogoi, hep-ph/0310260.
  * (23) D. K. Choudhury and R. Gogoi, hep-ph/0503047.
  * (24) A. Jahan and D. K. Choudhury, arXiv:1106.1145 [hep-ph].
  * (25) D. K. Choudhury and A. Jahan, Int. J. Mod. Phys. A 28, 1350079 (2013) doi:10.1142/S0217751X13500796 [arXiv:1305.6180 [hep-ph]].
  * (26) A. Jahan and D. K. Choudhury, J. Phys. Conf. Ser. 481, 012025 (2014). doi:10.1088/1742-6596/481/1/012025
  * (27) A. Jahan and D. K. Choudhury, Phys. Rev. D 89, no. 1, 014014 (2014) doi:10.1103/PhysRevD.89.014014 [arXiv:1401.4327 [hep-ph]].
  * (28) C. Adloff et al. [H1], Phys. Lett. B 520, 183-190 (2001) doi:10.1016/S0370-2693(01)01074-7 [arXiv:hep-ex/0108035 [hep-ex]].
  * (29) C. Adloff et al. [H1], Eur. Phys. J. C 21, 33-61 (2001) doi:10.1007/s100520100720 [arXiv:hep-ex/0012053 [hep-ex]].
  * (30) D. K. Choudhury, B. Saikia and K. Kalita, Int. J. Mod. Phys. A 32, no. 18, 1750107 (2017) doi:10.1142/S0217751X1750107X [arXiv:1608.02771 [hep-ph]].
  * (31) A. Deppman and E. Megias, EPJ Web Conf. 141, 01011 (2017) doi:10.1051/epjconf/201714101011 [arXiv:1610.09928 [hep-ph]].
  * (32) D. de Florian, M. Epele, R. J. Hernandez-Pinto, R. Sassot and M. Stratmann, Phys. Rev. D 95, no. 9, 094019 (2017) doi:10.1103/PhysRevD.95.094019 [arXiv:1702.06353 [hep-ph]].
  * (33) V. Bertone et al. [NNPDF Collaboration], Eur. Phys. J. C 77, no. 8, 516 (2017) doi:10.1140/epjc/s10052-017-5088-y [arXiv:1706.07049 [hep-ph]].
  * (34) D. P. Anderle, F. Ringer and M. Stratmann, Phys. Rev. D 92, no. 11, 114017 (2015) doi:10.1103/PhysRevD.92.114017 [arXiv:1510.05845 [hep-ph]].
  * (35) J. J. Ethier, N. Sato and W. Melnitchouk, Phys. Rev. Lett. 119, no. 13, 132001 (2017) doi:10.1103/PhysRevLett.119.132001 [arXiv:1705.05889 [hep-ph]].
  * (36) G. Altarelli, R. K. Ellis, G. Martinelli and S. Y. Pi, Nucl. Phys. B 160, 301 (1979). doi:10.1016/0550-3213(79)90062-2
  * (37) D. de Florian and R. Sassot, Phys. Rev. D 56, 426 (1997) doi:10.1103/PhysRevD.56.426 [hep-ph/9703228].
  * (38) A. Daleo and R. Sassot, Nucl. Phys. B 673, 357 (2003) doi:10.1016/j.nuclphysb.2003.09.007 [hep-ph/0309073].
  * (39) L. Trentadue and G. Veneziano, Phys. Lett. B 323, 201 (1994). doi:10.1016/0370-2693(94)90292-5
  * (40) L. Trentadue, Nucl. Phys. Proc. Suppl. 39BC, 50 (1995) doi:10.1016/0920-5632(95)00043-9 [hep-ph/9506324].
  * (41) C. Adloff et al. [H1 Collaboration], Eur. Phys. J. C 6, 587 (1999) doi:10.1007/s100529901072 [hep-ex/9811013].
  * (42) F. D. Aaron et al. [H1 Collaboration], Eur. Phys. J. C 68, 381 (2010) doi:10.1140/epjc/s10052-010-1369-4 [arXiv:1001.0532 [hep-ex]].
  * (43) S. Chekanov et al. [ZEUS Collaboration], JHEP 0906, 074 (2009) doi:10.1088/1126-6708/2009/06/074 [arXiv:0812.2416 [hep-ex]].
  * (44) S. Chekanov et al. [ZEUS Collaboration], Nucl. Phys. B 658, 3 (2003) doi:10.1016/S0550-3213(03)00152-4 [hep-ex/0210029].
  * (45) S. Chekanov et al. [ZEUS Collaboration], Nucl. Phys. B 637, 3 (2002) doi:10.1016/S0550-3213(02)00439-X [hep-ex/0205076].
  * (46) L. Rinaldi, doi:10.3204/DESY-THESIS-2006-028
  * (47) S. Shoeibi, H. Khanpour, F. Taghavi-Shahri and K. Javidan, Phys. Rev. D 95, no. 7, 074011 (2017) doi:10.1103/PhysRevD.95.074011 [arXiv:1703.04369 [hep-ph]].
  * (48) C. Adloff et al. [H1], Eur. Phys. J. C 6, 587-602 (1999) doi:10.1007/s100529901072 [arXiv:hep-ex/9811013 [hep-ex]].
  * (49) D. de Florian and R. Sassot, Phys. Rev. D 58, 054003 (1998) doi:10.1103/PhysRevD.58.054003 [hep-ph/9804240].
  * (50) F. A. Ceccopieri and L. Trentadue, Phys. Lett. B 655, 15 (2007) doi:10.1016/j.physletb.2007.07.074 [arXiv:0705.2326 [hep-ph]].
  * (51) G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977). doi:10.1016/0550-3213(77)90384-4
  * (52) J. A. M. Vermaseren, A. Vogt and S. Moch, Nucl. Phys. B 724, 3 (2005) doi:10.1016/j.nuclphysb.2005.06.020 [hep-ph/0504242].
  * (53) J. Benecke, T. T. Chou, C. N. Yang and E. Yen, Phys. Rev. 188, 2159 (1969). doi:10.1103/PhysRev.188.2159
  * (54) T. T. Chou and C. N. Yang, Phys. Rev. D 50, 590 (1994). doi:10.1103/PhysRevD.50.590
  * (55) F. Arash, Phys. Lett. B 557, 38 (2003) doi:10.1016/S0370-2693(03)00177-1 [hep-ph/0301260].
  * (56) J. R. McKenney, N. Sato, W. Melnitchouk and C. R. Ji, Phys. Rev. D 93, no. 5, 054011 (2016) doi:10.1103/PhysRevD.93.054011 [arXiv:1512.04459 [hep-ph]].
  * (57) P. C. Barry, N. Sato, W. Melnitchouk and C. R. Ji, arXiv:1804.01965 [hep-ph].
  * (58) J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai and W. K. Tung, Phys. Rev. D 65, 014013 (2001) doi:10.1103/PhysRevD.65.014013 [hep-ph/0101032].
  * (59) S. Shoeibi, F. Taghavi-Shahri, H. Khanpour and K. Javidan, Phys. Rev. D 97, no. 7, 074013 (2018) doi:10.1103/PhysRevD.97.074013 [arXiv:1710.06329 [hep-ph]].
  * (60) F. A. Ceccopieri, Eur. Phys. J. C 74, no. 8, 3029 (2014) doi:10.1140/epjc/s10052-014-3029-6 [arXiv:1406.0754 [hep-ph]].
  * (61) K. K. Sharma, A. C. Katoch and R. C. Verma, Z. Phys. C 75, 253 (1997) doi:10.1007/s002880050468 [hep-ph/9507358].
  * (62) M. Botje, Comput. Phys. Commun. 182, 490 (2011) doi:10.1016/j.cpc.2010.10.020 [arXiv:1005.1481 [hep-ph]].
  * (63) D. d’Enterria et al., arXiv:1512.05194 [hep-ph].
  * (64) H. Holtmann, G. Levman, N. N. Nikolaev, A. Szczurek and J. Speth, Phys. Lett. B 338, 363 (1994). doi:10.1016/0370-2693(94)91392-7
  * (65) B. Kopeliovich, B. Povh and I. Potashnikova, Z. Phys. C 73, 125 (1996) doi:10.1007/s002880050301 [hep-ph/9601291].
  * (66) F. James, CERN-D-506, CERN-D506.
  * (67) A. D. Martin, R. G. Roberts, W. J. Stirling and R. S. Thorne, Eur. Phys. J. C 35, 325 (2004) doi:10.1140/epjc/s2004-01825-2 [hep-ph/0308087].
  * (68) A. D. Martin, R. G. Roberts, W. J. Stirling and R. S. Thorne, Eur. Phys. J. C 28, 455 (2003) doi:10.1140/epjc/s2003-01196-2 [hep-ph/0211080].
  * (69) D. Stump, J. Pumplin, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai and W. K. Tung, Phys. Rev. D 65, 014012 (2001) doi:10.1103/PhysRevD.65.014012 [hep-ph/0101051].
  * (70) J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai and W. K. Tung, Phys. Rev. D 65, 014013 (2001) doi:10.1103/PhysRevD.65.014013 [hep-ph/0101032].
  * (71) A. Mohamaditabar, F. Taghavi-Shahri and S. Shoeibi, Eur. Phys. J. A 56, no.3, 77 (2020) doi:10.1140/epja/s10050-020-00087-y


11institutetext:  INAF – Osservatorio Astronomico di Roma, Via Frascati 33,
I-00040 Monte Porzio Catone (RM), Italy 22institutetext: ASI – Space Science
Data Center, Via del Politecnico snc, I-00133 Roma, Italy 33institutetext:
INFN – Sezione di Roma Tor Vergata, Via della Ricerca Scientifica 1, I-00133
Roma, Italy

# The quest for dark matter in dwarf spheroidal galaxies with the Cherenkov
Telescope Array

Francesco G. Saturni 1122<EMAIL_ADDRESS>Gonzalo Rodríguez-
Fernández 33 Aldo Morselli for the CTA Consortium 33

###### Abstract

Dwarf spheroidal galaxies are among the best environments that can be studied
with Cherenkov telescopes for indirect searches of $\gamma$-ray signals coming
from dark matter self-interaction (annihilation or decay), due to their
proximity and negligible background emission. We present new determinations of
the dark-matter amount – i.e. the astrophysical factors $J$ and $D$ – in
dwarf-galaxy halos obtained through the MCMC Jeans analysis of their
brightness and kinematic data. Such factors are of great importance to test
the performances of the next-generation $\gamma$-ray instruments such as the
Cherenkov Telescope Array in detecting dark-matter signals from astronomical
environments, or constraining the limits to dark-matter physics parameters
(particle mass and lifetime, annihilation cross section).

##  1 Introduction

The problem of establishing the nature of Dark Matter (DM; Zwi33 ) is one of
the major open challenges in modern astrophysics. Several efforts, both on the
side of elementary particles (e.g. Weakly Interacting Massive Particles,
WIMPs; Apr18 ) and macroscopic objects (Massive Compact Halo Objects, MACHOs;
Tis07 ), have been made to identify plausible DM candidates. However, the
parameter space covered by such candidates ranges over several orders of
magnitude in masses and cross sections Ros04 . The current framework for the
astronomical searches for DM signals is based on the possibility that DM
particles self-interact via annihilation or decay to produce final-state
$\gamma$-ray photons Cir11  whose detection is potentially at reach of next-
generation Cherenkov telescopes such as the Cherenkov Telescope Array (CTA;
Cta17 ). Here we present new estimates of the DM content in the dwarf
spheroidal galaxies (dSphs; Str08 ), which appear to be the best observational
targets for indirect DM searches in terms of proximity and lack of background
emission Eva04 ; Dor13 . In particular, we focus our analysis on the classical
dSph Draco I (Dra I) and the ultra-faint Coma Berenices (CBe), Canes Venatici
I (CVn I), Segue 1 (Seg 1), Reticulum II (Ret II) and Triangulum II (Tri II).
With respect to previous works, all of the results shown here are obtained
with a common methodology for all targets, and will be included, along with
other promising dSphs, in a future comprehensive publication on this research
topic.

Figure 1: Best-fit brightness profiles $\Sigma(R)$ of the analyzed dSphs as a
function of the object’s projected (2D) radial coordinate $R$ from the dSph
centroid. In each panel, the projected 3D profile resulting from the fit (red
line) is shown superimposed to the corresponding data set (black dots; Irw95 ;
Mar08 ; Mun10 ; Bec15 ; Lae15 ).

##  2 Indirect DM detection in the $\gamma$-ray domain

The $\gamma$-ray flux expected from DM self-interaction depends on the type of
DM reaction, the mass of the DM particle and the amount of DM along the line
of sight. In case of annihilation of two DM particles, the latter quantity is
summarized into the astrophysical factors for DM annihilation
$J(\Delta\Omega)$ or decay $D(\Delta\Omega)$:

$\displaystyle J(\Delta\Omega)=\int_{\Delta\Omega}d\Omega\int_{\rm
l.o.s.}\rho_{\rm DM}^{2}(\ell,\Omega)d\ell$ (1) $\displaystyle
D(\Delta\Omega)=\int_{\Delta\Omega}d\Omega\int_{\rm l.o.s.}\rho_{\rm
DM}(\ell,\Omega)d\ell$ (2)

Several attempts to get such values with increasing accuracy exist in the
literature (Acc10, ; Bon15, ; Ger15, ; Mar15, ; Hay16, ). Here we present new
determinations of the astrophysical factors for six of the best dSphs known to
date in the field of DM indirect searches, based on the procedure described in
Bon15  which exploits the capabilities of the Markov-Chain Monte Carlo (MCMC)
Jeans analysis performed by the CLUMPY code Hut18 . At variance with them, we
adopt a common framework for the treatment of surface-brightness and kinematic
data sets of the analyzed dSphs, including the most promising newly-discovered
targets.

Table 1: Top: best-fit parameters – heliocentric distance $d_{\odot}$ McC12 ;
3D brightness scale density $\rho^{*}_{\rm s}$, radius $r^{*}_{\rm s}$ and
profile indexes; dSph virial radius $R_{\rm vir}$; 3D DM density scale
$\rho_{\rm s}$, radius $r_{\rm s}$ and index $\alpha$ Bon15  – for the six
analyzed dSphs, classified for type (‘cls’: classical, ‘uft’: ultra-faint).
Bottom: logarithmic astrophysical factors $\log{J}$ and $\log{D}$ for DM self-
interaction, computed at their optimal angles ($\alpha_{J}=2r_{\rm
s}^{*}/d_{\odot}$ or $\alpha_{D}=r_{\rm s}^{*}/d_{\odot}$ respectively) and at
$0.5$ deg for a direct comparison with the literature Bon15 ; Bon15c ; Gen16 .

Name | Type | $d_{\odot}$ | $\rho_{\rm s}^{*}$ | $r_{\rm s}^{*}$ | $\alpha^{*}$ | $\beta^{*}$ | $\gamma^{*}$ | $R_{\rm vir}$ | $\rho_{\rm s}$ | $r_{\rm s}$ | $\alpha$
---|---|---|---|---|---|---|---|---|---|---|---
|  | (kpc) | ($10^{5}$ L⊙ kpc-3) | (kpc) |  |  |  | (kpc) | ($10^{6}$ M⊙ kpc-3) | (kpc) |
CBe | uft | $44\pm 4$ | $3.7\pm 1.7$ | $0.092\pm 0.008$ | 3.8 | 3.8 | 1.4 | $4.2^{+14.5}_{-2.7}$ | $37^{+54}_{-28}$ | $1.0^{+4.6}_{-0.7}$ | $0.64^{+0.25}_{-0.29}$
CVn I | uft | $218\pm 10$ | $2.15\pm 0.50$ | $0.568\pm 0.033$ | 3.2 | 3.8 | 0.2 | $12.2^{+3.2}_{-4.0}$ | $2.5^{+6.7}_{-1.4}$ | $1.9^{+2.8}_{-1.1}$ | $0.32^{+0.32}_{-0.45}$
Dra I | cls | $76\pm 6$ | $470\pm 140$ | $0.170\pm 0.016$ | 2.8 | 3.8 | 0.2 | $2.68^{+0.45}_{-0.32}$ | $85^{+68}_{-29}$ | $0.28^{+0.33}_{-0.20}$ | $0.41^{+0.23}_{-0.21}$
Ret II | uft | $30\pm 3$ | $157\pm 24$ | $0.023\pm 0.003$ | 3.2 | 3.8 | 0.4 | $1.7^{+5.9}_{-1.1}$ | $65^{+220}_{-50}$ | $0.41^{+2.36}_{-0.34}$ | $0.56^{+0.29}_{-0.27}$
Seg 1 | uft | $23\pm 2$ | $13\pm 11$ | $0.034\pm 0.007$ | 3.6 | 3.8 | 0.8 | $0.28^{+1.94}_{-0.23}$ | $1.6^{+22.7}_{-1.3}$ | $0.32^{+4.02}_{-0.25}$ | $0.53^{+0.33}_{-0.26}$
Tri II | uft | $30\pm 2$ | $45\pm 24$ | $0.020\pm 0.002$ | 3.6 | 3.2 | 0.2 | $6.2^{+15.0}_{-4.7}$ | $800^{+2070}_{-660}$ | $0.64^{+1.76}_{-0.55}$ | $0.52^{+0.32}_{-0.27}$
| $\alpha_{J}$ | $\log{J(\alpha_{J})}$ | $\log{J(0.5\mbox{ deg})}$ | $\alpha_{D}$ | $\log{D(\alpha_{D})}$ | $\log{D(0.5\mbox{ deg})}$
| (deg) | (GeV2 cm-5) | (deg) | (GeV cm-2)
CBe | 0.24 | $19.1^{+0.7}_{-0.6}$ | $19.4^{+0.9}_{-0.7}$ | 0.12 | $17.7^{+0.6}_{-0.4}$ | $18.7^{+0.7}_{-0.6}$
CVn I | 0.30 | $17.1^{+0.2}_{-0.2}$ | $17.2^{+0.2}_{-0.2}$ | 0.15 | $17.0^{+0.1}_{-0.1}$ | $17.6^{+0.2}_{-0.2}$
Dra I | 0.26 | $18.6^{+0.1}_{-0.2}$ | $18.6^{+0.1}_{-0.1}$ | 0.13 | $17.4^{+0.1}_{-0.1}$ | $18.0^{+0.1}_{-0.1}$
Ret II | 0.09 | $18.6^{+0.5}_{-0.5}$ | $19.3^{+0.9}_{-0.7}$ | 0.05 | $16.9^{+0.4}_{-0.3}$ | $18.5^{+0.7}_{-1.0}$
Seg 1 | 0.17 | $16.3^{+2.1}_{-2.1}$ | $16.7^{+2.3}_{-2.5}$ | 0.09 | $16.0^{+1.2}_{-1.2}$ | $17.0^{+1.4}_{-2.2}$
Tri II | 0.08 | $20.8^{+0.6}_{-0.5}$ | $21.6^{+1.0}_{-0.9}$ | 0.04 | $17.8^{+0.5}_{-0.4}$ | $19.6^{+0.8}_{-1.0}$

###  2.1 Input priors and data sets for the CLUMPY set-up

The CLUMPY software allows to perform the dynamical analysis of the DM halos
around dSphs assuming that such galaxies can be considered steady-state
collisionless systems in spherical symmetry and with negligible rotation, in
which the contribution of the stellar component to the total mass can be
neglected. In these conditions, their second-order Jeans equation Bin08  can
be solved to obtain the DM density profile $\rho_{\rm DM}(r)$ along the
object’s 3D radial coordinate $r$, once parametric forms of the stellar number
density, radial velocity dispersion and velocity anisotropy are given. We
follow the prescription by Bon15b , fitting the brightness density $n^{*}(r)$
of each galaxy with a 3D Navarro-Frenk-White (NFW) profile Nav97  projected
onto the corresponding surface brightness (see Fig. 1 and Tab. 1) and adopting
an Einasto profile Ein65  for $\rho_{\rm DM}$, so as to be directly comparable
with the majority of the values of $J$ and $D$ available in the literature.
For the dSph stellar kinematics, we use the same data sets of Bon15 , Bon15c
and Gen16 , also adopting the definition of stellar membership probabilities
$P$ by Bon15 .

###  2.2 Astrophysical factors from the MCMC spherical Jeans analysis of the
dSph halos

We run 10 MC chains with $10^{4}$ realizations each for every target, choosing
the unbinned likelihood analysis Bon15b . We then derive the distribution of
virial radii $R_{\rm vir}$ for each galaxy from the output profiles as made by
Bon15 . Finally, we compute the profiles of $J(\alpha_{\rm int})$ and
$D(\alpha_{\rm int})$ as functions of the integration angle $\alpha_{\rm int}$
by running the CLUMPY executable over the posterior distributions of the DM
profile parameters. In Tab. 1, we report the values of the astrophysical
factors obtained in this way for different integration angles, in order to
favor a direct comparison with the results already available in the
literature. In Fig. 2, we also show the profiles of $J$ and $D$ for the dSph
with the highest expected DM amount (Tri II) and that with the most robust DM
halo parameters (Dra I).

Figure 2: Left panels: DM-annihilation astrophysical factor profiles
$J(\alpha_{\rm int})$ as functions of the integration angle $\alpha_{\rm int}$
for the most massive halo (Tri II) and the halo with the most robust parameter
determination (Dra I). Right panels: equivalent profiles $D(\alpha_{\rm int})$
for DM decay.

##  3 Discussion

Our MCMC examination of the first six DM halos in a completely uniform
framework of input data allows us to produce $J$ and $D$ values for dSphs in
agreement with the literature. The most discrepant result is $J(0.5\mbox{
deg})$ for Dra I, which is offset by a factor of $\sim$4 from e.g. the value
given by Bon15 . Nevertheless, CLUMPY is only able to perform a spherical
Jeans analysis of the DM halos, thus not taking into account additional
uncertainties arising from a possible halo triaxiality ($\sim$0.4 dex). We
also note that, within the analyzed data sets, the dSph Tri II is potentially
observable with CTA as a single targets in $\sim$100 h integration. However,
the limited statistics available on its stellar kinematics (which is also
possibly contaminated by binary stars Kir17 ) prevents to draw any firm
conclusion, again emphasizing the importance of collecting good and clean
stellar data for such targets in order to evaluate them as viable candidates
for DM searches with Cherenkov telescopes.

## References

  * (1) F. Zwicky, Ac. H. Ph. 6, 110 (1933)
  * (2) E. Aprile, J. Aalbers, F. Agostini, et al., Ph. Rv. L. 121, 111302 (2018)
  * (3) P. Tisserand, L. Le Guillou, C. Afonso, et al., A&A 469, 387 (2007)
  * (4) L. Roszkowski, Pramana 62, 389 (2004)
  * (5) M. Cirelli, G. Corcella, A. Hektor, et al., JCAP 3, 051 (2011)
  * (6) CTA Consortium, _Science with the Cherenkov Telescope Array_ (World Scientific Pub., 2019)
  * (7) L.E. Strigari, J.S. Bullock, M. Kaplinghat, et al., Nature 454, 1096 (2008)
  * (8) N.W. Evans, F. Ferrer, S. Sarkar, Ph. Rv. D 69, 123501 (2004)
  * (9) M. Doro, J. Conrad, D. Emmanoulopoulos, et al., A. Ph. 43, 189 (2013)
  * (10) M. Irwin, D. Hatzidimitriou, MNRAS 277, 1354 (1995)
  * (11) N.F. Martin, J.T.A. de Jong, H.W. Rix, ApJ 684, 1075 (2008)
  * (12) R.R. Muñoz, M. Geha, B. Willman, AJ 140, 138 (2010)
  * (13) K. Bechtol, A. Drlica-Wagner, E. Balbinot, et al., ApJ 807, 50 (2015)
  * (14) B.P.M. Laevens, N.F. Martin, R.A. Ibata, et al., ApJL 802, L18 (2015)
  * (15) V.A. Acciari, T. Arlen, T. Aune, et al., ApJ 720, 1174 (2010)
  * (16) V. Bonnivard, C. Combet, M. Daniel, et al., MNRAS 453, 849 (2015)
  * (17) A. Geringer-Sameth, S.M. Koushiappas, M. Walker, ApJ 801, 74 (2015)
  * (18) G.D. Martinez, MNRAS 451, 2524 (2015)
  * (19) K. Hayashi, K. Ichikawa, S. Matsumoto, et al., MNRAS 461, 2914 (2016)
  * (20) M. Hütten, C. Combet, D. Maurin, Co. Ph. C. 235, 336 (2019)
  * (21) A.W. McConnachie, AJ 144, 4 (2012)
  * (22) V. Bonnivard, C. Combet, D. Maurin, et al., ApJL 808, L36 (2015)
  * (23) A. Genina, M. Fairbairn, MNRAS 463, 3630 (2016)
  * (24) J. Binney, S. Tremaine, _Galactic Dynamics: Second Edition_ (Princeton University Press, 2008)
  * (25) V. Bonnivard, C. Combet, D. Maurin, et al., MNRAS 446, 3002 (2015)
  * (26) J.F. Navarro, C.S. Frenk, S.D.M. White, ApJ 490, 493 (1997)
  * (27) J. Einasto, Tr. Alm. 5, 87 (1965)
  * (28) E.N. Kirby, J.G. Cohen, J.D. Simon, et al., ApJ 838, 83 (2017)


11institutetext: Michigan Robotics, University of Michigan, Ann Arbor, MI
48109, USA 11email<EMAIL_ADDRESS>
22institutetext: United Imaging Research Institute of Intelligent Imaging,
Beijing 100144, China
33institutetext: School of Electrical and Data Engineering, University of
Technology, NSW 2007, Australia
33email<EMAIL_ADDRESS>44institutetext: Hamlyn Centre for Robotic
Surgery, Imperial College London, London SW7 2AZ, UK
44email<EMAIL_ADDRESS>
# Bayesian dense inverse searching algorithm for real-time stereo matching in
minimally invasive surgery

Jingwei Song 1122 Qiuchen Zhu 33 Jianyu Lin 44 Maani Ghaffari 11

###### Abstract

This paper reports a CPU-level real-time stereo matching method for surgical
images ($10$ Hz on $640\times 480$ image with a single core of i5-9400). The
proposed method is built on the fast LK algorithm, which estimates the
disparity of the stereo images patch-wisely and in a coarse-to-fine manner. We
propose a Bayesian framework to evaluate the probability of the optimized
patch disparity at different scales. Moreover, we introduce a spatial Gaussian
mixed probability distribution to address the pixel-wise probability within
the patch. In-vivo and synthetic experiments show that our method can handle
ambiguities resulted from the textureless surfaces and the photometric
inconsistency caused by the non-Lambertian reflectance. Our Bayesian method
correctly balances the probability of the patch for stereo images at different
scales. Experiments indicate that the estimated depth has similar accuracy and
fewer outliers than the baseline methods in the surgical scenario with real-
time performance. The code and data set are available at
https://github.com/JingweiSong/BDIS.git.

###### Keywords:

Stereo matching Bayesian theory Posterior probability inference.

##  1 Introduction

Real-time 3D intra-operative tissue surface shape recovery from stereo images
is important in Computer Assisted Surgery (CAS). The reconstructed depth is a
crucial for dense Simultaneous Localization and Mapping (SLAM) [23, 24], AR
system [11, 28] and diseases diagnosis [18, 13]. All stereo matching
procedures follow the pinhole camera model [2] and conduct image
rectification, undistortion, and disparity estimation. The stereo matching
techniques are normally classified into two categories regarding disparity
estimation: prior-free and learning-based. Conventional prior-free methods
estimate the pixel-wise disparity using the image alignment techniques [12,
25, 8, 14, 20]. Based on the left-right image consistency assumption (photo-
metric or feature-metric), they either use corner feature registration, dense
direct pixel searching, or a combination. Differently, Deep Neural Network
(DNN) based techniques directly learn the disparity from the training image
pairs [30, 26, 29, 3, 16]. Although DNN methods are reported to be efficient,
the results may be invalidated with changing parameters such as focal length
and baseline or a large texture difference between the training and testing
data [19, 1]. Moreover, the DNN-based methods heavily depend on the size and
quality of the annotated training data, which are not accessible in many CAS
scenarios.

In the category of prior-free methods, ELAS [8] is still one of the most
widely used stereo matching algorithms due to its robustness and accuracy [23,
24, 33, 32]. It is also the most popular method in the industry [31, 4]. ELAS
uses Sobel descriptors to match sparse corners as the supporting points and
triangulate the pixel-wise disparity prior. Then, the optimal dense disparity
is retrieved with its proposed maximum a-posteriori algorithm. Its two-step
process requires around $0.25-1$ second on a single modern CPU core. This
paper aims for a faster CPU-based stereo matching method.

The Dense Inverse Searching (DIS) [14] shows the potential of dense direct
matching without the time-consuming sparse supporting points alignment. By
resizing the left and right images to several coarse scales, it adopts and
modifies the Lucas-Kanade (LK) optical flow algorithm [17] for fast estimating
the pixel-wise optimal disparity. [14] demonstrates that real-time computation
is possible with its patch-based coarse-to-fine dense matching, where patch
refers to an arbitrary squared image segment. However, DIS is strictly built
based on the photometric consistency and surface texture abundance
assumptions, which cannot always be satisfied in CAS. The two main challenges
are the textureless/dark surfaces and the serious non-Lambertian reflectance.
The weak/dark texture, which widely exists in CAS, leads to ambiguous
photometric consistency. Meanwhile, non-Lambertian reflectance brings uneven
disturbance on the surfaces, and it cannot be eliminated by just enforcing the
patch normalization [21].

In this paper, to deal with photometric inconsistency and non-Lambertian
reflectance in stereo matching, we propose a Bayesian Dense Inverse Searching
(BDIS) to quantify the posterior probability of each optimized patch. A
spatial Gaussian Mixture Model (GMM) is further adapted to quantify pixel-wise
confidence within the patch. The final pixel-wise disparity is the fusion of
multiple local overlapping patches, reducing the impact of those patches
suffering from the textureless/dark surfaces or the non-Lambertian
reflectance. In extreme cases, it is beneficial to give up the disparity
estimation of some patches identified as dubious. In particular, this work has
the following contributions:

  * •

A Bayesian approach is developed to quantify the posterior probability of the
patch.

  * •

A spatial GMM is introduced to quantify the pixels’ confidence within the
patch.

  * •

To our knowledge, BDIS is the first single core CPU based stereo matching
approach that achieves similar performance to the near real-time method ELAS.

##  2 Methodology

###  2.1 Multiscale DIS

Fig. 1 shows the DIS (based on fast LK) algorithm for stereo matching proposed
by [14]. It is a modified version of the LK algorithm. We use the fast DIS as
our base framework. Note that the variational refinement module in [14] is
abandoned because it has a small (less than $0.5\%$) contribution in promoting
the accuracy. The modified fast LK based DIS is achieved by minimizing the
following objective function:

$\Delta\mathbf{u}=\operatorname{argmin}_{\Delta\mathbf{u}^{\prime}}\sum_{x}\left[I_{r}\left(\mathbf{x}+\mathbf{u}\right)-I_{l}(\mathbf{x}+\Delta\mathbf{u}^{\prime})\right]^{2},$
(1)

where $\mathbf{x}$ is the processed location, $\mathbf{u}$ is the estimated
disparity in the loop, $I_{l}$ and $I_{r}$ are the left image patch and right
image, and $\Delta\mathbf{u}$ is the optimal update of $\mathbf{u}$ at one
loop. Different from authentic LK, $\Delta\mathbf{u}^{\prime}$ is moved from
the right image to the left image patch. The improvement avoids the expensive
re-evaluation of the Hessian on the right image. (1) is traversed on all
patches at different scales. The disparity at the fine-scale level is
initialized at the optimized coarse scale. The optimal disparity at the
location $\mathbf{x}$ is the weighted fusion with all covering patches using
inverse residual:

$\hat{\mathbf{u}}_{\mathbf{x}}=\sum_{k\in\Omega}\frac{1/\operatorname{max}(\lVert
I_{l}(\mathbf{x}+\mathbf{u}^{(k)})-I_{r}(\mathbf{x})\lVert^{2},1)}{\sum_{k\in\Omega}1/\operatorname{max}(\lVert
I_{l}(\mathbf{x}+\mathbf{u}^{(k)})-I_{r}(\mathbf{x})\lVert^{2},1)}\mathbf{u}^{(k)},$
(2)

where $\Omega$ is the set of patches covering the position $\mathbf{x}$,
$\mathbf{u}^{(k)}$ is the estimated disparity of the patch $k$ and
$\operatorname{max}(\cdot,\cdot)$ selectes the maximum value. The pixel-wise
disparity $\hat{\mathbf{u}}_{\mathbf{x}}$ is the weighted average of the
estimated disparities from all patches, wherein the weight is the inverse
residual of brightness.

(a)  Figure 1: The framework of the DIS algorithm [14]. It uses 3 scale levels
as an example.

###  2.2 The Bayesian patch-wise posterior probability

The residual-based weighted average fusion (2) suffers from the ambiguities
brought by the textureless/dark surface and non-Lambertian reflectance. The
textureless/dark surface leads to ambiguous local minima of the cost function
penalizing photometric inconsistency (1) and misleads the algorithm to be
over-confident on the estimation. Furthermore, the photometric consistency
presumption is seriously violated on the surface affected heavily by the non-
Lambertian reflectance. The affine lighting changes formulation in previous
large-scale SLAM studies [7] cannot fully tackle the complex and severe non-
Lambertian reflectance in CAS. In both situations, the weights retrieved from
the photometric residuals (2) are misleading. To overcome the difficulty in
defining the confidence of the estimated disparity, we propose a Bayesian
model to correctly estimate the confidence in the presence of textureless
surface and non-Lambertian reflectance. Since the uncertainty distribution of
both the left and right scenes is unclear, it is difficult to conduct the
direct inference of the posterior probability in terms of disparity. Thus, we
implicitly infer the probability with Bayesian modeling using Conditional
Random Fields (CRF) [27]. The posterior probability of the patch-wise
disparity $\mathbf{u}^{(k)}$ is

$\displaystyle
p(\mathbf{u}^{(k)}|I_{l},I_{r})\propto\frac{p(I_{r}|I_{l},\mathbf{u}^{(k)})}{p(I_{r},I_{l},\mathbf{u}^{(k)})}\propto\frac{p(I_{r}|I_{l},\mathbf{u}^{(k)})}{\Sigma_{\mathbf{u}_{i}^{(k)}\in\mathcal{P}}p(I_{r}|I_{l},\mathbf{u}_{i}^{(k)})}\propto\frac{p(I_{r}|I_{l},\mathbf{u}^{(k)})}{\Sigma_{\mathbf{u}_{i}^{(k)}\in\mathcal{P}^{\prime}}p(I_{r}|I_{l},\mathbf{u}_{i}^{(k)})}\mathrm{r},$
(3)

where $\mathcal{P}$ is the domain of all possible choice of
$\mathbf{u}_{i}^{(k)}$. To reduce computational load, $\mathrm{r}$ is applied
as the constant compensation ratio for all patches within the window.
$\mathcal{P}$ is reduced to a small window $\mathcal{P}^{\prime}$ assuming the
rest candidates are numerically trivial.

Equation (3) indicates that the posterior probability of the disparity can be
obtained by traversing the probability on all possible $\mathbf{u}_{i}^{(k)}$.
And the possible choice of disparity is equal to window size $\mathrm{s}$.
Even though the posterior probability suffers from the textureless surface and
non-Lambertian reflectance, the illumination consistency probability is
proportional to the residuals because the set of neighboring disparities is
within one patch, and the impact of the issues is consistent. Thus, we model
the illumination consistency probability
$p(I_{r}|I_{l},\mathbf{u}_{i}^{(k)},\mathbf{x})$ based on the Boltzmann
distribution [15] as

$p(I_{r}|I_{l},\mathbf{u}_{i}^{(k)})=\exp\left(-\frac{\lVert
I_{l}(\mathbf{u}_{i}^{(k)})-I_{r}(\mathbf{u}_{i}^{(k)})\lVert^{2}_{\mathrm{F}}}{2\sigma_{r}^{2}\mathrm{s}^{2}}\right),$
(4)

where $\lVert\cdot\lVert_{\mathrm{F}}$ is the Frobenius norm and $\sigma_{r}$
is the hyperparameter to describe the variance of the brightness. The relative
posterior probability can be obtained with (3) and (4). Generally, the
absolute exponential parameter of the Boltzmann distribution denotes the
entropy of the state. In our case, such entropy is defined as (4). Image with
abundant texture has more entropy loss. Hence, the entropy item is highly
related to the photometric inconsistency loss.

Fig. 2 shows the relationship between the illumination consistency probability
density function and the texture. The response is stronger on the textured
surface. The residuals are always small in the textureless surface, no matter
how the left and right images are aligned. (2) cannot correctly measure the
weights while (4) describes the relative probability of the estimation.
Moreover, it tests the local convergence to filter the Saddle point solutions.

(a)  Figure 2: The probability density function of the textureless region.

###  2.3 The prior spatial Gaussian probability

In addition to the patch-wise posterior probability of the disparity in the
last section, a spatial GMM is adopted to estimate pixel-wise probability
within the patch. Considering that medical images are natural images, a
multivariate Gaussian distribution is adopted to measure the confidence of the
pixel-wise probability using a Gaussian mask. In accordance with the
multivariate Gaussian distribution, the center of the patch has higher
confidence than the edge pixels since those central pixels preserve more
information for inference. Assuming all pixels in the patch are i.i.d, we have

$\displaystyle p(\mathbf{u}^{(k)}|I_{l},I_{r},\mathbf{x})\propto
p(\mathbf{u}^{(k)}|I_{l},I_{r})\exp\left(-\frac{\sum_{\mathbf{\xi}^{(k)}(\mathbf{x})}\lVert\mathbf{x}-\mathbf{\xi}^{(k)}(\mathbf{x})\lVert^{2}_{\mathrm{F}}}{2\sigma_{s}^{2}}\right),$
(5)

where $\mathbf{\xi}^{(k)}(\mathbf{x})$ is the set of all pixel positions
within the patch $k$ in image coordinate. $\sigma_{s}$ is the 2D spatial
variance of the probability. Note that (5) is independent of the patch and can
therefore be pre-computed before the process. Combining (3), (4) and (5), the
final pixel-wise posterior probability distribution can be represented as
follows,

$\displaystyle
p(\mathbf{u}^{(k)}\\!|\\!I_{l},I_{r},\mathbf{x})\\!\propto\\!\exp\left(\\!-\frac{\\!\sum_{\\!\mathbf{\xi}^{(k)}(\mathbf{x})}\\!\lVert\mathbf{x}-\mathbf{\xi}^{(k)}(\mathbf{x})\lVert^{2}_{\mathrm{F}}}{2\sigma_{s}^{2}}\right)\frac{\\!\exp\left(-\frac{\\!\lVert
I_{l}(\mathbf{u}^{(k)})-I_{r}(\mathbf{u}^{(k)})\\!\lVert^{2}_{\mathrm{F}}}{2\sigma_{r}^{2}\mathrm{s}^{2}}\\!\right)}{\\!\sum_{\mathbf{u}_{i}^{(k)}\\!\in\\!\mathcal{P}}\exp\left(-\frac{\lVert
I_{l}(\mathbf{u}_{i}^{(k)})-I_{r}(\mathbf{u}_{i}^{(k)})\lVert^{2}_{\mathrm{F}}}{2\sigma_{r}^{2}\mathrm{s}^{2}}\\!\right)}.$

Finally, it should be emphasized that (4) and (5) are not the cost functions
but probability/weight for each patch or pixel. Costly optimization steps are
avoided.

##  3 Results and discussion

BDIS was compared with DIS [14], SGBM [12] and ELAS [8] on the in-vivo and the
synthetic data sets 111Readers are encouraged to watch the attached video and
test the code.. The computations were implemented on a commercial desktop
(i5-9400) in C++. DNN-based methods PSMNet [5] and GwcNet [10] were also
compared for completeness and the computation was conducted on the GTX 1080ti
in PyTorch. The public in-vivo stereo videos from [9] were adopted which
contains 200 images with size $640\times 480$ and 200 images with size
$288\times 360$. All stereo images were rectified, undistorted, calibrated,
and vertically aligned with the provided intrinsic and extrinsic parameters.
We also provided a synthetic data set generated from an off-the-shelf virtual
phantom of a male’s digestive system. A virtual handheld colonoscope was
placed inside the colon and was manipulated to go through the colon to collect
the depth and stereo images. The 3D game engine Unity3D222https://unity.com/
was used to generate the sequential stereo and depth images with a pin-hole
camera in size $640\times 480$. The synthetic distortion-free data has
accurate intrinsic and extrinsic parameters. Both diffuse lighting ($100$
frames) and non-Lambertian reflectance ($100$ frames) were simulated. $\gamma$
was set to $0.75$ for $640\times 480$ and $0.25$ for $288\times 360$ data to
discard the patch without enough valid pixels. $\sigma_{r}$ and $\sigma_{s}$
were set to $4$; the sampling within one Bayesian window was 5; the
disturbance from the convergence was $0.5$ and $1$ pixel.

(a)
Figure 3: Sample reconstructions in Diffuse lighting and Lambertian reflectance scenarios. Table 1: The results on the synthetic data with diffused light and Lambertian reflectance. |  |  ELAS |  SGBM |  DIS |  BDIS |  GwcNet |  PSMnet
---|---|---|---|---|---|---|---
Diffuse light |  Median error |  0.178 |  0.512 |  0.251 |  0.161 |  0.542 |  0.417
|  Mean error |  0.220 |  1.113 |  0.753 |  0.320 |  0.809 |  0.641
|  Valid pixels ($1000$) |  166.77 |  103.92 |  288.41 |  208.44 |  100.00 |  301.42
Non-Lambertian reflectance |  Median error |  0.198 |  0.710 |  0.376 |  0.163 |  0.271 |  0.731
|  Mean error |  0.235 |  1.400 |  1.051 |  0.379 |  0.662 |  1.027
|  Valid pixels ($1000$) |  81.50 |  74.29 |  295.92 |  204.42 |  106.38 |  301.46

###  3.1 Quantitative comparisons on the synthetic data set

BDIS was compared quantitatively with the baseline methods ELAS, SGBM, DIS,
PSMNet, and GwcNet. The comparison between the prior-based DNN-based method
and BDIS is for completeness only. The default setting of PSMNet and GwcNet
were strictly followed. The pre-trained networks were adopted and finetuned
with the labeled 300 (training) and 50 (validation) synthetic images for
training and validation. Both were trained with Adam optimizer in 300 epochs.

Table 1 and Fig. 3 show the comparisons on the synthetic data set, which are
unaffected by distortion and inaccurate camera parameters. Considering the
median error, BDIS is the best and has $9.55\%$ and $17.68\%$ higher accuracy
than ELAS in diffuse lighting and non-Lambertian reflectance. The results
indicate that BDIS is more advantageous in the scenario of non-Lambertian
reflectance over ELAS, thus more robust in surgical scenarios. Results also
show that BDIS cannot handle the edges well. Fig. 3 and Table 1 reveal the bad
mean error comparison is attributed to the small group of far-out points on
the dark regions/edges. The number of valid prediction suggest BDIS produces
more predictions but suffers from inaccurate dark region predictions.

Readers may notice the bad performance of DNN, which contradicts the
conclusion from [1]. The reason is that the finetuning training process does
not yield satisfying model parameters. The synthetic data set for transfer
learning and the data used to pre-train the DNN are significantly different in
terms of textures. Studies [6, 22] indicate that the performance of the
convolutional DNN is heavily dependent on the image texture, and efforts were
devoted to bridging the domain gap [6, 34]. The bad training process indicates
its strong dependency on the training data set, which can be avoided using
prior-free methods. Further tests will be conducted on labeled in-vivo data
set.

###  3.2 Qualitative comparisons on the in-vivo dataset

We compared ELAS, BDIS, DIS, and SGBM on the in-vivo data sets. Since no
ground truth is provided, DNN-based methods cannot be implemented. We aim to
show that BDIS achieves similar accuracy as ELAS since near real-time ELAS is
widely used in the community. Based on the scope-to-surface distance, the
samples were categorized into five groups. Results show that BDIS achieves an
average $0.4-1.66mm$ (median error) and $0.65-2.32mm$ (mean error) deviation
from ELAS’s results.

The invalid/dark/bright pixels lead to photometric inconsistency in the stereo
matching process. Fig. 4 shows the qualitative comparisons of ELAS, DIS, SGBM,
and BDIS on the relatively well-textured images. Generally, BDIS achieves
similar performance as ELAS but better matches pixels at the image edge with
fewer outliers. DIS and SGBM suffer from the wrong edges. Invalid pixels
inevitably exist on the edges of the rectified image after the image
undistortion. Thus, in the coarse-level patch disparity estimation, patches
with more invalid pixels are more likely to fail in convergence or yield local
minima (abnormal depth) due to insufficient information. The dubious
predictions, however, substantially influence the prediction and the
initialization of the disparity at the finer-scale patch optimization (as in
(2)). BDIS solves the problem by quantifying the posterior probability,
discarding the patch that does not converge, and lowering the patches’
probabilities with invalid pixels. Although the discarded patch does not help
yield disparity, other patches compensate for the loss. If one pixel is not
covered by any patch, we follow ELAS not to optimize the pixel.

Another noticeable problem is the ambiguous local minima in the cost function,
which penalizes the photometric inconsistency. Fig. 4 shows BDIS has fewer
local minima than DIS and SGBM and is similar to ELAS. Fig. 4 (a-b) indicates
that the BDIS addresses the patchs’ probabilities with textured and alleviates
the ambiguous disparity from the textureless surface. Fig. 4 (c-e) show that
the ambiguities caused by the illumination have been greatly reduced. The
quantitative results also provide evidence on its side. It should be
emphasized that this work does not enforce any prior smoothness constraint in
the optimization process.

We additionally tested BDIS and ELAS on the surfaces with serious non-
Lambertian reflectance (Fig. 5). The photometric consistency of this data
deteriorates significantly. Fig. 5 shows that the center of the soft tissue is
exposed to intense lighting while the marginal region is dark. Fig. 5
indicates that ELAS suffers from the ambiguity on the marginal dark regions
while BDIS can ignore or estimate most dark pixels correctly.

###  3.3 Processing rate comparison

We compared the time consumption of ELAS, DIS, and BDIS on a single core of
CPU (i5-9400). BDIS runs $10Hz$ on $640\times 480$ image and $25Hz$ on
$360\times 288$ image while ELAS achieves $4Hz$ and $11Hz$. The two DNN
methods GwcNet and PSMnet run $3Hz$ and $5Hz$ on GTX 1080ti. BDIS consumes
double the time of DIS. The majority of the extra time of BDIS is devoted to
patch-wise window traversing. Since the sampling window size is 5 in the
experiment, 5 more times residual estimations are needed. In general, BDIS
achieves similar/better performance over ELAS but runs 2 times faster.

(a)  Figure 4: Sample recovered shapes of the 5 classes. Circles mark the
regions with large error.

(a)  Figure 5: The qualitative comparisons on the heavy Lambertian reflectance
and dark case.

##  4 Conclusion

We propose BDIS, the first CPU-level real-time stereo matching approach for
CAS. BDIS inherits the fast performance of DIS while being more robust to
textureless/dark surface and severe non-Lambertian reflectance. It achieves
similar or better performance in accuracy as the near real-time method ELAS. A
Bayesian approach and a spatial GMM are developed to describe the relative
confidence of the pixel-wise disparity to achieve the performance. Experiments
indicate that BDIS has fewer outliers than DIS and achieves a lower amount of
outlier predictions than the near real-time ELAS.

## References

  * [1] Allan, M., Mcleod, J., Wang, C.C., Rosenthal, J.C., Fu, K.X., Zeffiro, T., Xia, W., Zhanshi, Z., Luo, H., Zhang, X., et al.: Stereo correspondence and reconstruction of endoscopic data challenge. arXiv preprint arXiv:2101.01133 (2021)
  * [2] Andrew, A.M.: Multiple view geometry in computer vision. Kybernetes (2001)
  * [3] Brandao, P., Psychogyios, D., Mazomenos, E., Stoyanov, D., Janatka, M.: HAPNet: hierarchically aggregated pyramid network for real-time stereo matching. Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization pp. 1–6 (2020)
  * [4] Cartucho, J., Tukra, S., Li, Y., S. Elson, D., Giannarou, S.: VisionBlender: a tool to efficiently generate computer vision datasets for robotic surgery. Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization pp. 1–8 (2020)
  * [5] Chang, J.R., Chen, Y.S.: Pyramid stereo matching network. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 5410–5418 (2018)
  * [6] Chen, X., Wang, Y., Chen, X., Zeng, W.: S2R-DepthNet: Learning a generalizable depth-specific structural representation. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 3034–3043 (2021)
  * [7] Engel, J., Koltun, V., Cremers, D.: Direct sparse odometry. IEEE Transactions on Pattern Analysis and Machine Intelligence 40(3), 611–625 (2017)
  * [8] Geiger, A., Roser, M., Urtasun, R.: Efficient large-scale stereo matching. In: Asian conference on computer vision. pp. 25–38. Springer (2010)
  * [9] Giannarou, S., Visentini-Scarzanella, M., Yang, G.Z.: Probabilistic tracking of affine-invariant anisotropic regions. IEEE Transactions on Pattern Analysis and Machine Intelligence 35(1), 130–143 (2013)
  * [10] Guo, X., Yang, K., Yang, W., Wang, X., Li, H.: Group-wise correlation stereo network. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 3273–3282 (2019)
  * [11] Haouchine, N., Dequidt, J., Peterlik, I., Kerrien, E., Berger, M.O., Cotin, S.: Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery. In: 2013 IEEE international symposium on mixed and augmented reality (ISMAR). pp. 199–208. IEEE (2013)
  * [12] Hirschmuller, H.: Accurate and efficient stereo processing by semi-global matching and mutual information. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. vol. 2, pp. 807–814. IEEE (2005)
  * [13] Jia, X., Mai, X., Cui, Y., Yuan, Y., Xing, X., Seo, H., Xing, L., Meng, M.Q.H.: Automatic polyp recognition in colonoscopy images using deep learning and two-stage pyramidal feature prediction. IEEE Transactions on Automation Science and Engineering 17(3), 1570–1584 (2020)
  * [14] Kroeger, T., Timofte, R., Dai, D., Van Gool, L.: Fast optical flow using dense inverse search. In: Proceedings of the European Conference on Computer Vision. pp. 471–488. Springer (2016)
  * [15] Larochelle, H., Bengio, Y.: Classification using discriminative restricted boltzmann machines. pp. 536–543 (2008)
  * [16] Long, Y., Li, Z., Yee, C.H., Ng, C.F., Taylor, R.H., Unberath, M., Dou, Q.: E-DSSR: Efficient dynamic surgical scene reconstruction with transformer-based stereoscopic depth perception. arXiv preprint arXiv:2107.00229 (2021)
  * [17] Lucas, B.D., Kanade, T., et al.: An iterative image registration technique with an application to stereo vision. Vancouver, British Columbia (1981)
  * [18] Mahmood, F., Yang, Z., Chen, R., Borders, D., Xu, W., Durr, N.J.: Polyp segmentation and classification using predicted depth from monocular endoscopy. In: Medical Imaging 2019: Computer-Aided Diagnosis. vol. 10950, p. 1095011\. International Society for Optics and Photonics (2019)
  * [19] Pratt, P., Bergeles, C., Darzi, A., Yang, G.Z.: Practical intraoperative stereo camera calibration. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 667–675. Springer (2014)
  * [20] Rappel, J.K.: Surgical stereo vision systems and methods for microsurgery (May 3 2016), US Patent 9,330,477
  * [21] Shimasaki, Y., Iwahori, Y., Neog, D.R., Woodham, R.J., Bhuyan, M.: Generating Lambertian image with uniform reflectance for endoscope image. IWAIT2013 pp. 1–6 (2013)
  * [22] Song, J., Patel, M., Girgensohn, A., Kim, C.: Combining deep learning with geometric features for image-based localization in the gastrointestinal tract. Expert Systems with Applications p. 115631 (2021)
  * [23] Song, J., Wang, J., Zhao, L., Huang, S., Dissanayake, G.: Dynamic reconstruction of deformable soft-tissue with stereo scope in minimal invasive surgery. IEEE Robotics and Automation Letters 3(1), 155–162 (2017)
  * [24] Song, J., Wang, J., Zhao, L., Huang, S., Dissanayake, G.: MIS-SLAM: Real-time large-scale dense deformable SLAM system in minimal invasive surgery based on heterogeneous computing. IEEE Robotics and Automation Letters 3(4), 4068–4075 (2018)
  * [25] Stoyanov, D., Scarzanella, M.V., Pratt, P., Yang, G.Z.: Real-time stereo reconstruction in robotically assisted minimally invasive surgery. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 275–282. Springer (2010)
  * [26] Turan, M., Almalioglu, Y., Araujo, H., Konukoglu, E., Sitti, M.: Deep endovo: A recurrent convolutional neural network (RCNN) based visual odometry approach for endoscopic capsule robots. Neurocomputing 275, 1861 – 1870 (2018). https://doi.org/https://doi.org/10.1016/j.neucom.2017.10.014, http://www.sciencedirect.com/science/article/pii/S092523121731665X
  * [27] Uzunbas, M.G., Chen, C., Metaxas, D.: An efficient conditional random field approach for automatic and interactive neuron segmentation. Medical image analysis 27, 31–44 (2016)
  * [28] Widya, A.R., Monno, Y., Imahori, K., Okutomi, M., Suzuki, S., Gotoda, T., Miki, K.: 3D reconstruction of whole stomach from endoscope video using structure-from-motion. In: 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). pp. 3900–3904. IEEE (2019)
  * [29] Yang, G., Manela, J., Happold, M., Ramanan, D.: Hierarchical deep stereo matching on high-resolution images. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 5515–5524 (2019)
  * [30] Ye, M., Johns, E., Handa, A., Zhang, L., Pratt, P., Yang, G.Z.: Self-supervised siamese learning on stereo image pairs for depth estimation in robotic surgery. arXiv preprint arXiv:1705.08260 (2017)
  * [31] Zampokas, G., Tsiolis, K., Peleka, G., Mariolis, I., Malasiotis, S., Tzovaras, D.: Real-time 3D reconstruction in minimally invasive surgery with quasi-dense matching. In: 2018 IEEE International Conference on Imaging Systems and Techniques (IST). pp. 1–6. IEEE (2018)
  * [32] Zhan, J., Cartucho, J., Giannarou, S.: Autonomous tissue scanning under free-form motion for intraoperative tissue characterisation. In: Proceedings of the IEEE International Conference on Robotics and Automation. pp. 11147–11154. IEEE (2020)
  * [33] Zhang, L., Ye, M., Giataganas, P., Hughes, M., Yang, G.Z.: Autonomous scanning for endomicroscopic mosaicing and 3D fusion. In: Proceedings of the IEEE International Conference on Robotics and Automation. pp. 3587–3593. IEEE (2017)
  * [34] Zheng, C., Cham, T.J., Cai, J.: T2net: Synthetic-to-realistic translation for solving single-image depth estimation tasks. In: Proceedings of the European Conference on Computer Vision. pp. 767–783 (2018)