# Characteristics and Estimates of Double Parton Scattering at the Large
Hadron Collider

Edmond L. Berger<EMAIL_ADDRESS>High Energy Physics Division, Argonne
National Laboratory, Argonne, IL 60439  C. B. Jackson<EMAIL_ADDRESS>High
Energy Physics Division, Argonne National Laboratory, Argonne, IL 60439  Gabe
Shaughnessy<EMAIL_ADDRESS>High Energy Physics Division,
Argonne National Laboratory, Argonne, IL 60439  Department of Physics &
Astronomy, Northwestern University, Evanston, IL 60208

###### Abstract

We evaluate the kinematic distributions in phase space of 4-parton final-state
subprocesses produced by double parton scattering, and we contrast these with
the final-state distributions that originate from conventional single parton
scattering. Our goal is to establish the distinct topologies of events that
arise from these two sources and to provide a methodology for experimental
determination of the relative magnitude of the double parton and single parton
contributions at Large Hadron Collider energies. We examine two cases in
detail, the $b~{}\bar{b}~{}\rm{jet~{}jet}$ and the 4 jet final states. After
full parton-level simulations, we identify a few variables that separate the
two contributions remarkably well, and we suggest their use experimentally for
an empirical measurement of the relative cross section. We show that the
double parton contribution falls off significantly more rapidly with the
transverse momentum $p_{T}^{j1}$ of the leading jet, but, up to issues of the
relative normalization, may be dominant at modest values of $p_{T}^{j1}$ .

††preprint: ANL-HEP-PR-09-109, NU-HEP-TH/09-14

##  I Introduction

Double parton scattering (DPS) means that two short-distance subprocesses
occur in a given hadronic interaction, with two initial partons being active
from each of the incident protons in a collision at the Large Hadron Collider
(LHC). The concept is shown for illustrative purposes in Fig. 1, and it may be
contrasted with conventional single parton scattering (SPS) in which one
short-distance subprocess occurs, with one parton active from each initial
hadron. Since the probability of single parton scattering is itself small, it
is often expected that the chances are considerably suppressed for two or more
short-distance interactions in a given collision. However, expectations such
as these bear quantitative re-examination at the LHC where the high overall
center-of-mass energy provides access to very small values of the fractional
momentum $x$ carried by partons, a region in which parton densities grow
rapidly. A large contribution from double parton scattering could result in a
larger than otherwise anticipated rate for multi-jet production and produce
significant backgrounds in searches for signals of new phenomena. The high
energy of the LHC also provides an increased dynamic range of available phase
space for detailed investigations of DPS.

Figure 1: Sketch of a double-parton process in which the active partons are
$i$ and $k$ from one proton and $j$ and $l$ from the second proton. The two
hard scattering subprocess are $A(i~{}j\rightarrow a~{}b)$ and
$B(k~{}l\rightarrow c~{}d)$.

Investigations of double parton scattering have a long history theoretically
Goebel:1979mi ; Paver:1982yp ; Humpert:1983pw ; Mekhfi:1983az ; Humpert:1984ay
; Ametller:1985tp ; Halzen:1986ue ; Mangano:1988sq ; Godbole:1989ti ;
Drees:1996rw ; Eboli:1997sv ; Yuan:1997tr ; Calucci:1997uw ; DelFabbro:1999tf
; Kulesza:1999zh ; Korotkikh:2004bz ; Cattaruzza:2005nu ; Hussein:2006xr ;
Maina:2009sj ; Domdey:2009bg ; d'Enterria:2009hd ; Gaunt:2009re , and there is
evidence for their presence in collider data from the CERN Intersecting
Storage Rings Akesson:1986iv  and Fermilab Tevatron Abe:1997xk ; D0:2009 . A
significantly greater role for double-parton processes may be expected at the
LHC where higher luminosities are anticipated along with the higher collision
energies. Of substantial importance is to know empirically how large the
double parton contribution may be and its dependence on relevant kinematic
variables.

Our aim is to calculate characteristic final states at LHC energies in which
it may be straightforward to discern a double parton signal. We show in this
paper that double parton scattering produces an enhancement of events in
regions of phase space in which the “background” from single parton scattering
is relatively small. If such enhancements are observed experimentally, with
the kinematic dependence we predict, then we will have a direct empirical
means to measure the size of the double parton contribution. In addition to
its role in general LHC phenomenology, this measurement will have an impact on
the development of partonic models of hadrons, since the effective cross
section for double parton scattering measures the size in impact parameter
space of the incident hadron’s partonic hard core.

From the perspective of sensible rates and experimental tagging, a good
process to examine should be the 4 parton final state in which there are $2$
hadronic jets plus a $b$ quark and a $\bar{b}$ antiquark, viz.
$b~{}\bar{b}~{}j_{1}~{}j_{2}$. If the final state arises from double parton
scattering, then it is plausible that one subprocess produces the
$b~{}\bar{b}$ system and another subprocess produces the two jets. There are,
of course, many single parton scattering (2 to 4 parton) subprocesses that can
result in the $b~{}\bar{b}~{}j_{1}~{}j_{2}$ final state, and we look for
kinematic distributions that show notable separations of the two
contributions. As we show, the correlations in the final state are predicted
to be quite different between the double parton and the single parton
subprocesses. For example, the plane in which the $b~{}\bar{b}$ pair resides
is uncorrelated with the $j_{1}~{}j_{2}$ plane in double parton scattering,
but not in the single parton case.

The state-of-the-art of calculations of single parton scattering is well
developed whereas the phenomenology of double parton scattering is as yet much
less advanced. In the remainder of this Introduction, we first describe the
approach we adopt for the calculation of double parton scattering,
specializing to the proton-proton situation of the LHC. Then we outline the
paper and summarize our main results. Our calculations are done at leading-
order in perturbative QCD, adequate for the points we are trying to make.

Making the usual factorization assumption, we express the single-parton hard-
scattering differential cross section for $p~{}p\rightarrow a~{}b~{}X$ as

$\displaystyle d\sigma^{SPS}=\sum_{i,j}\int
f^{i}_{p}(x_{1},\mu)f^{j}_{p}(x_{1}^{\prime},\mu)d\hat{\sigma}_{(ij\rightarrow
ab)}(x_{1},x_{1}^{\prime},\mu)dx_{1}dx_{1}^{\prime}.$ (1)

Indices $i$ and $j$ run over the different parton species in each of the
incident protons. The parton-level subprocess cross sections
$d\hat{\sigma}_{(ij\rightarrow ab)}(x_{1},x_{1}^{\prime},\mu)$ are functions
of the fractional partonic longitudinal momenta $x_{1}$ and $x_{1}^{\prime}$
from each of the incident hadrons and of the partonic
factorization/renormalization scale $\mu$. The parton distribution functions
$f^{i}_{p}(x_{1},\mu)$ express the probability that parton $i$ is found with
fractional longitudinal momentum $x_{1}$ at scale $\mu$ in the proton; they
are integrated over the intrinsic transverse momentum (equivalently, impact
parameter) carried by the parton in the parent hadron.

A formal theoretical treatment of double parton scattering would begin with a
discussion of the hadronic matrix element of four field operators and an
explicit operator definition of two-parton correlation functions. This
procedure would lead to a decomposition of the hadronic matrix element into
non-perturbative two-parton distribution functions and the corresponding hard
partonic cross sections for $\hat{\sigma}(ijkl\rightarrow abcd)$. An operator
definition of two-parton correlation functions may be found in Ref.
Mueller:1985wy  where the two-parton correlation function is reduced to a
product of single parton distributions. An explicit operator definition of
two-parton distributions with different values of the two fractional momenta
$x_{1}$ and $x_{2}$ is presented in Ref. Guo:1997it , along with a model for
the two-parton distributions in terms of normal parton distributions. In this
paper, we follow a phenomenological approach along lines similar to Refs.
Goebel:1979mi ; Paver:1982yp ; Humpert:1983pw ; Mekhfi:1983az ; Humpert:1984ay
; Ametller:1985tp ; Halzen:1986ue ; Mangano:1988sq ; Godbole:1989ti ;
Drees:1996rw ; Eboli:1997sv ; Yuan:1997tr ; Calucci:1997uw ; DelFabbro:1999tf
; Kulesza:1999zh ; Korotkikh:2004bz ; Cattaruzza:2005nu ; Hussein:2006xr ;
Maina:2009sj ; Domdey:2009bg ; d'Enterria:2009hd ; Gaunt:2009re .

In a double parton process, partons $i$ and $k$ are both active in a given
incident proton. We require the joint probability that parton $k$ carries
fractional momentum $x_{2}$, given that parton $i$ carries fractional momentum
$x_{1}$. In general, this joint probability
$H^{i,k}(x_{1},x_{2},\mu_{A},\mu_{B})$ should also depend on the intrinsic
transverse momenta $k_{T,i}$ and $k_{T,k}$ of the two partons (or,
equivalently, their impact parameters). The hard scales $\mu_{A}$ and
$\mu_{B}$ are characteristic of the two hard subprocesses in which partons $i$
and $k$ participate. In the sections below, we discuss the choice we make of
the hard-scale and do not explore in this paper theoretical uncertainties
associated with higher-order perturbative contributions.

In contrast to single parton distributions functions $f^{i}_{p}(x_{1},\mu)$
for which global analyses have produced detailed information, very little is
known phenomenologically about the magnitude and functional dependences of
joint probabilities $H^{i,k}(x_{1},x_{2},\mu_{A},\mu_{B})$. A common
assumption made in estimates of double parton rates is to ignore possibly
strong correlations in longitudinal momentum and to use the approximation

$\displaystyle
H^{i,k}_{p}(x_{1},x_{2},\mu_{A},\mu_{B})=f^{i}_{p}(x_{1},\mu_{A})f^{k}_{p}(x_{2},\mu_{B}).$
(2)

For reasons of energy-momentum conservation, if not dynamics, the simple
factorized form of Eq. (2) cannot be true for all values of the fractional
momenta $x$. The values of $x_{2}$ available to the second interaction are
always limited by the values of $x_{1}$ in the initial interaction since
$x_{1}+x_{2}\leq 1$. The approximation certainly fails even at the kinematic
level if both partons carry a substantial fraction of the momentum of the
parent hadron. However, it may be adequate for applications in which the
values of $x_{1}$ and $x_{2}$ are small. We remark that the momentum integral

$\displaystyle\sum_{i,k}\int
x_{1}x_{2}H^{i,k}_{p}(x_{1},x_{2},\mu_{A},\mu_{B})dx_{1}dx_{2}=1,$ (3)

as long as we can run the upper limits of the $x_{1}$ and $x_{2}$ integrations
to $1$, independently. The large phase space at the LHC may make it possible
to explore dynamic correlations that break Eq. (2).

In Fig. 2, for the region of phase space of interest to us, we show the
contributions to the $b\bar{b}jj$ cross section as a function of $x$ from both
DPS and SPS, after minimal acceptance cuts are imposed (Sec. II). The center-
of-mass energy is $\sqrt{s}=10$ TeV. It is evident that the majority of DPS
events are associated with low $x$ values, in essence never exceeding $0.2$.
The momentum carried off by the beam remnant is $(1-x_{1}-x_{2})$ in DPS and
$(1-x)$ in SPS. The results in Fig. 2 show that this remnant momentum is not
too different in DPS and SPS. Thus, the use of Eq. (2) in calculations of
event rates at the LHC appears adequate as a good first approximation. While
available Tevatron data on double parton scattering Abe:1997xk ; D0:2009  are
insensitive to possible correlations in $x$, the greater dynamic range at the
LHC may make it possible to observe them. 111As emphasized in Refs.
Korotkikh:2004bz ; Gaunt:2009re , even if the approximation in Eq. (2) holds
at one hard scale, evolution of the parton densities with $\mu$ will induce
violations at larger scales.

Figure 2: Values of the parton longitudinal momentum fractions $x$ in the DPS
and SPS events. Most DPS events have low $x$ values. The events used for this
plot include the requirements $n_{\rm jet}=4$, $n_{\rm btag}=2$, and the
threshold cuts discussed in Sec. II.

Assuming next that the two subprocesses $A(i~{}j\rightarrow a~{}b)$ and
$B(k~{}l\rightarrow c~{}d)$ are dynamically uncorrelated, we express the
double parton scattering differential cross section as:

$\displaystyle d\sigma^{DPS}=\dfrac{m}{2\sigma_{\rm eff}}\sum_{i,j,k,l}\int
H^{ik}_{p}(x_{1},x_{2},\mu_{A},\mu_{B})H_{p}^{jl}(x_{1}^{\prime},x_{2}^{\prime},\mu_{A},\mu_{B})$
(4) $\displaystyle\times
d\hat{\sigma}^{A}_{ij}(x_{1},x_{1}^{\prime},\mu_{A})d\hat{\sigma}^{B}_{kl}(x_{2},x_{2}^{\prime},\mu_{B})dx_{1}dx_{2}dx_{1}^{\prime}dx_{2}^{\prime}.$

The symmetry factor $m$ is $1$ if the two hard-scattering subprocesses are
identical and is $2$ otherwise. In the denominator, there is a factor
$\sigma_{\rm eff}$ with the dimensions of a cross section. Given that one
hard-scatter has taken place, $\sigma_{\rm eff}$ measures the size of the
partonic core in which the flux of accompanying short-distance partons is
confined. It should be at most proportional to the transverse size of a
proton. For the first process of interest in this paper, $pp\rightarrow
b\bar{b}j_{1}j_{2}$, Eq. (4) reduces to

$\displaystyle d\sigma^{DPS}(pp\rightarrow
b\bar{b}j_{1}j_{2}X)=\dfrac{d\sigma^{SPS}(pp\rightarrow
b\bar{b}X)d\sigma^{SPS}(pp\rightarrow j_{1}j_{2}X)}{\sigma_{\rm eff}}.$ (5)

Tevatron collider data Abe:1997xk ; D0:2009  yield values in the range
$\sigma_{\rm eff}\sim 12$ mb. We use this value for the estimates we make, but
we emphasize that the goal should be to make an empirical determination of its
value at LHC energies.

In Sec. II, we present our calculation of the double parton and the single
parton contributions to $p~{}p\rightarrow b~{}\bar{b}~{}j_{1}~{}j_{2}~{}X$. We
identify variables that discriminate the two contributions quite well. In Sec.
III, we treat the double parton and the single parton contributions to $4$ jet
production, again finding that good separation is possible despite the
combinatorial uncertainty in the pairing of jets. We show in both cases that
the double parton contribution falls off significantly more rapidly with
$p_{T}^{j1}$, the transverse momentum of the leading jet. For the value of
$\sigma_{\rm eff}\sim 12$ mb and the cuts that we use, we find that, in the
region in which it is most identifiable, double parton scattering is dominant
for $p_{T}^{j1}<30$ GeV in $b~{}\bar{b}~{}j_{1}~{}j_{2}$ at LHC energies, and
$p_{T}^{j1}<50$ GeV in $4$ jet production. Our conclusions are found in Sec.
IV.

##  II Heavy quark pair and jet pair production in QCD.

In this section, we describe the calculation of the DPS and SPS event rates
for $b\bar{b}jj$ production at the LHC. For our purposes, light jets (denoted
by $j$) are assumed to originate only from gluons or one of the four lighter
quarks ($u,d,s$ or $c$) and, as stated above, we perform all calculations for
the LHC with a center-of-mass energy of $\sqrt{s}=10$ TeV. Event rates are
quoted for 10 pb-1 of data.

###  II.1 Outline of the method

The prediction for the DPS event rate is based on the assumption that the two
partonic interactions which produce the $b\bar{b}$ and $jj$ systems occur
independently (as expressed in Eq. (4)). At leading order, the only
contribution is:

$(ij\rightarrow b\bar{b})\otimes(kl\rightarrow jj)$ (6)

where the symbol $\otimes$ denotes the combination of one event each from the
$b\bar{b}$ and the $jj$ final states. In an attempt to model some of the
effects expected from initial- and final-state radiation, we also account for
the possibility of an additional jet which is undetected because it is either
too soft or outside of the accepted rapidity range. Thus, we include several
other contributions to the DPS event:

$\displaystyle b\bar{b}(j)\otimes
jj\,\,\,,\,\,\,b\bar{b}j\otimes(j)j\,\,\,,\,\,\,b\bar{b}j\otimes j(j)$ (7)
$\displaystyle b\bar{b}\otimes(j)jj\,\,\,,\,\,\,b\bar{b}\otimes
j(j)j\,\,\,,\,\,\,b\bar{b}\otimes jj(j)\,,$ (8)

where the parentheses surrounding a jet indicate that it is undetected. We
compute processes such as $jj(j)$ and $b\bar{b}(j)$ at LO as 3 parton final-
state processes.

The 2 to 3 parton amplitudes for $b\bar{b}(j)$ [and $jj(j)$] diverge as the
undetected jet $(j)$ becomes soft or collinear to one of the other final state
partons or to an initial parton. The divergences are removed in a full next-
to-leading order (NLO) treatment, in which real emission and virtual (loop)
contributions are incorporated, and the finite $b\bar{b}$, $b\bar{b}(j)$, and
$b\bar{b}j$ contributions are present with proper relative normalization. In
the LO parton level simulations done in this paper, we employ a cut at the
generator level to remove the divergences. All the final state objects in the
processes listed above are required to have transverse momentum $p_{T}\geq 20$
GeV. In this fashion, we model some aspects of the expected momentum imbalance
between the $b$ and $\bar{b}$ arising from the 2 to 3 process $ij\rightarrow
b\bar{b}j$, but we cannot claim to include the relative normalization between
the $b\bar{b}$ and $b\bar{b}j$ contributions that would result from a full NLO
treatment. We leave a complete NLO analysis for future work.

The SPS cross section is computed according to Eq. (1). It receives
contributions at lowest order from the 2 parton to 4 jet final state process:

$ij\rightarrow b\bar{b}jj\,,$ (9)

and, in the case where a jet is undetected, from the 5-jet final states
(computed at LO):

$b\bar{b}(j)jj\,\,\,,\,\,\,b\bar{b}j(j)j\,\,\,,\,\,\,b\bar{b}jj(j)\,.$ (10)

We also investigate the possibility of $jjjj$ and $jjjj(j)$ final state
contributions to the SPS cross section where two of the jets “fake” $b$ jets.
We find that the effects from these final states are subdominant compared to
the processes listed in Eqs. (9) and (10).

In our numerical analysis, we use the leading-order CTEQ6L1 parton
distribution functions (PDFs) Pumplin:2002vw  to compute both DPS and SPS
cross sections, and we evaluate all cross sections using one-loop evolution of
$\alpha_{s}(\mu)$. For the renormalization and factorization scales, we choose
the dynamic scale:

$\mu^{2}=\sum_{i}p_{T,i}^{2}+m_{i}^{2}\,,$ (11)

where $p_{T,i}$ is the transverse momentum of the $i^{th}$ jet and $m_{i}=0$
($m_{i}=4.7$ GeV) for light (bottom) jets. In the case of roughly equal values
of the transverse momenta $p_{T,i}$, Eq. (11) yields $\mu\sim 2p_{T}$ in SPS
and $\mu\sim\sqrt{2}p_{T}$ in DPS. At LO there is no obviously “right” hard
scale, and the choice in Eq. (11) seems as good as any other.

The DPS events are generated as two separate sets of events with
Madgraph/Madevent Maltoni:2002qb  and then combined as described above. For
example, at leading order, we generate events separately for $pp\to b\bar{b}X$
and $pp\to jjX$, and these events are then combined as indicated in Eq. (6).
To increase the speed of the simulations, the SPS events are generated with
Alpgen Mangano:2002ea  since the SPS processes of interest are hard-coded in
Alpgen, which contains more compact expressions for the squared-matrix-
elements than Madgraph.

The events accepted after generation are required to have 4 jets $n_{\rm
jet}=4$ with 2 of these tagged as $b$’s $n_{\rm btag}=2$. At the generator
level, all the final state objects in the processes listed in Eq. (6) through
Eq. (10) must have transverse momentum $p_{T}\geq 20$ GeV, as mentioned above.
Furthermore, at the analysis level, all events (DPS and SPS) are required to
pass the following acceptance cuts:

$\displaystyle p_{T,j}$ $\displaystyle\geq$ $\displaystyle
25\,\,\,\mbox{GeV},\,\,\,|\eta_{j}|\leq 2.5$ (12) $\displaystyle p_{T,b}$
$\displaystyle\geq$ $\displaystyle 25\,\,\,\mbox{GeV},\,\,\,|\eta_{b}|\leq
2.5$ (13) $\displaystyle\Delta R_{jj}$ $\displaystyle\geq$ $\displaystyle
0.4,\,\,\,\Delta R_{bb}\geq 0.4$ (14)

where $\eta_{i}$ is the jet’s pseudorapidity, and $\Delta R_{ij}$ is the
separation in the azimuthal angle ($\phi$) - pseudorapidity plane between jets
$i$ and $j$:

$\Delta R_{ij}=\sqrt{(\eta_{i}-\eta_{j})^{2}+(\phi_{i}-\phi_{j})^{2}}\,.$ (15)

We model detector resolution effects by smearing the final state energy
according to:

${\delta E\over E}={a\over\sqrt{E/\rm{GeV}}}\oplus b,$ (16)

where we take $a=50\%$ and $b=3\%$ for jets. To account for $b$ jet tagging
efficiencies, we assume a $b$-tagging rate of 60% for $b$-quarks with
$p_{T}>20\text{ GeV}$ and $|\eta_{b}|<2.0$. We also apply a mistagging rate
for charm-quarks as:

$\epsilon_{c\to b}=10\%\quad\quad\text{ for }p_{T}(c)>50\text{ GeV}\\\ $ (17)

while the mistagging rate for a light quark is:

$\displaystyle\epsilon_{u,d,s,g\to b}$ $\displaystyle=0.67\%\quad\quad\text{
for }$ $\displaystyle p_{T}(j)<100\text{ GeV}$ (18)
$\displaystyle\epsilon_{u,d,s,g\to b}$ $\displaystyle=2\%\quad\quad\quad\text{
for }$ $\displaystyle p_{T}(j)>250\text{ GeV}.$ (19)

Over the range $100\text{ GeV}<p_{T}(j)<250\text{ GeV}$, we linearly
interpolate the fake rates given above Baer:2007ya .

###  II.2 Properties of SPS and DPS in $b~{}\bar{b}~{}j~{}j$

Having detailed the calculation of the $b\bar{b}jj$ event rates from DPS and
SPS, we now discuss some of the distinguishing characteristics of the two
contributions. First, however, it is important to check that our simulations
of DPS events are not introducing an artificial correlation between the
$b\bar{b}$ and $jj$ final states. We do this by inspecting the angle $\Phi$
between the plane defined by the $b\bar{b}$ system and the plane defined by
the $jj$ system. If the two scattering processes $ij\rightarrow b\bar{b}$ and
$kl\rightarrow jj$ which produce the DPS final state are truly independent,
one would expect to see a flat distribution in the angle $\Phi$. By contrast,
many diagrams, including some with non-trivial spin correlations, contribute
to the 2 parton to 4 parton final state in SPS, and naively one would expect
some correlation between the two planes. To avoid possible effects from
boosting to the lab frame, we define the two planes in the partonic center-of-
mass frame.

We specify the planes by using the three-momenta of the outgoing jets. Then,
the angle between the two planes defined by the $jj$ and $b\bar{b}$ systems
is:

$\cos\Phi=\hat{n}_{3}(j_{1},j_{2})\cdot\hat{n}_{3}(b_{1},b_{2}),$ (20)

where $\hat{n}_{3}(x,y)$ is the unit three-vector normal to the plane defined
by the $x-y$ system.

The normal is undefined when $j_{1}$ and $j_{2}$ are back-to-back or $b_{1}$
and $b_{2}$ are back-to-back, as occurs in a large fraction of the DPS events.
Therefore, when $\cos\phi_{(x,y)}<-0.9$, we use a different procedure. We use
the three-momentum of one of the incoming partons along with the three-
momentum of one of the outgoing $b$ quarks to define the $b\bar{b}$ plane. Let
$q_{b}$ be the three-momentum of an incoming parton, and $p_{b}$ be the three-
momentum of the final-state $b$ (or $\bar{b}$) quark. We then define
$\phi_{p_{b},q_{b}}$ to be the azimuthal angle of the three-vector normal to
the $q_{b}-p_{b}$ plane. Note that we use $\phi$ here since the normal to any
three-vector and the beam-line will be transverse to the beam-line (not the
case in the SPS process). In this way, the jet which is not used to define the
plane is guaranteed to lie in the plane. The plane for the $jj$ system is
defined in an analogous manner. Finally, the angle between the planes is then:

$\Phi=|\phi_{p_{j},q_{j}}-\phi_{p_{b},q_{b}}|\,.$ (21)

In Fig. 3, we display the number of events as a function of the angle between
the two planes. There is an evident correlation between the two planes in SPS,
while the distribution is flat in DPS, indicative that the two planes are
uncorrelated.

Figure 3: Event rate as a function of the angle between the two planes defined
by the $b\bar{b}$ and $jj$ systems. In SPS events, there is a correlation
among the planes which is absent for DPS events.

Another interesting difference between DPS and SPS is the behavior of event
rates as a function of transverse momentum. As an example of this, in Fig. 4,
we show the transverse momentum distribution for the leading jet (either a $b$
or light $j$) for both DPS and SPS. Several characteristics are evident.
First, SPS produces a relatively hard spectrum, and for the value of
$\sigma_{\rm eff}$ and the cuts that we use, we see that SPS tends to dominate
over the full range of transverse momentum considered. On the other hand, DPS
produces a much softer spectrum which (up to issues of normalization in the
form of $\sigma_{\rm eff}$) can dominate at small values of transverse
momentum. The cross-over between the two contributions to the total event rate
is $\sim 30$ GeV for the acceptance cuts considered here. A smaller (larger)
value of $\sigma_{\rm eff}$ would move the cross-over to a larger (smaller)
value of the transverse momentum $p_{T}^{j1}$ of the leading jet.

Figure 4: The transverse momentum $p_{T}$ distribution of the leading jet in
$jjb\bar{b}$ after minimal cuts.

###  II.3 Distinguishing variables

We turn next to the search for variables that may allow for a clear separation
of the DPS and SPS contributions. Since the topology of the DPS events
includes two $2\to 2$ hard scattering events, the two pairs of jet objects are
roughly back-to-back. We expect the azimuthal angle between the pairs of jets
corresponding to each hard scattering event to be strongly peaked near
$\Delta\phi_{jj}\sim\Delta\phi_{bb}\sim\pi$. Real radiation of an additional
jet, where the extra jet is missed because it fails the threshold or
acceptance cuts, allows smaller values of $\Delta\phi_{jj}$. The relevant
distribution is shown for light jets (non $b$-tagged) in Fig. 5a. There is a
clear peak near $\Delta\phi_{jj}=\pi$ for DPS events, while the events are
more broadly distributed in SPS events. The secondary peak near small
$\Delta\phi_{jj}$ arises from gluon splitting which typically produces nearly
collinear jets. The suppression at still lower $\Delta\phi_{jj}$ comes from
the isolation cut $\Delta R_{jj}>0.4$.

Figure 5: (a) The difference $\Delta\phi$ in the azimuthal angles of light jet
pairs for DPS and both SPS+DPS events. The dijet pairs are back-to-back in DPS
events. (b) The variable $S_{\phi}$ for DPS and SPS+DPS events provides a
stronger separation of the underlying DPS events from the total sample when
compared to $\Delta\phi$ for any pair.

The separation of DPS events from SPS events becomes more pronounced if
information is used from both the $b\bar{b}$ and $jj$ systems. As an example,
we consider the distribution built from a combination of the azimuthal angle
separations of both $jj$ and $b\bar{b}$ pairs, using a variable adopted from
Ref. D0:2009 :

$S_{\phi}={1\over\sqrt{2}}\sqrt{\Delta\phi(b_{1},b_{2})^{2}+\Delta\phi(j_{1},j_{2})^{2}}.$
(22)

In Fig. 5b, we present a distribution in $S_{\phi}$ for both DPS and SPS+DPS
events. Again, as in the case of the $\Delta\phi$ distribution, we see that
the SPS events are broadly distributed across the allowed range of $S_{\phi}$.
However, the combined information from both the $b\bar{b}$ and $jj$ systems
shows that the DPS events produce a sharp and substantial peak near
$S_{\phi}\simeq\pi$ which is well-separated from the total sample.

The narrow peaks near $\Delta\phi_{jj}=\pi$ in Fig. 5a and near $S_{\phi}=1$
in Fig. 5b will be smeared somewhat once soft QCD radiation and other higher-
order terms are included in the calculation.

Another possibility for discerning DPS is the use of the total transverse
momentum of both the $b\bar{b}$ and $jj$ systems. At lowest order for a $2\to
2$ process, the vector sum of the transverse momenta of the final state pair
vanishes. In reality, radiation and momentum mismeasurement smear the expected
peak near zero. Nevertheless, we still expect DPS events to show a
distribution in the transverse momenta of the jet pairs that is reasonably
well-balanced. To encapsulate this expectation for both light jet pairs and
$b$-tagged pairs, we use the variable D0:2009 :

$S_{p_{T}}^{\prime}={1\over\sqrt{2}}\sqrt{\left({|p_{T}(b_{1},b_{2})|\over|p_{T}(b_{1})|+|p_{T}(b_{2})|}\right)^{2}+\left({|p_{T}(j_{1},j_{2})|\over|p_{T}(j_{1})|+|p_{T}(j_{2})|}\right)^{2}}.$
(23)

Here $p_{T}(b_{1},b_{2})$ is the vector sum of the transverse momenta of the
two final state $b$ jets, and $p_{T}(j_{1},j_{2})$ is the vector sum of the
transverse momenta of the two (non $b$) jets.

The distribution in $S_{p_{T}}^{\prime}$ is shown in Fig. 6. As expected, we
observe that the DPS events are peaked near $S_{p_{T}}^{\prime}\sim 0$ and are
well-separated from the total sample. The SPS events, on the other hand, tend
to be far from a back-to-back configuration and, in fact, are peaked near
$S_{p_{T}}^{\prime}\sim 1$. This behavior of the SPS events is presumably
related to the fact that a large number of the $b\bar{b}$ or $jj$ pairs arise
from gluon splitting which yields a large $p_{T}$ imbalance and, thus, larger
values of $S_{p_{T}}^{\prime}$.

Figure 6: Distribution of events in $S_{p_{T}}^{\prime}$ for the DPS and SPS
samples. Due to the back-to-back nature of the $2\to 2$ events in DPS
scattering, the transverse momenta of the jet pair and of the $b$-tagged jet
pair are small, resulting in a small value of $S_{p_{T}}^{\prime}$. In (a) we
show the $S_{p_{T}}^{\prime}$ distribution for our standard cuts, and in (b)
we increase the cut on the transverse momentum of the leading jet,
$p_{T}^{j1}>40$ GeV. The fraction of DPS events in the whole sample decreases
with increasing $p_{T}^{j1}$.

In this subsection, we find that extraction of the DPS “signal” for
$b\bar{b}jj$ production from the SPS “background” can be enhanced by combining
information from both $b\bar{b}$ and $jj$ systems. Our simulations suggest
that the variable $S_{p_{T}}^{\prime}$ may be a more effective discriminator
than $S_{\phi}$. However, given the leading order nature of our calculations
and the absence of smearing associated with initial state soft radiation, this
picture may change and a variable such as $S_{\phi}$ (or some other variable)
may become a clearer signal of DPS at the LHC. Realistically, it would be
valuable to study both distributions once LHC data are available in order to
determine which is more instructive. In the following, we use the clear
separation shown in Fig. 6 in our exploration of the distinct properties of
DPS and SPS events.

###  II.4 Two-dimensional distributions

The evidence in Fig. 5 and Fig. 6 for distinct regions of DPS dominance
prompts the search for greater discrimination in a plane represented by a two
dimensional distribution of one variable against another. We examined scatter
plots involving the inter-plane angle $\Phi$, the jet-jet azimuthal angle
difference $\Delta\phi_{jj}$, $S_{\phi}$, and $S^{\prime}_{p_{T}}$. Strong
kinematic correlations are evident in the plot of $S_{\phi}$ vs.
$S^{\prime}_{p_{T}}$ at the level of our leading order calculation, and we
observe no additional separation of DPS and SPS beyond that evident in Figs. 5
and 6. Likewise, there are strong correlations between $\Delta\phi_{jj}$ and
$S_{\phi}$.

One scatter plot with interesting features is displayed in Fig. 7. The DPS
events are seen to be clustered near $S^{\prime}_{p_{T}}=0$ and are uniformly
distributed in $\Phi$. The SPS events peak toward $S^{\prime}_{p_{T}}=1$ and
show a roughly $\sin\Phi$ character. While already evident in Figs. 3 and 6,
these two features are more apparent in the scatter plot Fig. 7. Moreover, the
scatter plot shows a valley of relatively low density between
$S^{\prime}_{p_{T}}\sim 0.1$ and $\sim 0.4$. In an experimental one-
dimensional $\Phi$ distribution such as Fig. 3, one would see the sum of the
DPS and SPS contributions. If structure is seen in data similar to that shown
in the scatter plot Fig. 7, one could make a cut at $S^{\prime}_{p_{T}}<0.1$
or $0.2$ and verify whether the experimental distribution in $\Phi$ is flat as
expected for DPS events.

Figure 7: Two-dimensional distribution of events in the variables $\Phi$ and
$S_{p_{T}}^{\prime}$ for the DPS and SPS samples.

In Fig. 4, we show that DPS produces a softer transverse momentum distribution
for the leading jet (either a $b$ or light $j$). In data one would see only
the sum of the DPS and SPS components in a plot like Fig. 4. A scatter plot of
$S_{p_{T}}^{\prime}$ vs. the transverse momentum of the leading jet motivates
an empirical separation of the two components. In Figs. 6(a) and 6(b) we
compare the $S_{p_{T}}^{\prime}$ distributions for two different selections on
the transverse momentum $p_{T}^{j1}$ of the leading jet in the $b\bar{b}jj$
sample. This comparison of the distributions confirms that events in the DPS
region, defined empirically by the region $S_{p_{T}}^{\prime}<0.1$ or $0.2$,
fall off more steeply with $p_{T}^{j1}$ than the rest of the sample. It will
be important and interesting to see whether the selection
$S_{p_{T}}^{\prime}<0.1$ or $0.2$ in LHC data also produces events that show a
more rapid decrease with $p_{T}^{j1}$.

The leading-jet transverse momentum distributions are shown in Figs. 8(a) and
8(b) for two different cuts on $S_{p_{T}}^{\prime}$. In both cases, we see
that the SPS sample has a broader distribution in $p_{T}^{j1}$ and that the
DPS sample dominates for small enough values of $p_{T}^{j1}$. For our chosen
value of $\sigma_{\rm eff}\sim 12$ mb, and for cuts we employ, the crossover
points are roughly $80$ GeV for $S_{p_{T}}^{\prime}<0.2$ and $40$ GeV for
$S_{p_{T}}^{\prime}<0.4$.

Figure 8: The distribution in the transverse momentum of the leading jet
$p_{T}^{j1}$ for (a) $S_{p_{T}}^{\prime}<0.2$ and (b)
$S_{p_{T}}^{\prime}<0.4$. As the signal region becomes more dominated by SPS
events (i.e. moving from (a) to (b)), the resulting distribution becomes
harder and shifts the SPS-DPS cross-over from $\sim 80$ GeV to $\sim 40$ GeV.

##  III Four Jet Production

In addition to $b\bar{b}jj$, we can also ask how important DPS can be for a
generic $4j$ final-state, where none of the jets are $b$-tagged. In this
section, we describe our calculation of the double parton scattering and the
single parton scattering contributions to the production of a $4j$ final
state, for which the cross section is larger. Our exposition can be brief
since we repeat the procedure described in some detail in Sec. II.

###  III.1 Outline of the method

The DPS process for $4j$ production is topologically equivalent to
$b\bar{b}jj$. However, in the $4j$ system, we lose the $b$-tagging ability
that reduces the combinatorial background in $b\bar{b}jj$, and the prospects
for isolating and measuring DPS over the SPS background may appear less
promising. Fortunately, in going from the $b\bar{b}$ subprocess to the $jj$
subprocess, a much larger DPS rate is possible due to the much larger cross
section for $jj$ production. As we show below, we find that the DPS signature
can be extracted in this $4j$ mode as well.

The DPS cross section for $4j$ production receives contributions from the
following sub-processes at the lowest order:

$jj\otimes jj\,\,\,,\,\,\,b\bar{b}\otimes jj,$ (24)

where both $b$-quarks fail the $b$-tag. We do not include the $b\bar{b}\otimes
b\bar{b}$ process due to its relatively small rate ($\sim 0.14$ nb). This rate
is further reduced by requiring no $b$-tags, yielding roughly 40 events in the
10 pb-1 of luminosity assumed here.

Following Sec. II, we account for the possibility of an additional jet which
is undetected because it is too soft or outside of the accepted rapidity
range. Thus, we include several other contributions to the DPS cross section:

$\displaystyle jjj\otimes(j)j\,\,\,,\,\,\,jj(j)\otimes jj\,,$ (25)
$\displaystyle b\bar{b}j\otimes j(j)\,\,\,,\,\,\,b\bar{b}(j)\otimes jj,$ (26)
$\displaystyle b\bar{b}\otimes j(j)j\,\,\,,\,\,\,b(\bar{b})\otimes
jjj\,\,\,,\,\,\,(b)\bar{b}\otimes jjj\,.$ (27)

where the parentheses surrounding a jet signify that it is not detected.

The SPS cross section receives contributions at lowest order from the final
state:

$jjjj\,\,\,,b\bar{b}jj\,,$ (28)

where both $b$-quarks fail the $b$-tag, and, in the case where a jet is not
detected, from the final states:

$b\bar{b}(j)jj\,\,\,,\,\,\,(b)\bar{b}jjj\,\,\,,\,\,\,b(\bar{b})jjj\,\,\,,\,\,\,(j)jjjj\,.$
(29)

We refer to Sec. II for the specification of acceptance cuts and detector
resolution, and for our treatment of the potential divergences present in the
amplitudes for the processes in Eqs. (24)-(29).

###  III.2 Results

Similar to the $b\bar{b}jj$ process, the leading jet in the $4j$ DPS sample is
typically softer than in the SPS channels (see Fig. 9). In this case, again
using $\sigma_{\rm eff}=12$ mb, we find that the cross-over between DPS and
SPS dominance occurs near $p_{T}\simeq 50$ GeV, higher than in the
$b\bar{b}jj$ case shown in Fig. 4.

Figure 9: As in Fig. 4, but for $4j$ events. Similar to the $b\bar{b}jj$
sample, the SPS sample exhibits a harder $p_{T}$ spectrum.

Improvement in the separation between DPS and SPS in the $4j$ case can be
achieved with an analogous version of the $S_{p_{T}}^{\prime}$ variable
introduced in Eq. (23):

$S_{p_{T}}^{\prime}={1\over\sqrt{2}}\sqrt{\left({|p_{T}(j_{a},j_{b})|\over|p_{T}(j_{a})|+|p_{T}(j_{b})|}\right)^{2}+\left({|p_{T}(j_{c},j_{d})|\over|p_{T}(j_{c})|+|p_{T}(j_{d})|}\right)^{2}}.$
(30)

Here $p_{T}(j_{a},j_{b})$ is the vector sum of the transverse momenta of two
final state jets, $a$ and $b$, chosen among the four. The remaining $c$ and
$d$ jets are then fixed. This choice is unique if a separation of the two hard
interactions is possible. In the $b\bar{b}jj$ system, the separation into the
$b\bar{b}$ and $jj$ subsystems via $b$-tagging removed most of the degeneracy
(some degeneracy still remained via tagging efficiencies or light jet
mistagging). In the $4j$ system, the degeneracy can at first glance be
problematic as there are 3 possible pairings of the four jets.

Figure 10: The democratic $S_{p_{T}}^{\prime}$ distribution for $4j$ events
shows much more combinatorial background than in the $b\bar{b}jj$ events. Even
after accepting two mis-matched jet pairs, we see that the DPS and SPS samples
can still be separated well.

One might be tempted to take the pairing of jets which minimizes the value of
$S_{p_{T}}^{\prime}$. Unfortunately, this choice places a bias on the
distribution that makes it potentially problematic to trust the
discrimination. Instead, to construct $S_{p_{T}}^{\prime}$ we take all three
combinations of pairings, which includes one “correct” pairing and two
incorrect pairings in the DPS process. This “democratic” $S_{p_{T}}^{\prime}$
distribution is shown in Fig. 10 and is re-weighted by 1/3 for proper
normalization. As in the $b\bar{b}jj$ case, we see that the DPS distribution
peaks near $S_{p_{T}}^{\prime}\sim 0$, indicative that two back-to-back hard
interactions are present. In addition to this expected feature, we also see a
continuum that extends above $S_{p_{T}}^{\prime}\sim 0.1$, associated with the
wrong combination taken in the democratic approach. In Fig. 10 we see that DPS
produces a secondary peak at $S_{p_{T}}^{\prime}\sim 1$, not present in the
$b\bar{b}jj$ case in Fig. 6. It appears to arise from the wrong pairings of
jets associated with the combinatorial background. In these instances, the
wrong combination of two jets that are close together in $\Delta R$, meaning
that their momenta are aligned, can maximize the value of
$S_{p_{T}}^{\prime}$. Overall, we see that the DPS peak near
$S_{p_{T}}^{\prime}=0$ provides a good means to separate DPS events from SPS
events.

Figure 11: As in Fig. 8, but for $4j$ events with (a) democratic
$S_{p_{T}}^{\prime}<0.2$ and (b) democratic $S_{p_{T}}^{\prime}<0.4$. As in
$b\bar{b}jj$ events, as one increases the cut on $S_{p_{T}}^{\prime}$, the SPS
fraction increases and the total distribution is harder.

As in the $b\bar{b}jj$ case, we inspect the distribution in the $p_{T}$ of the
leading jet after cuts on the $S_{p_{T}}^{\prime}$ variable. Since there are
three jet pairings per event, we now require that at least one of the three
pairings has $S_{p_{T}}^{\prime}$ in the given window. Due to this softer
constraint, the hardening of the $p_{T}$ spectrum of the leading jet is less
dramatic than in the $b\bar{b}jj$ case (e.g. compare Figs. 8 and 11). The
crossover of the SPS and DPS contributions occurs near $80$ GeV for
$S_{p_{T}}^{\prime}<0.2$ and near $50$ GeV for $S_{p_{T}}^{\prime}<0.4$

##  IV Discussion and conclusions

Our goal is to develop a method to search for a double parton scattering
contribution in the $b~{}\bar{b}~{}j~{}j$ and 4 jet final states at LHC
energies and to measure the magnitude of its contribution relative to the
single parton contribution to the same final states. Based on our parton level
simulations, we find that variables such as $S_{p_{T}}^{\prime}$ and
$S_{\phi}$ that take into account information from the entire final state,
thereby including both of the hard subprocesses in DPS, are more effective at
discrimination than variables such as $\Delta\phi_{jj}$ that reflect only a
subset of the final-state. The enhancement at low values of
$S_{p_{T}}^{\prime}$ shown in Figs. 6, 7 and 10 provides a good signature for
the presence of double parton scattering. We urge experimenters to search for
such a concentration of events in data at the LHC. Having found this
enhancement, we then suggest that the magnitude of this peak be examined as a
function of the transverse momentum $p_{T}^{j1}$ of the leading jet in the
event sample. The double parton scattering contribution in the peak region
should fall off more rapidly with $p_{T}^{j1}$ than the rest of the sample.
The distribution of events in the region of small values of
$S_{p_{T}}^{\prime}$ should also be examined as a function of the inter-plane
angle $\Phi$ to see whether the flat behavior is seen, as expected for two
independent production processes. Once these characteristics of double parton
scattering are established, the data can be used to determine the effective
normalization $\sigma_{\rm eff}$, defined and discussed in the Introduction.
It will be interesting to see whether the values extracted for $\sigma_{\rm
eff}$ are about the same in the $b~{}\bar{b}~{}j~{}j$ and 4 jet final states
and how they compare with values measured at the Fermilab Tevatron.

Once double parton scattering is established in data, and $\sigma_{\rm eff}$
is determined, in a relatively clean process such as $b\bar{b}jj$, double
parton contributions to a wide range of other processes can be computed with
more certainty about their expected rates at LHC energies. To be sure, given
the approximations described in the Introduction, some variation in the values
of $\sigma_{\rm eff}$ might be expected and appropriate for different
processes and in different kinematic regions. The connection of $\sigma_{\rm
eff}$ with the effective size of the hard-scattering core of the proton may
mean that $\sigma_{\rm eff}$ will have different values for $gg$, $qq$, and
$q\bar{q}$ scattering.

There are several avenues for future work. Of great importance is the proper
inclusion of next-to-leading order contributions NLO . They are needed to make
more robust predictions of the relative normalization of the DPS and SPS
contributions, of the shape of the $p_{T}$ distribution of the leading jet,
and for proper softening of the sharp peaks seen near $S_{p_{T}}^{\prime}=1$
in Figs. 6 and 10, and near $S_{\phi}=\pi$ in Fig. 5b.

It will also be important to develop joint probabilities
$H^{i,k}(x_{1},x_{2},\mu_{A},\mu_{B})$ that are more sophisticated
theoretically than the first approximation represented by Eq. (2) in which
parton-parton correlations are absent. A valuable development in this
direction are the studies presented in Refs. Korotkikh:2004bz ; Gaunt:2009re .

Double parton contributions are potentially relevant for a wide range of
standard model processes, many already considered in the literature
Goebel:1979mi ; Paver:1982yp ; Humpert:1983pw ; Mekhfi:1983az ; Humpert:1984ay
; Ametller:1985tp ; Halzen:1986ue ; Mangano:1988sq ; Godbole:1989ti ;
Drees:1996rw ; Eboli:1997sv ; Yuan:1997tr ; Calucci:1997uw ; DelFabbro:1999tf
; Kulesza:1999zh ; Cattaruzza:2005nu ; Hussein:2006xr ; Maina:2009sj ;
Domdey:2009bg ; d'Enterria:2009hd ; Akesson:1986iv ; Abe:1997xk ; D0:2009 ,
and they may also feed pertinent standard model backgrounds to new physics
processes Sullivan:2008ki . They could be an issue in studies of Higgs boson
production in weak-boson-fusion since the “forward” jets could come from a
second hard subprocess.

##  V Acknowledgments

We benefited greatly from discussions with Dr. Thomas LeCompte and from
communications with Dr. John Campbell during the early development of this
project. We also thank Tom, John, and Professor Jianwei Qiu for valuable
comments and suggestions on an earlier draft of this paper. Research in the
High Energy Physics Division at Argonne is supported by the U. S. Department
of Energy under Contract No. DE-AC02-06CH11357. The research of GS at
Northwestern is supported by the U. S. Department of Energy under Contract No.
DE-FG02-91ER40684.

## References

  * (1) C. Goebel, F. Halzen and D. M. Scott, Phys. Rev. D 22, 2789 (1980).
  * (2) N. Paver and D. Treleani, Nuovo Cim. A 70, 215 (1982); Phys. Lett. B 146, 252 (1984); Z. Phys. C 28, 187 (1985).
  * (3) B. Humpert, Phys. Lett. B 131, 461 (1983).
  * (4) M. Mekhfi, Phys. Rev. D 32, 2371 (1985); Phys. Rev. D 32, 2380 (1985).
  * (5) B. Humpert and R. Odorico, Phys. Lett. B 154, 211 (1985).
  * (6) L. Ametller, N. Paver and D. Treleani, Phys. Lett. B 169, 289 (1986).
  * (7) F. Halzen, P. Hoyer and W. J. Stirling, Phys. Lett. B 188, 375 (1987).
  * (8) M. L. Mangano, Z. Phys. C 42, 331 (1989).
  * (9) R. M. Godbole, S. Gupta and J. Lindfors, Z. Phys. C 47, 69 (1990).
  * (10) M. Drees and T. Han, Phys. Rev. Lett. 77, 4142 (1996) [arXiv:hep-ph/9605430].
  * (11) O. J. P. Eboli, F. Halzen and J. K. Mizukoshi, Phys. Rev. D 57, 1730 (1998) [arXiv:hep-ph/9710443].
  * (12) F. Yuan and K. T. Chao, J. Phys. G 24, 1105 (1998) [arXiv:hep-ph/9706293].
  * (13) G. Calucci and D. Treleani, Nucl. Phys. Proc. Suppl. 71, 392 (1999) [arXiv:hep-ph/9711225]; Phys. Rev. D 80, 054025 (2009) [arXiv:0907.4772 [hep-ph]]; Phys. Rev. D 79, 074013 (2009) [arXiv:0901.3089 [hep-ph]].
  * (14) A. Del Fabbro and D. Treleani, Phys. Rev. D 61, 077502 (2000) [arXiv:hep-ph/9911358]; Phys. Rev. D 66, 074012 (2002) [arXiv:hep-ph/0207311].
  * (15) A. Kulesza and W. J. Stirling, Phys. Lett. B 475, 168 (2000) [arXiv:hep-ph/9912232].
  * (16) V. L. Korotkikh and A. M. Snigirev, Phys. Lett. B 594, 171 (2004) [arXiv:hep-ph/0404155].
  * (17) E. Cattaruzza, A. Del Fabbro and D. Treleani, Phys. Rev. D 72, 034022 (2005) [arXiv:hep-ph/0507052].
  * (18) M. Y. Hussein, Nucl. Phys. Proc. Suppl. 174, 55 (2007) [arXiv:hep-ph/0610207]; Nucl. Phys. Proc. Suppl. 152, 296 (2006); arXiv:0710.0203 [hep-ph].
  * (19) E. Maina, JHEP 0909, 081 (2009) [arXiv:0909.1586 [hep-ph]]; JHEP 0904, 098 (2009) [arXiv:0904.2682 [hep-ph]].
  * (20) S. Domdey, H. J. Pirner and U. A. Wiedemann, arXiv:0906.4335 [hep-ph].
  * (21) D. d’Enterria, G. K. Eyyubova, V. L. Korotkikh, I. P. Lokhtin, S. V. Petrushanko, L. I. Sarycheva and A. M. Snigirev, arXiv:0910.3029 [hep-ph].
  * (22) J. R. Gaunt and W. J. Stirling, arXiv:0910.4347 [hep-ph].
  * (23) T. Akesson et al. [Axial Field Spectrometer Collaboration], Z. Phys. C 34, 163 (1987).
  * (24) F. Abe et al. [CDF Collaboration], Phys. Rev. D 56, 3811 (1997); Phys. Rev. Lett. 79, 584 (1997).
  * (25) D0 Collaboration, “Double parton interactions in $\gamma$ \+ 3 jet events in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV in D0”, D0 Conference Note 5910-CONF (April 24, 2009); http://www-d0.fnal.gov/Run2Physics/WWW/results/prelim/QCD/Q13/Q13.pdf.
  * (26) A. H. Mueller and J. W. Qiu, Nucl. Phys. B 268, 427 (1986).
  * (27) X. F. Guo, Phys. Rev. D 58, 036001 (1998) [arXiv:hep-ph/9711453]; Nucl. Phys. A 638, 539C (1998) [arXiv:hep-ph/9801450].
  * (28) J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky and W. K. Tung, JHEP 0207, 012 (2002) [arXiv:hep-ph/0201195].
  * (29) F. Maltoni and T. Stelzer, JHEP 0302, 027 (2003) [arXiv:hep-ph/0208156].
  * (30) M. L. Mangano, M. Moretti, F. Piccinini, R. Pittau and A. D. Polosa, JHEP 0307, 001 (2003) [arXiv:hep-ph/0206293].
  * (31) H. Baer, V. Barger, G. Shaughnessy, H. Summy and L. T. Wang, Phys. Rev. D 75, 095010 (2007) [arXiv:hep-ph/0703289].
  * (32) Many next-to-leading order results can be obtained from the Monte Carlo program MCFM, J. M. Campbell and R. K. Ellis, Phys. Rev. D 60, 113006 (1999). Other NLO publications include L. J. Dixon, Z. Kunszt, and A. Signer, Phys. Rev. D 60, 114037 (1999) and Nucl. Phys. B 531, 3 (1998).
  * (33) Z. Sullivan and E. L. Berger, Phys. Rev. D 78, 034030 (2008) [arXiv:0805.3720 [hep-ph]].


# Understanding the Mechanisms of Deep Transfer Learning for Medical Images

Hariharan Ravishankar, Prasad Sudhakar,
Rahul Venkataramani, Sheshadri Thiruvenkadam, Pavan Annangi,
Narayanan Babu and Vivek Vaidya

_GE Global Research, Bangalore, India_

###### Abstract

The ability to automatically learn task specific feature representations has
led to a huge success of deep learning methods. When large training data is
scarce, such as in medical imaging problems, transfer learning has been very
effective. In this paper, we systematically investigate the process of
transferring a Convolutional Neural Network, trained on ImageNet images to
perform image classification, to kidney detection problem in ultrasound
images. We study how the detection performance depends on the extent of
transfer. We show that a transferred and tuned CNN can outperform a state-of-
the-art feature engineered pipeline and a hybridization of these two
techniques achieves 20% higher performance. We also investigate how the
evolution of intermediate response images from our network. Finally, we
compare these responses to state-of-the-art image processing filters in order
to gain greater insight into how transfer learning is able to effectively
manage widely varying imaging regimes.

##  1 Introduction

Automated organ localization and segmentation from ultrasound images is a
challenging problem because of specular noise, low soft tissue contrast and
wide variability of data from patient to patient. In such difficult problem
settings, data driven machine learning methods, and especially deep learning
methods in recent times, have found quite a bit of success. Usually, a large
amount of labeled data is needed to train machine learning models and a
careful feature engineering is required for each problem. The question of how
much data is needed for satisfactory performance of these methods is still
unanswered, with some recent works in this direction [7]. However,  _transfer
learning_ has been successfully employed in data scarce situations, with model
knowledge being effectively transferred across (possibly unrelated)
tasks/domains. It is fascinating that a model, learnt for an unrelated problem
setting can actually solve a problem at hand with minimal retraining. In this
paper, we have attempted to demonstrate and understand the effectiveness and
mechanism of transfer learning a CNN, originally learnt on camera images for
image recognition, to solve the problem of automated kidney localization from
ultrasound B-mode images.

Kidney detection is challenging due to wide variability in its shape, size and
orientation. Depending upon the acquisition scan plane, inconsistency in
appearance of internal regions (renal sinus) and presence of adjacent
structures like diaphragm, liver boundaries, etc. pose additional challenges.
This is also a clinically relevant problem as kidney morphology measurements
are essential in assessing renal abnormalities [13], planning and monitoring
radiation therapy, and renal transplant.

There have been semi-automated and automated kidney detection approaches
reported in literature. In [22], a texture model is built by an expectation
maximization algorithm using features inferred from a bank of Gabor filters,
followed by iterative segmentation to combine texture measures into parametric
shape model. In [16], Markov random fields and active contour methods have
been used to detect kidney boundaries in 3D ultrasound images. Recently,
machine learning approaches [1, 17] based on kidney texture analysis have
proven successful for segmentation of kidney regions from 2-D and 3-D
ultrasound images.

##  2 State of the art

CNNs [15] provide effective models for vision learning tasks by incorporating
spatial context and weight sharing between pixels. A typical deep CNN for a
learning task has, as input, $N$ channel image patches $P_{k}$ of size
$n_{1}\times n_{2}$, where
$P_{k}:\\{1,2,\cdots,n_{1}\\}\times\\{1,2,\cdots,n_{2}\\}\rightarrow
D\subset\mathbb{R}$, $k=1,2,\cdots,N$. The output is $M$ feature maps,
$G_{j}\in\mathbb{R}^{m_{1}\times m_{2}},j=1,2,\cdots,M$, defined as
convolutions using $MN$ filters $v^{j}_{k},(j=1,2,\cdots,M)$,
$(k=1,2,\cdots,N)$ of size $S=s_{1}\times s_{2}$, and $M$ scalars
$b^{j},j=1,2,\cdots,M$. We then have:

$G_{j}=\mathcal{S}_{\downarrow}\left(\sigma\left(\sum_{k}P_{k}*v^{j}_{k}+b^{j}\right)\right),j=1,2,\cdots,M.$
(1)

Here, $*$ denotes convolution, $\sigma$ is a non-linear function (sigmoid or a
linear cutoff (ReLU)). $\mathcal{S}_{\downarrow}$ is a down sampling operator.
The number of feature maps, filter size, and size of the feature maps are
hyperparameters in the above expression, with a total of $M(NS+1)$ parameters
that one has to optimize for a learning task. A deep CNN architecture is
multi-layered, with the above expression being hierarchically stitched
together, given the number of input/output maps, sizes of filters and maps for
each layer, resulting in a huge number of parameters to be optimized. When
data is scarce, the learning problem is under-determined and therefore
transferring CNN parameters from a pre-learned model helps.

For medical image problems, transfer learning is additionally attractive due
to the heterogeneity of data types (modalities, anatomies, etc.) and clinical
challenges. In [4], the authors perform breast image classification using a
CNN model trained on ImageNet. Shie et al. [18] employ the CaffeNet, trained
on ImageNet, to extract features and classify Otitis Media images. In [6], a
pre-trained CNN is used to extract features on ultrasound images to localize a
certain standard plane that is important for diagnosis.

Studies on transferability of features across CNNs include [23] and more
specifically [19],[21] for medical images. While our work demonstrates yet
another success of transfer learning for medical imaging and the tuning
aspects of transfer learning, we

  1. 1.

Reason out the effectiveness of transfer learning by methodically comparing
the response maps from various layers of transfer learnt network with
traditional image processing filters.

  2. 2.

Investigate the effect of level of tuning on performance. We demonstrate that
full network adaptation leads to learning problem specific features and also
establishes the superiority over off-the-shelf image processing filters.

  3. 3.

Re-establish the relevance and complementary advantages of state-of-the-art,
hand-crafted features and merits of hybridisation approaches with CNNs, to
help us achieve next level performance improvement [24].

##  3 Methods

From a set of training images, we build classifiers to differentiate between
kidney and non-kidney regions. On a test image, the maximum likelihood
detection problem of finding the best kidney region of interest (ROI) $S^{*}$
from a set of candidate ROIs $\\{S\\}$ is split into two steps, similar to [3,
17]. The entire set $\\{S\\}$ is passed through our classifier models and the
candidates with positive class labels ($Y$) are retained (Eq. (2)). The ROI
with highest likelihood ($L$) from the set $\\{S^{+}\\}$ is selected as the
detected kidney region (Eq. (3))

${\\{Y,L\\}=MLClassifier(S)\ \mathrm{and}\ \\{S^{+}\in S\mid Y=1\\}},$ (2)
$S^{*}=\arg\max(L^{+}),\ \mathrm{where}\ \\{L^{+}=L(S^{+})\\}.$ (3)

We propose to employ CNNs as feature extractors similar to [10] to facilitate
comparisons with traditional texture features. We also propose to use a well-
known machine learning classifier, to evaluate performance of different
feature sets, thereby eliminating the effects of having soft-max layer for
CNNs and a different classifier on traditional features as our likelihood
functions.

###  3.1 Dataset and Training

We considered a total of 90 long axis kidney images acquired on GE Healthcare
LOGIQ E9 scanner, split into two equal and distinct sets, for training and
validation. The images contained kidney of different sizes with lengths
varying between 7.5cm and 14cm and widths varying between 3.5cm and 7cm,
demonstrating wide variability in the dataset. The orientation of the kidneys
varied between -25∘ and +15∘. The images were acquired at varying depths of
ultrasound acquisition ranging between 9cm and 16cm. Accurate rectangular
ground truth kidney ROIs were manually marked by a clinical expert.

To build our binary classification models from training images, we swept the
field of view (FOV) to generate many overlapping patches of varying sizes (see
Fig. 1) that satisfied clinical guidelines on average adult kidney dimensions
and aspect ratio [8]. We downsampled these ROIs to a common size and were
further binned into two classes based on their overlap with ground truth
annotations. We used Dice similarity coefficient (DSC) as the metric and a
threshold of 0.8 (based on visual and clinical feedback) was used to generate
positive and negative class samples. This was followed by feature extraction
and model building.

(a)

(b)

(c)

Figure 1: (a) Sample Generation (b) Positive patches (c) Negative patches

###  3.2 Transfer Learned Features

Our study on transfer learning was based on adapting the popular CaffeNet [11]
architecture built on ImageNet database to ultrasound kidney detection, whose
simplified schematic is in Fig. 2. We extracted features after the ‘fc7’ layer
from all the updated nets, resulting in 4096 features. The features extracted
were:

  1. 1.

Full Network adaptation (CaffeNet_FA) \- Initialized with weights from
CaffeNet parameters, the entire network weights were updated by training on
kidney image samples from Section 3.1. The experiment settings were:
stochastic gradient descent update with a batch size of 100, momentum of 0.5
and weight decay of $5\times 10^{-4}$.

  2. 2.

Partial Network adaptation (CaffeNet_PA) \- To understand the performance
difference based on level of tuning, we froze the weights of ‘conv1’ and
‘conv2’ layers, while updating the weights of other layers. The reasoning
behind freezing the first two layers was to evaluate how sharable were the
low-level features and also to help us in interpret-ability (Sec. 5). The
experiment settings were same as those for full network adaptation.

  3. 3.

Zero Network adaptation (CaffeNet_NA) \- Finally, we also extracted features
from the original CaffeNet model without modifying the weights.

###  3.3 Traditional Texture Features

Some of the well-studied texture features used for ultrasound images include
(i) Haar features [3] for fetal anatomy studies, (ii) Gray Level Co-Occurrence
Matrix (GLCM) [20], (iii) Histogram of oriented gradient (HoG) (for automatic
view classification of echocardiogram images [12]).

Haar features have been reported to have the best performance for kidney
detection n [17]. For our study, we extracted Haar features similar to [3],
yielding a total of $\sim$2000 features.

###  3.4 Gradient Boosting Machine (GBM)

Ensemble classifiers have been shown to be successful in ultrasound organ
detection problems. In [3], authors have used probabilistic boosting tree
classifier for fetal anatomy detection. In [17], it has been noted that
gradient boosting machine (GBM) have outperformed adaboost classfiers. In an
empirical comparison study of supervised learning algorithms [5] comparing
random forests and boosted decision trees, calibrated boosted trees had the
best overall performance with random forests being close second. Motivated by
these successes, we have chosen to use Gradient boosting tree as our
classifier model. We build GBM classifiers for all the feature sets explained
in Section 3.2 and 3.3 using GBM implementation inspired by [2], with
parameters: shrinkage factor and sampling factor set to 0.5, maximum tree
depth = 2 and number of iterations = 200.

###  3.5 Hybrid approach

Investigation of the failure modes of baseline method (Haar + GBM) and
CaffeNet_FA revealed that they had failed on different images (Section 4). To
exploit the complementary advantages, we propose a simple scheme of averaging
the spatial likelihood maps from GBMs of these two approaches and employing it
in (2), which yields dramatic improvement.

Figure 2: Simple schematic of CaffeNet architecture [40]

##  4 Results

To quantitatively evaluate the performance on 45 validation images, we used
two metrics: (1) Number of localization failures - the number of images for
which the dice similarity coefficient between detected kidney ROI and ground
truth annotation was $<$ 0.80. (2) Detection accuracy - average dice overlap
across 45 images between detection results and ground truth, which . From
Table 1, we see that CaffeNet features without any adaptation outperformed
baseline by 2% in average detection accuracy with same number of failures.
This improvement is consistent with other results reported in literature [18],
where CaffeNet features outperform state-of-the-art pipeline. However, by
allowing these network weights to get adapted to the kidney data, we achieved
a performance boost of 4% over the baseline method, with number of failure
cases reducing to 10 from 12. Interestingly, tuning with the first two
convolutional layers frozen yielded intermediate performance, suggesting that
multiple levels of feature adaptation are important to the problem.

Fig. 3LABEL:sub@subfig:a and LABEL:sub@subfig:b shows a case in which the
baseline method was affected by the presence of diaphragm, kidney and liver
boundaries creating a texture similar to renal-sinus portion, while CaffeNet
had excellent localization. Fig. 3LABEL:sub@subfig:c and LABEL:sub@subfig:d
illustrate a case where CaffeNet resulted in over-segmentation containing the
diaphragm, clearly illustrating that in limited data problems careful feature-
engineering incorporating domain knowledge still carries a lot of relevance.
Finally, we achieved a best performance of 86% average detection accuracy
using the hybrid approach (Section 3.5). More importantly, the number of
failures of the hybrid approach was 3/45, which is 20% better than either of
the methods.

Method |  Haar Features |  CaffeNet_NA |  CaffeNet_PA |  CaffeNet_FA |  Haar + CaffeNet_FA
---|---|---|---|---|---
Average Dice overlap |  0.793 |  0.825 |  0.831 |  0.842 |  0.857
# of failures |  12/45 |  12/45 |  11/45 |  10/45 |  3/45
Table 1: Performance Comparison on unseen 45 validation images

(a) Haar

(b) CaffeNet_FA

(c) Haar

(d) CaffeNet_FA

Figure 3: Visual comparison of Baseline Method with CaffeNet Transfer

##  5 Discussion

(a) Sample Kidney patch

(b) Phase Congruency

(c) Frangi Vesselness

(d) CaffeNet_NA L1_1

(e) CaffeNet_NA L1_2

(f) CaffeNet_NA L2_1

(g) CaffeNet_FA L1_1

(h) CaffeNet_FA L1_2

(i) CaffeNet_FA L2_1

Figure 4: Filter responses for learned & hand-engineered features for a sample
patch (CaffeNet_FA(NA) LX_Y denotes response image from $Y^{th}$ filter of
layer $X$)

It is indeed very interesting to see that features learnt on camera images
were able to outperform careful feature engineering on sharply different
detection problems, in modalities whose acquisition physics are distinctly
different. Fig. 4 compares some of the response images generated from layers 1
and 2 of the learned network with traditional image processing outputs like
Phase Congruency [14] and Frangi vesselness filter [9] for an example patch.

Here, we would like to highlight two main points: (1) Visually, we find the
output has intriguing similarities with the outputs of hand crafted feature
extractors optimized for Ultrasound. The response maps of Fig.
4LABEL:sub@subfig:CNFA_L11 and LABEL:sub@subfig:CNFA_L2 are similar to
4LABEL:sub@subfig:PC and LABEL:sub@subfig:frangi. This is very encouraging
because of the fact that CNNs learns features that are equivalent to some of
these widely used non-linear feature extractors. (2) The second important
observation here is the reduction in speckle noise on CaffeNet_FAL1_1,
compared to CaffeNet_PAL1_1. By carefully tuning CaffetNet features on
ultrasound data, the model was able to learn the underlying noise
characteristics, while preserving edges, and this resulted in a much improved
response map as shown in Fig. 4LABEL:sub@subfig:CNFA_L11 and
LABEL:sub@subfig:CNFA_L2.

|  |  |  Layer |  |
---|---|---|---|---|---
|  Conv1 |  Conv2 |  Conv3 |  Conv4 |  Conv5
# of filters with $\geq$ 40% change |  0 |  5 |  125 |  22 |  62
Table 2: No. of filters in each layer that changed by more than 40% in
$\ell_{2}$ norm

Further, we quantitatively analyzed changes (% in $\ell_{2}$ norm) in filter
weights in each layer to identify significant trends. Table 2 shows a large
number of filters have significantly changed in the 3rd layer, with filters in
the 1st and 2nd layer showing minimal change. This is possibly due to the
lower level features being fairly the same for both natural and ultrasound
images. We also noted that the use of ReLU as the activation function also
avoided the vanishing gradient problem, resulting in this skew in distribution
of weight changes across layers. The response images past layer 2 proved to be
difficult to interpret, and may require more intensive techniques. Our
quantitative results and the literature in the field show that a great deal of
the power of deep networks lies in these layers, and so we feel this is an
important area for our future investigation.

In a clinical context, the interpretability of models is crucial and we feel
this insight into why the deep CNN was able to outperform hand-crafted
features is as important as the results demonstrated in Sec. 4. We also see
this as opening up new ways of understanding and utilizing deep networks for
medical problems.

## References

  * [1] Roberto Ardon, Remi Cuingnet, Ketan Bacchuwar, and Vincent Auvray.  Fast kidney detection and segmentation with learned kernel convolution and model deformation in 3d ultrasound images.  In Proc. of ISBI, pages 268–271, 2015.
  * [2] Carlos Becker, Roberto Rigamonti, Vincent Lepetit, and Pascal Fua.  Supervised feature learning for curvilinear structure segmentation.  In Proc. of MICCAI 2013, pages 526–533, 2013.
  * [3] G. Carneiro, B. Georgescu, S. Good, and D. Comaniciu.  Detection and measurement of fetal anatomies from ultrasound images using a constrained probabilistic boosting tree.  IEEE Trans. on Med. Imag., 27(9):1342–1355, Sept 2008.
  * [4] Gustavo Carneiro, Jacinto Nascimento, and Andrew P. Bradley.  Unregistered multiview mammogram analysis with pre-trained deep learning models.  In Proc. of MICCAI, pages 652–660, 2015.
  * [5] Rich Caruana and Alexandru Niculescu-Mizil.  An empirical comparison of supervised learning algorithms.  In Proc. of ICML, pages 161–168, 2006.
  * [6] H. Chen, D. Ni, J. Qin, S. Li, X. Yang, T. Wang, and P. A. Heng.  Standard plane localization in fetal ultrasound via domain transferred deep neural networks.  IEEE Journal of Biomedical and Health Informatics, 19(5):1627–1636, Sept 2015.
  * [7] Junghwan Cho, Kyewook Lee, Ellie Shin, Garry Choy, and Synho Do.  How much data is needed to train a medical image deep learning system to achieve necessary high accuracy?  In ICLR, 2016.
  * [8] S. A. Emamian, M. B. Nielsen, J. F. Pedersen, and L. Ytte.  Kidney dimensions at sonography: correlation with age, sex, and habitus in 665 adult volunteers.  American Journal of Roentgenology, 160(1):83–86, Jan 1993.
  * [9] A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever.  Multiscale vessel enhancement filtering.  In Proc. of MICCAI, pages 130–137, 1998.
  * [10] Ross Girshick, Jeff Donahue, Trevor Darrell, and Jitendra Malik.  Rich feature hierarchies for accurate object detection and semantic segmentation.  2014\.
  * [11] Yangqing Jia, Evan Shelhamer, Jeff Donahue, Sergey Karayev, Jonathan Long, Ross Girshick, Sergio Guadarrama, and Trevor Darrell.  Caffe: Convolutional architecture for fast feature embedding.  In Proc. of the ACM Int. Conf. on Multimedia, pages 675–678, 2014\.
  * [12] Eystratios G. Keramidas, Dimitris K. Iakovidis, Dimitris Maroulis, and Stavros Karkanis.  Efficient and effective ultrasound image analysis scheme for thyroid nodule detection.  In Proc. of ICIAR, pages 1052–1060, 2007.
  * [13] Arpana M. Kop and Ravindra Hegadi.  Kidney segmentation from ultrasound images using gradient vector force.  IJCA, Special Issue on RTIPPR, 2:104–109, 2010.
  * [14] Peter Kovesi.  Phase congruency detects corners and edges.  In Proc. of The Australian Pattern Recognition Society Conference: DICTA, pages 309–318, 2003.
  * [15] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton.  Imagenet classification with deep convolutional neural networks.  In Proc. of Adv. in NIPS, pages 1106–1114, 2012.
  * [16] Marcos Martin-Fernandez and Carlos Alberola-Lopez.  An approach for contour detection of human kidneys from ultrasound images using markov random fields and active contours.  Medical Image Analysis, 9(1):1 – 23, 2005.
  * [17] H Ravishankar and Pavan Annangi.  Automated kidney morphology measurements from ultrasound images using texture and edge analysis.  In SPIE Med. Imag., 2016.
  * [18] Chuen-Kai Shie, Chung-Hisang Chuang, Chun-Nan Chou, Meng-Hsi Wu, and E. Y. Chang.  Transfer representation learning for medical image analysis.  In EMBC, pages 711–714, 2015.
  * [19] Hoo-Chang Shin, Holger R. Roth, Mingchen Gao, Le Lu, Ziyue Xu, Isabella Nogues, Jianhua Yao, Daniel J. Mollura, and Ronald M. Summers.  Deep convolutional neural networks for computer-aided detection: CNN architectures, dataset characteristics and transfer learning.  IEEE Trans. Med. Imaging, 35(5):1285–1298, 2016.
  * [20] A. S. M. Sohail, M. M. Rahman, P. Bhattacharya, S. Krishnamurthy, and S. P. Mudur.  Retrieval and classification of ultrasound images of ovarian cysts combining texture features and histogram moments.  In Proc. of ISBI, pages 288–291, April 2010.
  * [21] Nima Tajbakhsh, Jae Y. Shin, Suryakanth R. Gurudu, R. Todd Hurst, Christopher B. Kendall, Michael B. Gotway, and Jianming Liang.  Convolutional neural networks for medical image analysis: Full training or fine tuning?  IEEE Trans. Med. Imaging, 35(5):1299–1312, 2016.
  * [22] Jun Xie, Yifeng Jiang, and Hung-Tat Tsui.  Segmentation of kidney from ultrasound images based on texture and shape priors.  IEEE Trans. on Med. Imag., 24(1):45–57, 2005.
  * [23] Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson.  How transferable are features in deep neural networks?  In Proc. of Adv. in NIPS, pages 3320–3328, 2014.
  * [24] Yefeng Zheng, David Liu, Bogdan Georgescu, Hien Nguyen, and Dorin Comaniciu.  3d deep learning for efficient and robust landmark detection in volumetric data.  In Proc. of MICCAI, pages 565–572, 2015.