# Landau Level Single-Electron Pumping

E. Pyurbeeva<EMAIL_ADDRESS>Fritz Haber Center for
Molecular Dynamics, Institute of Chemistry, The Hebrew University of
Jerusalem, Jerusalem 91904, Israel  School of Physical and Chemical Sciences,
Queen Mary University of London, Mile End Road, London E1 4NS, UK  M.D.
Blumenthal<EMAIL_ADDRESS>Department of Physics, University of
Cape Town, Rondebosch, Cape Town, 7700, South Africa  J.A. Mol  School of
Physical and Chemical Sciences, Queen Mary University of London, Mile End
Road, London E1 4NS, UK  H. Howe  University College London, London Centre for
Nanotechnology and Department of Electrical Engineering, London, WC1A 0AH, UK
H. E. Beere  University of Cambridge, Physics Department, Cavendish
Laboratory, Cambridge, CB3 0HE, UK  T. Mitchell  University of Cambridge,
Physics Department, Cavendish Laboratory, Cambridge, CB3 0HE, UK  D. A.
Ritchie  University of Cambridge, Physics Department, Cavendish Laboratory,
Cambridge, CB3 0HE, UK  M. Pepper  University College London, London Centre
for Nanotechnology and Department of Electrical Engineering, London, WC1A 0AH,
UK

###### Abstract

We present the first detailed study of the effect of a strong magnetic field
on single-electron pumping in a device utilising a finger-gate split-gate
configuration. In the quantum Hall regime, we demonstrate electron pumping
from Landau levels in the leads, where the measurements exhibit pronounced
oscillations in the lengths of the pumping plateaus with the magnetic field,
reminiscent of Shubnikov-de Haas oscillations. This similarity indicates that
the pumping process is dependent on the density of states of the 2D electron
gas over a narrow energy window. Based on these observations, we develop a new
theoretical description of the operation of single-electron pumps which for
the first time allows for the determination of the physical parameters of the
experiment; such as the capture energy of the electrons, the broadening of the
quantised Landau levels in the leads, and the quantum lifetime of the
electrons.

††preprint: APS/123-QED

##  I Introduction

Technological applications, such as quantum information processing,
nanoelectronics, and electron quantum optics require the control of individual
electrons with a high degree of accuracyPekola  _et al._ (2013); Keller  _et
al._ (1997); Bäuerle  _et al._ (2018); Kaneko, Nakamura, and Okazaki (2016);
Dupont-Ferrier  _et al._ (2013); DiVincenzo and IBM. (2000); Wright  _et al._
(2009); Maisi  _et al._ (2009); Giblin  _et al._ (2019); Yamamoto  _et al._
(2012); Ji  _et al._ (2003); Ubbelohde  _et al._ (2015); Gumbs  _et al._
(2009); Fletcher  _et al._ (2013); Blumenthal  _et al._ (2023); Ubbelohde  _et
al._ (2023); Brange and Flindt (2023). AlGaAs high frequency single-electron
pumpsBlumenthal  _et al._ (2007); Kaestner  _et al._ (2008) capable of
delivering quantised charges, have been extensively studied for this purpose
in the past decades, with considerable advancements madeKaestner and
Kashcheyevs (2015). In silicon based electronYamahata, Nishiguchi, and
Fujiwara (2014) and hole pumpsYamahata and Fujiwara (2024), there has been
growing interest in single-electron pumps that utilise localised states, such
as dopants implanted in a channel, instead of electrically defined
islandsMoraru  _et al._ (2007); Roche  _et al._ (2013); Lansbergen, Ono, and
Fujiwara (2012). Notably, the high-speed operation of a single donor pump has
been demonstratedTettamanzi, Wacquez, and Rogge (2014), highlighting the
potential of this approach.

In the characterisation of deterministic single-electron pumps, the Universal
Decay Cascade (UDC) model presented in the seminal work by _Kashcheyevs et al_
Kashcheyevs and Kaestner (2010) remains the only paradigm. It gives good
agreement with experimental data and allows for the determination of a set of
parameters (a “fingerprint”) that quantifies the individual device’s accuracy
and dynamics. The model, however, is limited in its ability to explain
experimentally observed physical phenomena. This limitation needs to be
overcome as more complex devices are utilised as electron pumpsRossi  _et al._
(2018).

We present experimental observation of a split-gate finger-gate (SFG) electron
pumpHowe  _et al._ (2021) operating in a perpendicular magnetic field. The
evolution of the pumped current is studied as the magnetic field is swept from
$1$ T to $9$ T. We find that the pumping plateaus do not lengthen
monotonically with field, as has been previously suggestedWright  _et al._
(2008); Kaestner  _et al._ (2009); Leicht  _et al._ (2010); Hanief  _et al._
(2024), but instead their evolution is closely correlated with the
Shubnikov–de Haas (SdH) oscillations in the same device. This agreement allows
us to introduce a new physical theoretical model, an extension of the UDC
model, separating the individual properties of the device from the effect of
the magnetic field. It gives new insight into the dynamics of the pumping
process, describes the data qualitatively and allows us to determine certain
numerical parameters of the system, such as the capture energy of the pump and
the electrons’ quantum lifetime.

##  II Experimental Data

Figure 1: a) A false colour SEM image of the SFG device with experimental
setup. b) 2D colour map showing the dependence of the pumped current on the
magnetic field and exit voltage $V_{\rm{exit}}$. c) A set of schematics
illustrating a single pump cycle in a SFG pump. (i) At low momentary values of
$V_{\rm{ent}}(t)$ the potential profile presents a single maximum, preventing
direct current from the source to the drain. (ii) As $V_{\rm{ent}}(t)$
increases, a confined state forms in the potential profile. (iii) Further
increase of $V_{\rm{ent}}(t)$ isolates the confined electron from both the
source and the drain. (iv) The confined state is destroyed and electron is
deposited into the drain. d) Pumping traces in $V_{\rm{exit}}$ for several
values of $B$. The lengthening of the plateaus can be observed initially, in
agreement with previous work Howe  _et al._ (2021), but it is not monotonic in
B. e) The starting voltage of the first plateau, $\delta_{1}$ (solid purple)
vs. the Shubnikov-de Haas effect in the same sample with RF excitation on the
finger-gate (dashed turquoise) and without (dashed orange).

The SEM image of the pump and the schematic of the experiment are shown in
Fig.1a, with further details given in the Methods section. The principle of
device operation is illustrated in Fig.1c.

The time evolution panels (i) - (iv) show a cross-section of the spatial 3D
potential created by the gates. At the beginning of the cycle (i), the left
entrance barrier is low and the dot potential profile has a single maximum. As
the entrance barrier potential increases (ii), a local minimum appears in the
potential, creating confined states, which can be occupied by electrons from
the source (capture stage). With further increase (iii), the confined states
become isolated from the source and move with the evolution of the potential
profile. Finally (iv), the confined states are destroyed and electrons are
ejected into the drain.

The number of pumped electrons per cycle can be fine-tuned at constant RF
frequency and amplitude by both the DC gate voltages. Fig.1d shows four
pumping traces at different magnetic fields for a fixed entrance gate offset
voltage and RF parameters. The pumped current evolves as a function of the
exit gate voltage $V_{\rm{exit}}$, with the integer number of pumped electrons
increasing as the DC voltage on the exit barrier is made more positive (the
barrier is lowered).

Fig.1b shows a pumpmap of the current as a function of $V_{\rm{exit}}$ and
magnetic field $B$Wright  _et al._ (2010) with the entrance gate DC voltage
fixed at $V_{\text{ent}}=-600$ mV, the source-drain bias $V_{\text{SD}}=100$
mV, and the RF amplitude and frequency $V_{\text{Amp}}=300$ mV and
$f_{\text{RF}}=180$ MHz. On observation, the plateau lengths depend strongly
on the magnetic field, and the general trend over the field $1-3$ T is towards
plateau lengthening, in agreement with previous studiesWright  _et al._
(2008); Kaestner  _et al._ (2009). The plateau length dependence is however
non-monotonic, but instead oscillatory, including a dramatic resonance peak at
around $7.2$ T.

These oscillations are reminiscent of magnetic field dependencies seen in the
Shubnikov-de Haas effect. We study the SdH effect in the pump by measuring the
longitudinal resistance of the device as a function of magnetic field, first
with the finger gate grounded (Fig.1e, orange) and then with an RF signal on
the finger gate (Fig.1e, turquoise) with the parameters of the RF signal set
as used during pumping. These SdH plots in the quantum Hall regime are
compared to the exit gate voltage corresponding to the start of the first
pumping plateau (current is at $1/e$ of the expected plateau value, see Eq.1)
as a function of magnetic field (Fig.1e, blue).

The SdH plot without the RF signal differs from assessment measurements in the
same wafer with a standard Hall bar. This difference can be attributed to the
geometry of the mesa in the pump deviating substantially from a standard Hall
bar geometry. The inclusion of Ti/Au gates as well as the fabrication process
itself also affect the wafer, altering the electron density and mobility. The
application of the RF signal further influences the magnetic field dependence
of the device resistance, however both retain the general peak structure of
the signature SdH.

As the magnetic field oscillations in the SdH effect originate in the Landau
quantisation of the electron states, we make the conjecture that similarly,
the pumping of electrons depends on the population of the density of states in
the source side of the pump. The rest of the work expands on this assumption.

##  III Theoretical model

Figure 2: a) The onset of pumping at 1/e of the true integer plateau of the
first four plateaus ($n=1...4$) linearly scaled to coincide with each other.
b) The calculated scaling factors of $\delta_{1}$–$\delta_{4}$. c) Normalised
pumped current vs. $V_{\rm{exit}}$ for varying magnetic fields from $T=1$ T to
$T=9$ T in $0.5$ T steps. Red dashed line shows the UDC fit.

###  III.1 Motivations and criteria for the model

A single-electron pump is a complex open quantum system with many aspects
which make it a serious challenge for theoretical descriptionPyurbeeva, Mol,
and Gehring (2022): both particle and energy exchange being present, a time-
dependent Hamiltonian, multiple competing time-scales, out-of-equilibrium
conditions and multi-electron processes. Combinations of these are studied in
topics such as driven quantum dots, for applications such quantum heat engines
or refrigeratorsJuergens  _et al._ (2013); Hino and Hayakawa (2021); Monsel
_et al._ (2022); Monsel, Schulenborg, and Splettstoesser (2023), or far-from-
equilibrium thermodynamics in nanodevicesPyurbeeva, Thomas, and Mol (2023).

However, single-electron pumps exhibit additional complexity due to
qualitative changes, such as the creation and destruction of confined states,
and strong system-bath coupling associated with it, alongside quantitative
changes like cycling of addition energy. Recent researchSchulenborg  _et al._
(2023) has focused on dynamical analysis from an open quantum system
perspective, including a study on two-electron emission processes from a
quantum dot, however such approaches are computationally involved. Our
investigation, focusing on the capture side of the cycle and incorporating a
magnetic field, takes a more qualitative and intuitive approach.

The main result of the UDC model is a double-exponential fitting equation for
the pumping curvesHowe  _et al._ (2021):

$I=ef\sum_{n}\exp\left[-\exp\left(-\alpha_{n}(V_{\rm{exit}}-\delta_{n})\right)\right]$
(1)

where subscript $n$ is the plateau number, $\alpha_{n}$ is the plateau’s lever
arm, and $\delta_{n}$ – its position in $V_{\rm{exit}}$. These parameters act
as a fingerprint of the pump device, allowing for the pump comparison and
optimisationKashcheyevs and Kaestner (2010).

Equation 1 gives a good fit to the experimental data in a wide variety of
electron pump geometriesHowe  _et al._ (2021), but despite later
expansionsKashcheyevs and Timoshenko (2012, 2014) on the model, it does not
link explicitly to any key physical parameters such a temperature, B field, RF
signal amplitude, and frequency.

Our experimental data (Fig.1b) provides the first systematic high-precision
study of the dependence of the pumping on an external parameter, the magnetic
field. We now construct a model to describe our findings.

In Fig.2a the starting voltages (corresponding to the current reaching $1/e$
of the true integer plateau) for $n=1$ (blue) to $n=4$ (red) are given as a
function of $B$. These line scans were taken from Fig.1b and linearly scaled
by a factor of ($\eta$), shown in Fig.2b. The scaling factors follow a power
law in $n$ with the power $a=1.58$.

The fact that $\delta_{1}$–$\delta_{4}$ coincide up to a linear transformation
(Fig.2a), indicates that there is a single parameter $\delta(B)$ governing the
overall length of the pumping trace, rather than the start of each plateau
having separate magnetic field dependencies. Fig.2c shows the pumping traces
at different magnetic fields stretched by the length factor determined in
Fig.2a, further supporting this claim, and showing a “true” pumping trace
inherent to the device.

The length function $\delta(B)$, bears a strong resemblance to the SdH
oscillations of longitudinal resistance of the 2DEG which is proportional to
the density of states at the Fermi-level as discussed previously (Fig.1e). We
can conclude that the pumping process not only depends on the density of
states in the source electrode, but on the density of states at a single
energy – if $\delta(B)$ was governed by electron exchange processes at
different energies, it would be a convolution of $n(B,E)$ at different values
of $E$, and the contributions from different Landau levels would be washed
out. The sharp Landau level signatures observed indicate that the entire
pumping effect is determined by exchange processes confined to a very narrow
energy window.

###  III.2 The 0-DIP model and dot parameters

We now propose a physical model based on our observations. We utilise a 1D
system and describe the electron energy profiles created by the application of
the voltages $V_{\rm{ent}}$ and $V_{\rm{exit}}$ to the entrance and exit
electrodes as smooth differentiable peak functions $\varphi_{\rm{ent}}$ and
$\varphi_{\rm{exit}}$ with continuous first and second derivatives. The
maximal values of the electron energies are $E_{\rm{ent}}$ and
$E_{\rm{exit}}$, which depend linearly on the voltages applied (Fig.3a).

If the peaks are sufficiently close to one another, for every value of
$E_{\rm{exit}}$ there exist two values of $E_{\rm{ent}}$ for which the total
potential profile exhibits an inflection point with a zero derivative, with
the inflection points on the source side for low $E_{\rm{ent}}$ and drain side
for high $E_{\rm{ent}}$. We call the former configuration the _0-DIP
potential_ (Fig.3b). The value of $E_{\rm{ent}}$ at the 0-DIP potential is
proportional to $E_{\rm{exit}}$, where we define the proportionality
coefficient as $\alpha$.

The significance of the 0-DIP potential is that it approximately separates two
configurations of the dynamically defined dot – that with no confined electron
states ($E_{\rm{ent}}<\alpha E_{\rm{exit}}$), and that other containing one or
multiple confined states. We assume that confined states form soon after
$E_{\rm{ent}}$ passes the 0-DIP potential. A crucial point in our model is
that a further increase of $E_{\rm{ent}}$ with the application of the RF
signal leads to a rapid reduction in the tunnel coupling between the now-
formed dot and the source (Fig.3c). After a small increment $dE$ in entry
barrier height, the dot is completely isolated, not only from the drain, but
also from the source (grey shaded profile in Fig.3b).

The value of $\differential E$ is assumed small in comparison to the
characteristic broadened width $\Gamma$ of the Landau levels. This assumption
contributes to the explanation of the oscillatory behaviour of the pumped
current vs magnetic field in Fig.1b and is a requirement of our model.

To demonstrate the validity of this 0-DIP model, we investigate the evolution
of the geometry of the quantum dot during the pump cycle. We define two
fundamental dimensions of the dot: the width of the entrance potential barrier
$W_{d}$ and the depth of the dot $D_{d}$ (Fig.3d) and model the change in
these dimensions as a function of $\differential E$ in a one-dimensional
system. Analytically making use of Taylor expansion (see Methods) we find the
following two expressions for the dot parameters:

$\begin{dcases}W_{d}=C_{W}\sqrt{\frac{\differential
E}{E_{\rm{exit}}}}\\\\[5.0pt] D_{d}=C_{D}\sqrt{\frac{(\differential
E)^{3}}{E_{\rm{exit}}}}\end{dcases}$ (2)

While the exact form of the dependence will be different for a two-dimensional
case and will depend on the geometry of the electrodes, qualitatively it will
follow the 1D case, with the barrier height and width increasing with
$\differential E$ and decreasing with $E_{\rm{exit}}$. The effect of dot size
on pumping can be seen in the well-known universal feature of the pumping
traces – less negative $V_{\rm{exit}}$ (lower $E_{\rm{exit}}$) leading to more
electrons being pumped per cycle, which intuitively agrees with a
geometrically larger quantum dot. Additionally, the increase of $W_{d}$ with
$\differential E$ supports the initial assumption of the dot only being
coupled to the source in a narrow energy window, as the coupling to the source
$\Gamma_{S}$ depends exponentially on $W_{d}$ (Fig.3c).

##  IV Data analysis – density of states and capture energy

Figure 3: a) A labelled schematic of the theoretical model. b) The device
potential profile in the 0-DIP configuration. The grey shaded profile shows
the quantum dot at the cut-off point where it becomes isolated from the source
electrode. c) The schematic showing the change in tunnel rate $\Gamma$ and
source electron density $n(E)$ as the energy of the entrance barrier is
changed ($\differential E$). d) Schematic showing the zoomed-in region of the
minima in the QD potential with the length scale dimensions ($W_{d}$ width,
and $D_{d}$ depth) indicated. The grey area indicates the 0-DIP configuration.
e) A schematic of the effect of $E_{\rm{exit}}$ on the pumping. Lowering the
energy of the exit barrier leads to an increase in the size of the dot at the
0-DIP potential configuration, which lowers energy of capture on
$V_{\rm{ent}}$ given by $\alpha\Delta E_{\rm{exit}}$. In the case of a very
narrow band in $n(E)$ (shown in dark grey) the change in capture energy could
lead to a reduction of the pumping current despite an increase in the size of
the dot.

We now apply the developed 0-DIP model to the experimental data. This allows
us to determine certain numerical properties of the pumping mechanism in high
fields and validate our model.

Figure 4: a) Schematics of the evolution of the spin-split first Landau level
for $B$ values chosen around the main resonance peak centred at $B=7.3$T. The
black shape markers indicate the chosen $B$ values on the resonance peak in
panel b. The Fermi energy is indicated by the red vertical lines, the
turquoise dotted line indicates a fixed capture energy and the green line a
changing capture energy, with $B$. b) Onset of 1st plateau pumping
(experimental data) in red taken directly from Fig.1e. Blue and turquoise
plots indicate the simulated density of states for changing capture energy and
fixed capture energy respectively vs $B$ field. $\Gamma$ shows from where in
the experimental data (red) plot the Landau broadening was determined.

The electron pump, defined in a quasi one-dimensional restriction etched into
a 2DEG AlGaAs mesa with multiple Ohmic contacts, allowed for routine Quantum
Hall and SdH measurements to determine the electron density of the pump mesa.
Initially, all gates were grounded, and then measurements were repeated with
the finger gate oscillating at the same frequency and amplitude used during
the pump’s operation. These measurements revealed an electron density of
$n_{D}=1.53\times 10^{15}$ m-2, which corresponds favourably with what was
measured by the growers ($n_{D}=1.87\times 10^{15}$ m-2) in the Cavendish
Laboratory on the stock wafer before any fabrication on the wafer.

From the SdH oscillations and measured electron density we were in a position
to determine the filling factors $\nu$. We observe that the final minimum in
the SdH data (Fig.1e orange), where the Fermi level is in a localised state
corresponds with all the electrons in the first Landau level ($\nu=2$,
including spin). This aligns with the minimum in the pumping data at $6.7$,T,
occurring just before the onset of the significant resonance peak at $7.3$,T
(Fig.1e). The last minimum in the SdH data shows a stronger coincidence with
the minima in the pumped data if the RF gate when measuring the SdH was
switched on.
Using the above, we then simulated using code written in Matlab, the movement
of the Landau levels, and the Fermi energy governed by the determined density
of electrons above as the magnetic field was changed. We focused initially on
the region of the resonance peak at B values of $6.0$ T, $6.6$ T, $7.3$ T,
$8.2$ T and $9.0$ T.

First we determined the Landau level broadening (fwhm), a requirement for
observing the quantum Hall effectPotts  _et al._ (1996), which is given by
$\Gamma=\frac{\hbar}{\tau_{i}}$, where $\tau_{i}$ is the quantum lifetime.
$\Gamma$ is determined from the resonance peak (red line, Fig.4b), which as
previously discussed in the 0-DIP model, corresponds to the dot probing the
Landau level with filling factor $\leq 2$ (including spin) as the capture
energy changes monotonically as B is swept (Fig.4a). We obtained
$\Gamma=1.35\times 10^{-22}$,J, corresponding to a lifetime of
$\tau_{i}=0.78$,ps. This lifetime, the mean time a carrier remains in a
momentum eigenstate before scatteringQian  _et al._ (2017), is an order of
magnitude shorter than typical time scales from standard quantum Hall
measurementsMartin  _et al._ (1988). In the pumping regime, introducing the RF
signal disrupts the system, increasing scattering and shortening the quantum
lifetime. Further investigations will explore the RF frequency and amplitude’s
impact on $\tau_{i}$.

The movement of the Landau levels (black peaks) and Fermi level (red lines)
are shown in the selection of plots in Fig.4a with the width of the Landau
peaks as determined above. At the resonance peak at $7.3$ T the Fermi level
despite not being centred in the lowest Landau level is contained within it.
With this finding we then determined the capture energy of the electrons. We
start with the assumption that at the resonance peak at $7.3$ T, the capture
energy given by the green line is at the centre of the first Landau level for
$b=7.3$ T(Fig.4a $\blacksquare$). It is also below the Fermi level given in
red, which is essential as there are no electrons for pump capture at energies
above the source lead Fermi level.

We then modelled two different scenarios. First, we assumed that the capture
energy of the electrons during the pump cycle does not change as a function of
magnetic field and we pinned the capture energy at the centre of the lowest
Landau level $E_{c}=1.51\times 10^{-21}$ J for all B as shown in Fig. 4a by
the dotted turquoise line. We then monitored the change in the density of
states at this fixed capture energy as a function of B which is then plotted
as the dotted turquoise line in Fig.4b. Second, we assumed that the capture
energy is changing as a function of magnetic field linearly. We scaled the
capture energy at $B=7.2$ T ($E_{c}=1.51\times 10^{-21}$ J) over the full B
field from $1$ T to $9$ T. The capture energy therefore changed by $\Delta
E_{C}=3.9\times 10^{-22}$ J over the full range of B field. The plot of the
density of states as the capture energy was changed from $E_{c}=1.21\times
10^{-21}$ J to $E_{c}=1.60\times 10^{-21}$ J in step with the change in B from
$1$ T to $9$ T, in multiples of $\Delta B=9.0\times 10^{-4}$ T is shown by the
dark blue solid trace in Fig.4b).

From these two plots we see a much stronger coincidence between the dark blue
and red onset pumping curve allowing us to conclude the capture energy of the
electrons changes, albeit slightly, over the full magnetic field sweep.

##  V Conclusion

Our experimental data of single-electron pumping through a split-gate finger-
gate quantum dot at high magnetic fields demonstrates oscillatory behaviour
reminiscent of the SdH effect. The inclusion of a split-gate on the exit side
of the quantum dot allows for an increase in gate resolution due to a lower
lever-arm term. This, with the change in potential profile on the exit side of
the quantum dot from a simple Gaussian to a saddle-point potential, has
allowed for further probing of the dynamics of these on-demand electron
sources. Our analysis together with a new pumping model based on the
foundations of the UDC model has allowed for the first reported direct
measurement of the capture energy, the broadening of Landau levels in such
pump devices and therefore the quantum lifetime of the pumped electrons. This
demonstration of single electron pumping in the quantum Hall regime, where
electrons are captured from specific Landau levels, shows the capture energy
to be independent of the exit gate DC voltage for fixed B field, but a change
in the capture energy is observed as the B field itself is varied.

Future work will explore further the frequency and temperature dependence of
pumps operating in this regime as well the ejection energy of the pumped
electrons.

## Acknowledgements

This work was supported by the UK Engineering Physical Sciences Research
Council, grants EP/K040359/1 and EP/R029075/1, the National Research
Foundation of South Africa. J.A.M. acknowledges a UKRI Future Leaders
Fellowship, grant no. MR/X023125/1.

## Author contributions

E.P developed the theoretical model, produced data plots, interpreted the data
and contributed to the drafting of the paper. M.D.B ran the experiment,
carried out the simulation work, interpreted the data and contributed to the
drafting of the paper. Both E.P and M.D.B contributed equally. H.H. designed
and processed the devices. J.A.M contributed to data interpretation. T.M
carried out the e-beam lithography. H.E.B and D.A.R developed and grew the
GaAs/AlGaAs heterostructures used. M.P. provided support, references, and
useful discussions. All authors reviewed the manuscript.

## Methods

###  V.1 Fabrication and Experiment Setup

Our device comprises a single-electron pump with two Ti/Au gates in a finger-
gate split-gate geometry (red and green in Fig.1a respectively). Only the
entrance gate, in our the case the finger-gate, is coupled directly to the
Agilent RF signal generator. The pumped current is measured via a Keithley
6430 Source-Measure unit capable of measuring at the femto-ampere level. DC
voltages are applied to the gates via multiple channel 16 bit DAQs. Unlike the
earlier pump configurations, in which oscillating voltages have been applied
to both the exit and entrance gateFujiwara  _et al._ (2004); Blumenthal  _et
al._ (2007), we follow a later protocolFujiwara, Nishiguchi, and Ono (2008);
Kaestner  _et al._ (2008), in which the potential of the exit gate is held
constant during the pump cycle, and the entrance gate oscillates due to the RF
signal added to the constant potential via a bias T. The gates were fabricated
using electron beam lithography with Ti/Au thermally evaporated. The entrance
finger gate (red) and split gate (green) are used to define the QD entrance
and exit gate respectively. DC voltages were applied to all gates using a
NI2969 cDAQ with the RF signal coupled to the exit gate via a bias tee. The RF
source used was a HP E4400B driving a simple sinusoidal wave.

The device was fabricated on a two-dimensional high mobility GaAs/AlxGa1-xAs
Si-doped electron gas system grown using molecular beam epitaxy. The 2DEG was
formed 90 nm below the surface (10 nm GaAs cap, 40 nm Si-doped
GaAs/AlxGa1-xAs, 40 nm GaAs/AlxGa1-xAs spacer, and GaAs substrate with a
carrier density $n=1.9\times 10^{11}$ cm-2 and mobility $\mu=1.014\times
10^{6}$ cm2/Vs). The 2DEG channel pattern was defined using electron-beam
lithography (EBL) and etched to a depth of 40 nm using wet chemistry. The
split gate has a width and gap of $400$ nm, whilst the finger gate a width of
$150$ nm. The pitch of the gates are $200$ nm. The sample was loaded into a
Leiden Cryogenics dilution fridge with a base temperature of $7$ mK and a $10$
T superconducting magnet. The pumped current was measured using a Keithley
6430 source measure unit (SMU), which was connected to the drain side of the
pump. The source side was grounded.

###  V.2 Dot parameters – analytical consideration

We find the dependence of the quantum dot parameters $W_{d}$ and $D_{d}$ –
width of the potential barrier and the dot depth in a 1D case – shortly after
the system has passed the 0-DIP configuration. The potential profile $E(x)$ is
the sum of the potential profiles of both electrodes:

$E(x)=E_{\rm{ent}}\varphi_{\rm{ent}}+E_{\rm{exit}}\varphi_{\rm{exit}}$ (3)

where $\varphi_{\rm{ent}}$ and $\varphi_{\rm{exit}}$ are the normalised
functions of the potentials creates by each of the electrodes independently,
we assume they are smooth, differentiable, with a simple peak, and
$\max(\varphi_{\rm{ent}})=\max(\varphi_{\rm{exit}})=1$.

In the 0-DIP configuration, $E_{\rm{ent}}=\alpha E_{\rm{exit}}$. If $x_{0}$ is
the inflection point, then the first and second derivatives of the potential
profile in it are equal to zero:

$\begin{dcases}E^{\prime}(x_{0})=\alpha
E_{\rm{exit}}\varphi_{\rm{ent}}^{\prime}(x_{0})+E_{\rm{exit}}\varphi_{\rm{exit}}^{\prime}(x_{0})=0\\\\[5.0pt]
E^{\prime\prime}(x_{0})=\alpha
E_{\rm{exit}}\varphi_{\rm{ent}}^{\prime\prime}(x_{0})+E_{\rm{exit}}\varphi_{\rm{exit}}^{\prime\prime}(x_{0})=0\end{dcases}$
(4)

Shortly past the 0-DIP configuration, $E_{\rm{ent}}=\alpha
E_{\rm{exit}}+\differential E$, the derivatives at the same point are:

$\begin{dcases}E^{\prime}(x_{0})=\differential
E\varphi_{\rm{ent}}^{\prime}(x_{0})\\\\[5.0pt]
E^{\prime\prime}(x_{0})=\differential
E\varphi_{\rm{ent}}^{\prime\prime}(x_{0})\end{dcases}$ (5)

In this case (see Fig.3d), the potential has three special points – the new
inflection point $x_{i}$, the local maximum $x_{+}$ and the local minimum
$x_{-}$. Assuming that $\differential E$ is small and the special points lie
within a small vicinity of the former inflection point $x_{0}$ we find the
distance between the new and the former inflection points $\differential
x_{i}$ in the Taylor expansion. Postulating
$E^{\prime\prime}(x_{0}+\differential x_{i})=0$, we arrive at:

$\differential x_{i}=-\frac{\differential
E\varphi_{\rm{ent}}^{\prime\prime}(x_{0})}{\alpha
E_{\rm{exit}}\varphi_{\rm{ent}}^{\prime\prime\prime}(x_{0})+E_{\rm{exit}}\varphi_{\rm{exit}}^{\prime\prime\prime}(x_{0})}$
(6)

At the new inflection point, to the first order,

$\begin{dcases}E^{\prime}(x_{i})=\differential
E\varphi_{\rm{ent}}^{\prime}(x_{0})\\\\[5.0pt]
E^{\prime\prime}(x_{i})=0\\\\[5.0pt] E^{\prime\prime\prime}(x_{i})=\alpha
E_{\rm{exit}}\varphi_{\rm{ent}}^{\prime\prime\prime}(x_{0})+E_{\rm{exit}}\varphi_{\rm{exit}}^{\prime\prime\prime}(x_{0})\end{dcases}$
(7)

Finally, we find the distance between the new inflection point $x_{i}$ and the
maximum and minimum $x_{+}$ and $x_{-}$. Setting

$E^{\prime}(x_{i}+\differential
x)=E^{\prime}(x_{i})+E^{\prime\prime}(x_{i})\differential
x+E^{\prime\prime\prime}(x_{i})\frac{\differential x^{2}}{2}=0$ (8)

we arrive at:

$\differential x=\sqrt{-\frac{2\differential
E\varphi_{\rm{ent}}^{\prime}(x_{0})}{\alpha
E_{\rm{exit}}\varphi_{\rm{ent}}^{\prime\prime\prime}(x_{0})+E_{\rm{exit}}\varphi_{\rm{exit}}^{\prime\prime\prime}(x_{0})}}$
(9)

We set $2\differential x$, the distance between the local minimum and maximum
as the characteristic size of the dot, which is equal to the width of the
potential barrier; and the depth of the dot, or the height of the barrier can
be found as:

$D_{d}=2E^{\prime}(x_{i})\differential
x+2E^{\prime\prime\prime}(x_{0})\frac{\differential x^{3}}{6}$ (10)

which yields:

$D_{d}=\frac{4}{3}\differential
E\varphi_{\rm{ent}}^{\prime}(x_{0})\sqrt{-\frac{2\differential
E\varphi_{\rm{ent}}^{\prime}(x_{0})}{\alpha
E_{\rm{exit}}\varphi_{\rm{ent}}^{\prime\prime\prime}(x_{0})+E_{\rm{exit}}\varphi_{\rm{exit}}^{\prime\prime\prime}(x_{0})}}$
(11)

The final dependence on the external parameters has the form:

$\begin{dcases}W_{d}=C_{W}\sqrt{\frac{\differential
E}{E_{\rm{exit}}}}\\\\[5.0pt] D_{d}=C_{D}\sqrt{\frac{(\differential
E)^{3}}{E_{\rm{exit}}}}\end{dcases}$ (12)

.

## References

  * Pekola  _et al._ (2013) J. P. Pekola, O.-P. Saira, V. F. Maisi, A. Kemppinen, M. Möttönen, Y. A. Pashkin, and D. V. Averin, Rev. Mod. Phys. 85, 1421 (2013).
  * Keller  _et al._ (1997) M. Keller, J. Martinis, A. Steinbach, and N. Zimmerman, Instrumentation and Measurement, IEEE Transactions on 46, 307 (1997).
  * Bäuerle  _et al._ (2018) C. Bäuerle, D. Christian Glattli, T. Meunier, F. Portier, P. Roche, P. Roulleau, S. Takada, and X. Waintal, Reports Prog. Phys. 81, 056503 (2018), arXiv:1801.07497 .
  * Kaneko, Nakamura, and Okazaki (2016) N.-H. Kaneko, S. Nakamura, and Y. Okazaki, Measurement Science and Technology 27, 032001 (2016).
  * Dupont-Ferrier  _et al._ (2013) E. Dupont-Ferrier, B. Roche, B. Voisin, X. Jehl, R. Wacquez, M. Vinet, M. Sanquer, and S. De Franceschi, Phys. Rev. Lett. 110, 136802 (2013).
  * DiVincenzo and IBM. (2000) D. DiVincenzo and IBM., Fortschritte der Physik 48 (2000), 10.1002/1521-3978(200009)48:9/113.0.CO;2-E.
  * Wright  _et al._ (2009) S. J. Wright, M. D. Blumenthal, M. Pepper, D. Anderson, G. A. C. Jones, C. A. Nicoll, and D. A. Ritchie, Phys. Rev. B 80, 113303 (2009).
  * Maisi  _et al._ (2009) V. Maisi, Y. Pashkin, S. Kafanov, J.-S. Tsai, and J. Pekola, New Journal of Physics 11 (2009), 10.1088/1367-2630/11/11/113057.
  * Giblin  _et al._ (2019) S. P. Giblin, A. Fujiwara, G. Yamahata, M.-H. Bae, N. Kim, A. Rossi, M. Möttönen, and M. Kataoka, Metrologia 56, 044004 (2019).
  * Yamamoto  _et al._ (2012) M. Yamamoto, S. Takada, C. Bäuerle, K. Watanabe, A. D. Wieck, and S. Tarucha, Nature Nanotechnology 7, 247 (2012).
  * Ji  _et al._ (2003) Y. Ji, Y. Chung, D. Sprinzak, M. Heiblum, D. Mahalu, and H. Shtrikman, Nature 422, 415 (2003).
  * Ubbelohde  _et al._ (2015) N. Ubbelohde, F. Hohls, V. Kashcheyevs, T. Wagner, L. Fricke, B. Kästner, K. Pierz, H. W. Schumacher, and R. J. Haug, Nature Nanotechnology 10, 46 (2015).
  * Gumbs  _et al._ (2009) G. Gumbs, D. Huang, M. Blumenthal, S. J. Wright, M. Pepper, and Y. Abranyos, Semiconductor Science and Technology 24 (2009), 10.1088/0268-1242/24/11/115001.
  * Fletcher  _et al._ (2013) J. D. Fletcher, P. See, H. Howe, M. Pepper, S. P. Giblin, J. P. Griffiths, G. A. C. Jones, I. Farrer, D. A. Ritchie, T. J. B. M. Janssen, and M. Kataoka, Phys. Rev. Lett. 111, 216807 (2013).
  * Blumenthal  _et al._ (2023) M. D. Blumenthal, D. Mahony, S. Ahmad, D. Gouveia, H. Howe, H. E. Beere, T. Mitchel, D. A. Ritchie, and M. Pepper, EPJ Quantum Technol. 10, 46 (2023).
  * Ubbelohde  _et al._ (2023) N. Ubbelohde, L. Freise, E. Pavlovska, P. G. Silvestrov, P. Recher, M. Kokainis, G. Barinovs, F. Hohls, T. Weimann, K. Pierz, and V. Kashcheyevs, Nature Nanotechnology 18, 733 (2023).
  * Brange and Flindt (2023) F. Brange and C. Flindt, Nature Nanotechnology 18, 696 (2023).
  * Blumenthal  _et al._ (2007) M. D. Blumenthal, B. Kaestner, L. Li, S. Giblin, T. J. B. M. Janssen, M. Pepper, D. Anderson, G. Jones, and D. A. Ritchie, Nat. Phys. 3, 343 (2007).
  * Kaestner  _et al._ (2008) B. Kaestner, V. Kashcheyevs, S. Amakawa, M. D. Blumenthal, L. Li, T. J. B. M. Janssen, G. Hein, K. Pierz, T. Weimann, U. Siegner, and H. W. Schumacher, Phys. Rev. B 77, 153301 (2008), arXiv:0707.0993 .
  * Kaestner and Kashcheyevs (2015) B. Kaestner and V. Kashcheyevs, Reports Prog. Phys. 78, 103901 (2015), arXiv:1412.7150 .
  * Yamahata, Nishiguchi, and Fujiwara (2014) G. Yamahata, K. Nishiguchi, and A. Fujiwara, Nature communications 5, 5038 (2014).
  * Yamahata and Fujiwara (2024) G. Yamahata and A. Fujiwara, Journal of Applied Physics 135 (2024), 10.1063/5.0179374.
  * Moraru  _et al._ (2007) D. Moraru, Y. Ono, H. Inokawa, and M. Tabe, Phys. Rev. B - Condens. Matter Mater. Phys. 76, 1 (2007).
  * Roche  _et al._ (2013) B. Roche, R.-P. Riwar, B. Voisin, E. Dupont-Ferrier, R. Wacquez, M. Vinet, M. Sanquer, J. Splettstoesser, and X. Jehl, Nature Communications 4, 1581 (2013).
  * Lansbergen, Ono, and Fujiwara (2012) G. P. Lansbergen, Y. Ono, and A. Fujiwara, Nano Letters 12, 763 (2012).
  * Tettamanzi, Wacquez, and Rogge (2014) G. Tettamanzi, R. Wacquez, and S. Rogge, New Journal of Physics 16, 063036 (2014).
  * Kashcheyevs and Kaestner (2010) V. Kashcheyevs and B. Kaestner, Phys. Rev. Lett. 104, 1 (2010), arXiv:0901.4102 .
  * Rossi  _et al._ (2018) A. Rossi, J. Klochan, J. Timoshenko, F. E. Hudson, M. Möttönen, S. Rogge, A. S. Dzurak, V. Kashcheyevs, and G. C. Tettamanzi, Nano Lett. 18, 4141 (2018), arXiv:1803.00791 .
  * Howe  _et al._ (2021) H. Howe, M. Blumenthal, H. E. Beere, T. Mitchell, D. A. Ritchie, and M. Pepper, Appl. Phys. Lett. 119, 153102 (2021).
  * Wright  _et al._ (2008) S. J. Wright, M. D. Blumenthal, G. Gumbs, A. L. Thorn, M. Pepper, T. J. Janssen, S. N. Holmes, D. Anderson, G. A. Jones, C. A. Nicoll, and D. A. Ritchie, Phys. Rev. B - Condens. Matter Mater. Phys. 78, 1 (2008), arXiv:0811.0494 .
  * Kaestner  _et al._ (2009) B. Kaestner, C. Leicht, V. Kashcheyevs, K. Pierz, U. Siegner, and H. W. Schumacher, Appl. Phys. Lett. 94, 012106 (2009), arXiv:0811.1121 .
  * Leicht  _et al._ (2010) C. Leicht, B. Kaestner, V. Kashcheyevs, P. Mirovsky, T. Weimann, K. Pierz, and H. W. Schumacher, Phys. E Low-Dimensional Syst. Nanostructures 42, 911 (2010), arXiv:0909.2778 .
  * Hanief  _et al._ (2024) N. Hanief, M. Blumenthal, H. Howe, H. Beere, T. Mitchell, D. Ritchie, and M. Pepper, Scientific African 24, e02150 (2024).
  * Wright  _et al._ (2010) S. Wright, A. Thorn, M. Blumenthal, S. Giblin, M. Pepper, T. Janssen, M. Kataoka, J. Fletcher, G. Jones, C. Nicoll, G. Gumbs, and D. Ritchie, Journal of Applied Physics - J APPL PHYS 109 (2010), 10.1063/1.3578685.
  * Pyurbeeva, Mol, and Gehring (2022) E. Pyurbeeva, J. A. Mol, and P. Gehring, Chem. Phys. Rev. 3, 041308 (2022), arXiv:2206.05793 .
  * Juergens  _et al._ (2013) S. Juergens, F. Haupt, M. Moskalets, and J. Splettstoesser, Phys. Rev. B - Condens. Matter Mater. Phys. 87, 1 (2013), arXiv:1303.5225 .
  * Hino and Hayakawa (2021) Y. Hino and H. Hayakawa, Phys. Rev. Res. 3, 013187 (2021), arXiv:2003.05567 .
  * Monsel  _et al._ (2022) J. Monsel, J. Schulenborg, T. Baquet, and J. Splettstoesser, Phys. Rev. B 106, 035405 (2022), arXiv:2202.12221 .
  * Monsel, Schulenborg, and Splettstoesser (2023) J. Monsel, J. Schulenborg, and J. Splettstoesser, Eur. Phys. J. Spec. Top. 123 (2023), 10.1140/epjs/s11734-023-00969-4.
  * Pyurbeeva, Thomas, and Mol (2023) E. Pyurbeeva, J. O. Thomas, and J. A. Mol, Mater. Quantum Technol. 3, 025003 (2023), arXiv:2212.10600 .
  * Schulenborg  _et al._ (2023) J. Schulenborg, J. D. Fletcher, M. Kataoka, and J. Splettstoesser,  1, 1 (2023), arXiv:2303.15436 .
  * Kashcheyevs and Timoshenko (2012) V. Kashcheyevs and J. Timoshenko, Phys. Rev. Lett. 109, 216801 (2012), arXiv:1205.3497 .
  * Kashcheyevs and Timoshenko (2014) V. Kashcheyevs and J. Timoshenko, CPEM Dig. (Conference Precis. Electromagn. Meas. , 536 (2014), arXiv:1402.0637 .
  * Potts  _et al._ (1996) A. Potts, R. Shepherd, W. G. Herrenden-Harker, M. Elliott, C. L. Jones, A. Usher, G. A. C. Jones, D. A. Ritchie, E. H. Linfield, and M. Grimshaw, Journal of Physics: Condensed Matter 8, 5189 (1996).
  * Qian  _et al._ (2017) Q. Qian, J. Nakamura, S. Fallahi, G. C. Gardner, J. D. Watson, S. Lüscher, J. A. Folk, G. A. Csáthy, and M. J. Manfra, Phys. Rev. B 96, 035309 (2017).
  * Martin  _et al._ (1988) K. Martin, R. Higgins, J. Rascol, H. Yoo, and J. R. Arthur, Surface Science 196, 323 (1988).
  * Fujiwara  _et al._ (2004) A. Fujiwara, N. M. Zimmerman, Y. Ono, and Y. Takahashi, Appl. Phys. Lett. 84, 1323 (2004).
  * Fujiwara, Nishiguchi, and Ono (2008) A. Fujiwara, K. Nishiguchi, and Y. Ono, Appl. Phys. Lett. 92 (2008), 10.1063/1.2837544.


# Momentum sum rule and factorization of double parton distributions

Krzysztof Golec-Biernat<EMAIL_ADDRESS>Institute of Nuclear Physics Polish
Academy of Sciences, 31-342 Cracow, Poland  Anna M. Staśto<EMAIL_ADDRESS>Penn
State University, University Park, Pennsylvania 16802, United States

###### Abstract

We show that the momentum sum rule is a necessary condition for factorization
of double parton distributions into a product of two single parton
distributions for small values of the parton momentum fractions $x$ and large
enough values of the evolution scale $Q^{2}$. This is a somewhat surprising
result since the momentum sum rule involves integration over all values of the
momentum fraction. In essence, the momentum sum rule provides a proper
relation between the double and single parton distributions, which is
necessary for the small $x$ factorization at large $Q^{2}$.

quantum chromodynamics, multi-parton interactions, parton distributions,
evolution equations, sum rules

##  I Introduction

Multiparton interactions play an important role in high energy scattering of
hadrons Kirschner:1979im ; Shelest:1982dg ; Zinovev:1982be ; Ellis:1982cd ;
Paver:1982yp ; Paver:1983hi ; Paver:1984ux ; Bukhvostov:1985rn ;
Kulesza:1999zh ; Snigirev:2003cq ; Korotkikh:2004bz ; Gaunt:2009re ;
Blok:2010ge ; Ceccopieri:2010kg ; Diehl:2011tt ; Gaunt:2011xd ; Ryskin:2011kk
; Blok:2011bu ; Kom:2011nu ; Diehl:2011yj ; Luszczak:2011zp ; Manohar:2012jr ;
Ryskin:2012qx ; Gaunt:2012dd ; Blok:2013bpa ; Broniowski:2013xba ;
Diehl:2014vaa ; vanHameren:2014ava ; Golec-Biernat:2014bva ;
Ceccopieri:2014ufa ; Maciula:2014pla ; Snigirev:2014eua ; Golec-
Biernat:2014nsa ; Gaunt:2014rua ; Harland-Lang:2014efa ; Blok:2014rza ;
Maciula:2015vza ; Golec-Biernat:2015aza ; Diehl:2015bca ; Broniowski:2016trx ;
Blok:2016lmd ; Rinaldi:2016jvu ; Golec-Biernat:2016vbt ; Ceccopieri:2017oqe ;
Diehl:2017kgu ; Buffing:2017mqm ; Blok:2017alw ; Diehl:2017wew ; Elias:2017flu
; Rinaldi:2018slz ; Rinaldi:2018bsf ; Gaunt:2018eix ; Diehl:2018wfy ;
Diehl:2019rdh ; Broniowski:2019rmu ; Diehl:2020xyg ; Bali:2020mij ;
Cabouat:2020ssr ; Diehl:2021wpp ; Diehl:2021wvd . In particular, the double
parton scattering (DPS) in which two pairs of partons from colliding hadrons
take part in a hard scattering process is of special importance. The DPS
processes were first observed at the Tevatron Akesson:1986iv ; Abe:1997bp ;
Abe:1997xk ; Abazov:2009gc  and are presently studied at the Large Hadron
Collider by the ATLAS Aad:2013bjm ; Aad:2014rua ; ATLAS:2014ofp ;
ATLAS:2016rnd ; ATLAS:2016ydt ; ATLAS:2018zbr , CMS Chatrchyan:2013xxa ;
CMS:2016liw ; CMS:2019jcb ; CMS:2022pio  and LHCb LHCb:2011kri ; LHCb:2012aiv
; LHCb:2015wvu ; LHCb:2016wuo  collaborations.

The computation of DPS cross sections within the collinear framework makes use
of the double parton distribution functions (DPDFs) which obey QCD evolution
equations Kirschner:1979im ; Shelest:1982dg ; Zinovev:1982be ; Snigirev:2003cq
; Korotkikh:2004bz ; Ceccopieri:2010kg ; Ceccopieri:2014ufa ; Gaunt:2009re ,
similar to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution
equations for the single parton distribution functions (PDFs), see also
Ryskin:2012qx  for a pedagogical presentation and Diehl:2017wew  for an
overview. The evolution equations for DPDFs conserve sum rules which relate
the double and single parton distribution functions. It means that once these
rules are imposed on initial conditions for the evolution equations at an
initial scale $Q^{2}_{0}$, they are also obeyed by evolved distributions at
the scale $Q^{2}\geq Q^{2}_{0}$.

All the attempts to construct conditions which exactly satisfy the new sum
rules were rather unsuccessful until now; see, e.g., Refs. Korotkikh:2004bz ;
Gaunt:2011xd ; Golec-Biernat:2014bva  with an exception of the analysis
Broniowski:2013xba  for valence quarks only. Also in a pure gluonic case, the
double gluon distribution $D_{gg}$ was proposed in Ref. Golec-Biernat:2015aza
which obeys the momentum sum rule relating $D_{gg}$ and the single gluon
distribution $D_{g}$. This was achieved due to a particular form of $D_{g}$,
used in global fits to hard scattering data as an initial condition for the
DGLAP evolution equations. The parameters of $D_{g}$ fully determined the
parameters of the initial $D_{gg}$ such that the momentum sum rule is
fulfilled.

In most practical applications of DPDFs it is usually assumed that for small
parton momentum fractions, $x_{1},x_{2}\ll 1$, the DPDFs factorize into a
product of single PDFs, e.g. $D_{gg}(x_{1},x_{2},Q^{2})\approx
D_{g}(x_{1},Q^{2})D_{g}(x_{2},Q^{2})$. However, the example of the initial
condition for $D_{gg}$ from Golec-Biernat:2015aza  shows that the small $x$
factorization can be strongly violated at the initial scale
$Q_{0}^{2}=1~{}{\rm GeV}^{2}$ but is restored after the evolution to a
sufficiently large $Q^{2}$.

Therefore, it is a main goal of this paper to elucidate the issue of the small
$x$ factorization in the framework of the QCD evolution equations for DPDFs.
In particular, we will show that the momentum sum rule, which in principle
must be imposed on the initial distributions for the QCD evolution equations,
is a necessary condition for the small $x$ factorization, i.e., without the
momentum sum rule the small $x$ factorization is violated. This result has
important phenomenological consequences for the current prescriptions for the
initial conditions for DPDFs which use the information on the single PDFs.
Such parton distributions usually fulfill the momentum sum rule only
approximately and must be used with some care in precise studies based on the
QCD evolution equations.

It should be mentioned that the issue of the small $x$ factorization is
studied here for the DPDFs with transverse momentum ${\mathbf{q}}=0$, when the
momentum sum rule is valid. However, there are strong indications that for
${\mathbf{q}}\neq 0$ there is no factorization at all due to correlations
between partons in the impact parameter space with the variable
${\mathbf{b}}$, being Fourier conjugate to ${\mathbf{q}}$ Diehl:2014vaa ;
Rinaldi:2016jvu ; Diehl:2020xyg . Nevertheless, our studies are important
since the form of the DPDFs at ${\mathbf{q}}=0$ is essential for the
generalization to the ${\mathbf{q}}\neq 0$ case.

The paper is organized as follows. In Sec. II we present main information
about the DPDFs, their QCD evolution equations and sum rules. In Sec. III we
restrict ourselves to the pure gluon case with the double gluon distribution
$D_{gg}$ for a detailed numerical analysis presented in the next sections. In
Sec. IV, based on a simple ansatz for the single and double gluon
distributions, we present two numerical examples showing that the small $x$
factorization of $D_{gg}$ is conserved or restored in the QCD evolution. In
Sec. V we illustrate on numerical examples that the momentum sum rule is a
necessary condition for the latter conclusions. In Sec. VI we provide an
analytical insight to these results, using the Mellin moment representation.
Finally, in Sec. VII we present the current situation concerning the relation
between the double and single parton distributions, including also quarks into
the considerations. In Appendix A the proof of the momentum sum rule
conservation by the evolution equations is presented, while in Appendix B the
asymptotic solution to the evolution equations is discussed in the Mellin
moment space.

##  II Evolution equations and sum rules

The double parton distributions
$\bar{D}_{f_{1}f_{2}}(x_{1},x_{2},{\mathbf{b}})$ have probabilistic
interpretation of the number density of pairs of partons with longitudinal
momentum fractions $x_{1}$ and $x_{2}$ (for which $0<x_{1}+x_{2}\leq 1$) at a
relative transverse vector ${\mathbf{b}}$ Diehl:2011tt ; Diehl:2011yj . Parton
flavours (including gluon) are denoted by $f_{1}$ and $f_{2}$. In this paper,
we are interested in the double parton distributions integrated over
${\mathbf{b}}$, which are equal to the Fourier transform,

$D_{f_{1}f_{2}}(x_{1},x_{2},{\mathbf{q}})=\int
d^{2}{\mathbf{b}}\,{\rm{e}}^{i{\mathbf{q}}\cdot{\mathbf{b}}}\,\bar{D}_{f_{1}f_{2}}(x_{1},x_{2},{\mathbf{b}})\;,$
(1)

taken at ${\mathbf{q}}=0$,

$D_{f_{1}f_{2}}(x_{1},x_{2},{\mathbf{q}}=0)=\int
d^{2}{\mathbf{b}}\,\bar{D}_{f_{1}f_{2}}(x_{1},x_{2},{\mathbf{b}})\;.$ (2)

For the physical meaning of the transverse momentum ${\mathbf{q}}$, see
Diehl:2011yj . The DPDFs (2) obey QCD evolution equations with respect to two
hard scales, $Q_{1}^{2}$ and $Q^{2}_{2}$, which are also present due to
renormalization. We consider the DPDFs with equal hard scales,
$Q_{1}^{2}=Q_{2}^{2}\equiv Q^{2}$, i.e.

$D_{f_{1}f_{2}}(x_{1},x_{2},Q^{2})\,\equiv\,D_{f_{1}f_{2}}(x_{1},x_{2},{\mathbf{q}}=0,Q^{2},Q^{2})\;.$
(3)

Introducing the evolution parameter

$t=t(Q^{2})=\int^{Q^{2}}_{Q_{0}^{2}}\frac{\alpha_{s}(\mu^{2})}{2\pi}\frac{d\mu^{2}}{\mu^{2}}=\frac{6}{33-2n_{f}}\ln\frac{\ln(Q^{2}/\Lambda_{\rm
QCD}^{2})}{\ln(Q_{0}^{2}/\Lambda_{\rm QCD}^{2})}\;,$ (4)

where $\alpha_{s}(\mu^{2})$ is the running strong coupling constant in the
leading order (LO) approximation, the evolution equations read Shelest:1982dg
; Zinovev:1982be ; Snigirev:2003cq ; Gaunt:2009re ; Ceccopieri:2010kg ;
Ryskin:2011kk

$\displaystyle\frac{\partial}{\partial
t}\,D_{f_{1}f_{2}}(x_{1},x_{2},t)=\sum_{f^{\prime}}\Bigg{\\{}$
$\displaystyle\int^{1-x_{2}}_{x_{1}}\frac{du}{u}\,{{P}}_{f_{1}f^{\prime}}\\!\left(\frac{x_{1}}{u}\right)D_{f^{\prime}f_{2}}(u,x_{2},t)+\int_{x_{2}}^{1-x_{1}}\frac{du}{u}\,{{P}}_{f_{2}f^{\prime}}\\!\left(\frac{x_{2}}{u}\right)D_{f_{1}f^{\prime}}(x_{1},u,t)$
$\displaystyle+\frac{1}{x_{1}+x_{2}}\,{P}_{f^{\prime}\to
f_{1}f_{2}}\\!\left(\frac{x_{1}}{x_{1}+x_{2}}\right)D_{f^{\prime}}(x_{1}+x_{2},t)\Bigg{\\}}\;.$
(5)

Here, the functions ${P}$ on the rhs are the LO Altarelli-Parisi splitting
functions (with virtual corrections included) and the summation is performed
over quark/antiquark flavors and gluon. The third, splitting term on the rhs
corresponds to the splitting of one parton into two daughter partons,
described by the LO Altarelli-Parisi splitting function for the real emission,
$P_{f^{\prime}\to f_{1}f_{2}}=P^{R}_{f^{\prime}f_{1}}$. It also contains the
single PDFs; thus Eq. (II) has to be solved together with the ordinary DGLAP
equations

$\frac{\partial}{\partial
t}\,D_{f}(x,t)=\sum_{f^{\prime}}\int^{1}_{x}\frac{du}{u}\,{{P}}_{ff^{\prime}}\\!\left(\frac{x}{u}\right)D_{f^{\prime}}(u,t)\;.$
(6)

To formulate the next-to-leading order (NLO) evolution equations for the
DPDFs, the two-loop splitting functions ${P}_{f^{\prime}\to f_{1}f_{2}}$ were
calculated in Diehl:2019rdh ; Diehl:2021wpp . The NLO formulation, however, is
beyond the scope of this paper.

The DPDFs for ${\mathbf{q}}=0$ obey sum rules which can be derived starting
from the definition of DPDFs based on the light-cone nucleon wave function
Gaunt:thesis . What is most important, these sum rules are consistent with the
evolution equations, which means that once the sum rules are assumed at the
initial scale $t=0$, they are conserved during the QCD evolution to any $t$
Gaunt:2009re ; Diehl:2011tt ; Gaunt:2011xd . These are: the momentum sum rule

$\displaystyle\sum_{f_{1}}\int_{0}^{1-x_{2}}dx_{1}\,x_{1}D_{f_{1}f_{2}}(x_{1},x_{2},t)=(1-x_{2})D_{f_{2}}(x_{2},t)\;,$
(7)

and the valence quark number sum rule

$\displaystyle\int_{0}^{1-x_{2}}dx_{1}\\!\left\\{D_{q_{i}f_{2}}(x_{1},x_{2},t)-D_{{\overline{q}}_{i}f_{2}}(x_{1},x_{2},t)\right\\}=(N_{q_{i}}-\delta_{q_{i}f_{2}}+\delta_{{\overline{q}}_{i}f_{2}})D_{f_{2}}(x_{2},t)\,,$
(8)

where $q_{i}=u,d,s$ and $N_{u}=2,N_{d}=1,N_{s}=0$ are the corresponding
valence quark numbers. Their form can be readily understood from the
probability theory, treating the ratios of the double to single distributions
as conditional probabilities, see Gaunt:2009re  for details. Analogous sum
rules also hold with respect to the second parton momentum $x_{2}$. To ensure
this property, the parton exchange symmetry must be imposed on the initial
conditions, which is kept by the evolution to any $t$,

$D_{f_{1}f_{2}}(x_{1},x_{2},0)=D_{f_{2}f_{1}}(x_{2},x_{1},0)\;.$ (9)

Relations (7) and (8) should be considered together with the momentum sum rule
for single PDFs,

$\displaystyle\sum_{f}\int_{0}^{1}dx\,xD_{f}(x,t)=1\;,$ (10)

and the valence quark number sum rule

$\displaystyle\int_{0}^{1}dx\,\left\\{D_{q_{i}}(x,t)-D_{{\overline{q}}_{i}}(x,t)\right\\}=N_{q_{i}}\,,$
(11)

which are also conserved by the DGLAP evolution equations. In Appendix A we
present the proof of the conservation of the momentum sum rule (7) by the
evolution equations (II) and (6), in which the presence of the splitting term
in Eq. (II) plays the crucial role.

##  III Pure gluonic case

In the following, we shall consider only the gluonic case for the purposes of
simplicity of the presentation. In this case, we only deal with the double
gluon distribution $D_{gg}$ and the single gluon distribution $D_{g}$, which
obey the evolution equations (II) reduced to the gluon sector,

$\displaystyle\frac{\partial}{\partial t}\,D_{gg}(x_{1},x_{2},t)=$
$\displaystyle\int^{1-x_{2}}_{x_{1}}\frac{du}{u}\,{{P}}_{gg}\\!\left(\frac{x_{1}}{u}\right)D_{gg}(u,x_{2},t)+\int_{x_{2}}^{1-x_{1}}\frac{du}{u}\,{{P}}_{gg}\\!\left(\frac{x_{2}}{u}\right)D_{gg}(x_{1},u,t)$
$\displaystyle+\frac{1}{x_{1}+x_{2}}\,P^{R}_{gg}\\!\left(\frac{x_{1}}{x_{1}+x_{2}}\right)D_{g}(x_{1}+x_{2},t)\;,$
(12)

and

$\frac{\partial}{\partial
t}\,D_{g}(x,t)=\int_{x}^{1}\frac{dz}{z}\,{P}_{gg}\\!\left(\frac{x}{z}\right)D_{g}(z,t)\;.$
(13)

The momentum sum rules (7) and (10), conserved by the above evolution
equations, take the following form

$\int_{0}^{1-x_{2}}dx_{1}\,x_{1}D_{gg}(x_{1},x_{2},t)=(1-x_{2})D_{g}(x_{2},t)\;,$
(14)

and

$\int_{0}^{1}dx\,xD_{g}(x,t)=1\;.$ (15)

With the normalization to one in the last condition, we assume that gluons
carry all the proton longitudinal momentum. The momentum momentum sum rules
should be valid for any $t$, including $t=0$ where the initial conditions for
the evolution equations (III) and (13) are specified. As in the general case,
we assume the gluon exchange symmetry for the initial conditions to ensure the
momentum sum rule (14) with respect to the second momentum fraction,

$D_{gg}(x_{1},x_{2},0)=D_{gg}(x_{2},x_{1},0)\;.$ (16)

The conservation of the momentum sum rules by the evolution equations means
that the single gluon distribution obtained from (14),

$D_{g}(x_{2},t)=\frac{1}{1-x_{2}}\int_{0}^{1-x_{2}}dx_{1}\,x_{1}D_{gg}(x_{1},x_{2},t)\;,$
(17)

obeys equation (13) for any $t$, including $t=0$. This imposes strong
constraint on the double gluon distribution alone which results from the
momentum sum rule (15),

$\int_{0}^{1}dx_{2}\int_{0}^{1-x_{2}}dx_{1}\,\frac{x_{1}x_{2}}{1-x_{2}}\,D_{gg}(x_{1},x_{2},t)=1\;.$
(18)

In the next section, we will consider a simple ansatz for the initial gluon
distributions which satisfies the momentum sum rules (14) and (15), which will
allow for detailed studies of the relation between the momentum sum rule and
the small $x$ factorization of $D_{gg}(x_{1},x_{2},t)$.

##  IV A simple example

Figure 1: In the upper plots: $x_{1}x_{2}D_{gg}(x_{1},x_{2})$ as a function of
$x_{1}$ at the initial scale $Q^{2}=1~{}{\rm GeV}^{2}$ (left panels) and final
scale $Q^{2}=100~{}{\rm GeV}^{2}$ (right panels) and fixed $x_{2}=10^{-2}$ for
model (23) (solid lines), Gaunt model (27) (dashed lines) and fully factorized
ansatz (28) (dash-dotted lines). The parameters in (23) are: $\alpha_{g}=-1$
and $\beta_{g}=2.5$. In the bottom plots: the corresponding ratios (26). Both
the momentum sum rule (14) and the small $x$ factorization (25) hold true for
our model.

A general construction for the initial conditions for the single and double
gluon distributions which obey the sum rules was proposed in Golec-
Biernat:2015aza . This framework is based on sums over Dirichlet
distributions, where relations between the powers and the normalizations can
be found to ensure that the sum rules for double and single distributions are
simultaneously satisfied. In the simple example below we shall follow the
construction of Ref. Golec-Biernat:2015aza  for the simple case of the
distributions having one term in the sum.

Let us consider the single gluon distribution of the form

$D_{g}(x,0)=A_{g}\,x^{\alpha_{g}}(1-x)^{\beta_{g}}\;,$ (19)

where the normalization constant is determined from the momentum sum rule (15)

$A_{g}\int_{0}^{1}dx\,x^{\alpha_{g}+1}(1-x)^{\beta_{g}}=1\;.$ (20)

Using the well-known formula relating the Euler beta function to the gamma
functions,

$\int_{0}^{1}dx\,x^{\alpha}(1-x)^{\beta}=\frac{\Gamma(\alpha+1)\Gamma(\beta+1)}{\Gamma(\alpha+\beta+2)}\;,$
(21)

we obtain

$A_{g}=\frac{\Gamma(\alpha_{g}+\beta_{g}+3)}{\Gamma(\alpha_{g}+2)\Gamma(\beta_{g}+1)}\;.$
(22)

It can be shown Golec-Biernat:2015aza  that the double gluon distribution
which obeys the momentum sum rule (14) with $D_{g}$ of the form (19) is given
by

$D_{gg}(x_{1},x_{2},0)=A_{gg}\,(x_{1}x_{2})^{\alpha_{g}}(1-x_{1}-x_{2})^{\beta_{g}-\alpha_{g}-1}\;,$
(23)

where the normalization constant

$A_{gg}=\frac{\Gamma(\beta_{g}+2)}{\Gamma(\alpha_{g}+2)\Gamma(\beta_{g}-\alpha_{g})}\,A_{g}\;.$
(24)

This is a remarkable result since the parameters of the double distribution
$D_{gg}$ are completely determined by the parameters of the single
distribution $D_{g}$ through the momentum sum rule (14). Note also, that the
powers governing the small $x_{1,2}$ behavior are the same in Eqs. (19) and
(23).

Figure 2: The same as in Fig. 1 but for $\alpha_{g}=-0.5$ and $\beta_{g}=2.5$
in (23). The momentum sum rule is valid but the small $x$ factorization is
violated by the initial conditions (solid line in the left bottom plot). After
the evolution to large $Q^{2}=100\,\rm GeV^{2}$, the small $x$ factorization
is to a good approximation restored (solid line on the right bottom plot).

For $\alpha_{g}=-1$ we have $A_{gg}=A_{g}^{2}$ and the double gluon
distribution factorizes for $x_{1},x_{2}\ll 1$,

$D_{gg}(x_{1},x_{2},0)\approx A_{g}^{2}\,(x_{1}x_{2})^{\alpha_{g}}\,\approx
D_{g}(x_{1},0)D_{g}(x_{2},0)\;.$ (25)

We illustrate this case in Fig. 1 choosing $\alpha_{g}=-1$ and $\beta_{g}=2.5$
in Eq. (19) for the single gluon distribution defined at the initial scale
$Q_{0}^{2}=1~{}{\rm GeV}^{2}$. The choice of $\alpha_{g}=-1$ corresponds to
“flat” initial distribution $xD_{g}(x,0)$ at small values of $x$. The
corresponding initial double gluon distribution was computed from Eq. (23). In
the two upper plots we show the distribution $x_{1}x_{2}D_{gg}(x_{1},x_{2},t)$
at the initial scale and final scale equal to $Q^{2}=100~{}{\rm GeV}^{2}$
(solid lines). Both plots are presented in the way appropriate for the studies
of the small $x$ limit of the double gluon distributions. Namely, we plot them
as functions of $x_{1}$ for fixed $x_{2}=10^{-2}$ to study the limit $x_{1}\to
0$. The two bottom plots show the ratio

${\rm Ratio}=\frac{D_{gg}(x_{1},x_{2},t)}{D_{g}(x_{1},t)D_{g}(x_{2},t)}\;,$
(26)

at the corresponding values of $Q^{2}$ [related to $t$ by Eq. (4)]. We see
that both before and after the evolution, the small $x$ factorization holds
for the double gluon distribution to a good approximation up to
$x_{1}=10^{-1}$.

We also show the results obtained for the Gaunt-Stirling (GS) prescription
Gaunt:2009re  for the initial condition (dashed lines),

$D_{gg}(x_{1},x_{2},0)=D_{g}(x_{1},0)D_{g}(x_{2},0)\frac{(1-x_{1}-x_{2})^{2}}{(1-x_{1})^{2}(1-x_{2})^{2}}\;,$
(27)

as well as for the fully factorized ansatz (dash-dotted lines)

$D_{gg}(x_{1},x_{2},0)=D_{g}(x_{1},0)\,D_{g}(x_{2},0)\,\theta(1-x_{1}-x_{2})\;,$
(28)

where the single gluon distribution (19) with the assumed parameters was used.
Both prescriptions give the approximate factorization after the evolution to
$Q^{2}=100~{}{\rm GeV}^{2}$, although the result for the fully factorized
ansatz is slightly worse. This agreement can be attributed to a rather weak
violation of the momentum sum rule by these distributions. In particular, the
integral (18) computed for the ansatz (27) gives $1.04$ while for the ansatz
(28) we find $1.3$ (instead of 1 as in our model).

For $\alpha_{g}\neq-1$, the small $x$ factorization of the initial double
gluon distribution (23) is no longer true,

$D_{gg}(x_{1},x_{2},0)\neq D_{g}(x_{1},0)D_{g}(x_{2},0)\;,$ (29)

but the momentum sum rule is still fulfilled. This means that the momentum sum
rule is neither necessary nor sufficient condition for the small $x$
factorization of the initial double gluon distribution. We illustrate this
fact in Fig. 2 with the initial double gluon distribution (23) with the
parameters $\alpha_{g}=-0.5$ and $\beta_{g}=2.5$ (solid line on the left
bottom plot). However, it is very interesting that the small $x$ factorization
of $D_{gg}(x_{1},x_{2})$ is to a good approximation restored by the evolution
to $Q^{2}=100~{}{\rm GeV}^{2}$ since the ratio (26) is close to one (solid
line on the right bottom plot).

The same effect is also clearly visible in the analysis in Ref. Golec-
Biernat:2015aza  with a realistic ansatz based on a particular form of the
single PDFs in the MSTW08 parametrization Martin:2009iq ,

$D_{f}(x,0)=\sum_{i=1}^{3}A_{f}^{i}\,x^{\alpha_{f}^{i}}(1-x)^{\beta_{f}^{i}}\;,$
(30)

where $f$ denotes quark/antiquark flavor or gluon. In the pure gluonic case,
the double gluon distribution is fully determined by the parameters of the
single gluon distribution $D_{g}(x,0)$, using the momentum sum rule (14),

$D_{gg}(x_{1},x_{2},0)=\sum_{i=1}^{3}A_{gg}^{i}\,(x_{1}x_{2})^{\alpha_{g}^{i}}(1-x_{1}-x_{2})^{\beta_{g}^{i}-\alpha_{g}^{i}-1}\;,$
(31)

where $A_{gg}^{i}$ is also given in terms of the parameters $\alpha_{g}^{i}$
and $\beta_{g}^{i}$. The initial condition (31) strongly violates the small
$x$ factorization at the scale $Q_{0}^{2}=1~{}{\rm GeV}^{2}$. However, after
the evolution to $Q^{2}=100~{}{\rm GeV}^{2}$ the small $x$ factorization is
restored, see in Fig. 1 in Ref. Golec-Biernat:2015aza . In the next section,
we will explore this phenomenon in more detail.

##  V Momentum sum rule violation

Figure 3: The same as in Fig. 1 but for $\eta=5$ in the initial condition (33)
(solid lines). The momentum sum rule is no longer valid but the small $x$
factorization holds true (solid line on the left bottom plot), while it is
violated after the evolution (solid line on the right bottom plot).

In order to understand the role of the momentum sum rule conservation for the
restoration of the small $x$ factorization, we consider the single gluon
distribution (19) with normalization (22),

$D_{g}(x,0)=A_{g}\,x^{\alpha_{g}}(1-x)^{\beta_{g}}\;,$ (32)

while the double gluon distribution (23) with normalization (24) is modified
by changing the parameter $\eta=\beta_{g}-\alpha_{g}-1$, which controls the
large $x$ behaviour, i.e.,

$D^{(m)}_{gg}(x_{1},x_{2},0)=A_{gg}\,(x_{1}x_{2})^{\alpha_{g}}(1-x_{1}-x_{2})^{\eta}\;,$
(33)

where $\eta$ is now arbitrary. The single gluon distribution still obeys the
momentum sum rule (15) but relation (14) is violated. Indeed, by computing the
gluon distribution from relation (17),

$D_{g}^{(m)}(x_{2},0)=\frac{1}{1-x_{2}}\int_{0}^{1-x_{2}}dx_{1}\,x_{1}D_{gg}^{(m)}(x_{1},x_{2},0)\;,$
(34)

we obtain

$D^{(m)}_{g}(x,0)=A_{g}\,\frac{\Gamma(\beta_{g}+2)\Gamma(\eta+1)}{\Gamma(\alpha_{g}+\eta+3)\Gamma(\beta_{g}-\alpha_{g})}\,x^{\alpha_{g}}\,(1-x)^{\alpha_{g}+\eta+1}\;,$
(35)

Thus, only for $\eta=\beta_{g}-\alpha_{g}-1$ we find
$D^{(m)}_{g}(x,0)=D_{g}(x,0)$ and the momentum sum rule (14) is fulfilled. For
any other $\eta$, this rule is violated by the initial distributions
$D_{gg}^{(m)}(x_{1},x_{2},0)$ and $D_{g}(x,0)$.

We illustrate this situation in Fig. 3 for the choice of the parameters
$\alpha_{g}=-1$, $\beta_{g}=2.5$ and $\eta=5$. Since
$\beta_{g}-\alpha_{g}-1=2.5\neq 5$ we deal with the momentum sum rule
violation. Nevertheless, the small $x$ factorization holds at the initial
scale, which is shown on the left bottom plot in Fig. 3 by the solid line,
being close to one for $x_{1}<10^{-2}$. However, the factorization is strongly
violated after the evolution to $Q^{2}=100~{}{\rm GeV}^{2}$ (solid line on the
right bottom plot) in the same range of $x_{1}$,

$\frac{D_{gg}^{(m)}(x_{1},x_{2},t)}{D_{g}(x_{1},t)D_{g}(x_{2},t)}\approx
0.6-0.7\;.$ (36)

This can be attributed to the mismatch between the single gluon distributions,
$D_{g}(x,t)$ and $D_{g}^{(m)}(x,t)$. We numerically checked that if
$D_{g}^{(m)}(x,t)$ is substituted to the ratio (36), after the additional
rescaling of $D^{(m)}_{gg}(x_{1},x_{2},0)$ such that the momentum sum rule
(15) is fulfilled by $D_{g}^{(m)}(x,0)$, we find factorization at the scale
$Q^{2}=100~{}{\rm GeV}^{2}$. The reason for the additional rescaling will
become clear from the analytical insight presented in the next section.

To summarize the presented example, the violation of the momentum sum rule
(14) by the initial conditions leads to a mismatch between the single and
double gluon distributions, which is manifested in the violation of the small
$x$ factorization even at high scales. On the other hand, when the momentum
sum rules (14) and (15) are fulfilled, then the approximate restoration of the
small $x$ factorization is observed after the evolution to high scales.

##  VI Mellin moment formulation

In order to give an analytical insight into the relation between the momentum
sum rule and small $x$ factorization, let us introduce the Mellin moments for
a single gluon distribution,

$\tilde{D}_{g}(n,t)=\int_{0}^{1}dx\,x^{n-1}D_{g}(x,t)\;,$ (37)

and for a double gluon distribution,

$\tilde{D}_{gg}(n_{1},n_{2},t)=\int_{0}^{1}dx_{1}\int_{0}^{1}dx_{2}\,x_{1}^{n_{1}-1}x_{2}^{n_{2}-1}\,\theta(1-x_{1}-x_{2})\,D_{gg}(x_{1},x_{2},t)\;.$
(38)

The momentum sum rule (14) in terms of the Mellin moments is given by

$\tilde{D}_{gg}(n_{1},2,t)=\tilde{D}_{g}(n_{1},t)-\tilde{D}_{g}(n_{1}+1,t)\;,$
(39)

while the momentum sum (15) reads as

$\tilde{D}_{g}(2,t)=1\;.$ (40)

It was shown in Golec-Biernat:2014nsa  that the solution to the evolution
equations (III) and (13) in the Mellin moment space reads as

$\displaystyle\tilde{D}_{gg}(n_{1},n_{2},t)$
$\displaystyle={\rm{e}}^{\gamma(n_{1})t+\gamma(n_{2})t}\,\tilde{D}_{gg}(n_{1},n_{2},0)+\int_{0}^{t}dt^{\prime}{\rm{e}}^{\gamma(n_{1})(t-t^{\prime})+\gamma(n_{2})(t-t^{\prime})}\,\tilde{\gamma}(n_{1},n_{2})\,\tilde{D}_{g}(n_{1}+n_{2}-1,t^{\prime})\;,$
(41)

where the gluonic anomalous dimension equals

$\gamma(n)=\int_{0}^{1}dx\,x^{n-1}\,{{P}}_{gg}(x)\;,$ (42)

while

$\tilde{\gamma}(n_{1},n_{2})=\int_{0}^{1}dx\,x^{n_{1}-1}(1-x)^{n_{2}-1}\,P^{R}_{gg}(x)\;.$
(43)

In Appendix B we show that in the small $x$ limit, when $x_{1}\to 0$ and
$x_{2}={\rm fixed}$ and small (e.g. $x_{2}=10^{-2}$ as in the presented
analysis), the solution for $t\to\infty$ is given in terms of the Mellin
moments by

$\tilde{D}_{gg}(n_{1},n_{2},t)\simeq
e^{\gamma(n_{1})t+\gamma(n_{2})t}\left(\tilde{D}_{gg}(n_{1},n_{2},0)+\tilde{D}_{g}(n_{1}+n_{2}-1,0)\right)\;.$
(44)

where $(n_{1}-1)\to 0$ and $(n_{2}-1)$ is finite.

This needs to be compared with the numerically found result which gives us
approximate factorization at low $x$

$\tilde{D}_{gg}(n_{1},n_{2},t)\simeq
e^{\gamma(n_{1})t+\gamma(n_{2})t}\;\tilde{D}_{g}(n_{1},0)\,\tilde{D}_{g}(n_{2},0)\;.$
(45)

So the small $x$ factorization holds if we have

$\tilde{D}_{gg}(n_{1},n_{2},0)+\tilde{D}_{g}(n_{1}+n_{2}-1,0)=\tilde{D}_{g}(n_{1},0)\,\tilde{D}_{g}(n_{2},0)\;.$
(46)

For $n_{2}=2$ we find the following relation

$\tilde{D}_{gg}(n_{1},2,0)+\tilde{D}_{g}(n_{1}+1,0)=\tilde{D}_{g}(n_{1},0)\,\tilde{D}_{g}(2,0)\;,$
(47)

Assuming the momentum sum rule for the single gluon distribution,
$D_{g}(2,0)=1$, we obtain

$\tilde{D}_{gg}(n_{1},2,0)=\tilde{D}_{g}(n_{1},0)-\tilde{D}_{g}(n_{1}+1,0)\;,$
(48)

which is the momentum sum rule (39) to be satisfied by the initial condition.
This is a necessary condition for the small $x$ factorization at sufficiently
large $t$. Note that to arrive at this conclusion, the momentum sum rule (40)
has to be fulfilled, which explains the numerical observation presented in the
previous section.

##  VII General case

Figure 4: In the upper plots: the double gluon distributions
$x_{1}x_{2}D_{gg}(x_{1},x_{2})$ as a function of $x_{1}$ at
$Q^{2}=1,100~{}{\rm GeV}^{2}$ and fixed $x_{2}=10^{-2}$ for the fully
factorized prescription (51) (dash-dotted lines) and the GS ansatz (50)
(dashed lines) with the initial conditions given by the MSTW08 PDFs. The
corresponding ratios (26) are shown in the bottom plots. The evolution
equations with quarks were used for the results at $Q^{2}=100~{}{\rm
GeV}^{2}$.

Let us consider the case with quarks governed by the evolution equations and
sum rules presented in Sec. II. The momentum sum rules (7) and (10) lead to
the following consistency condition which should be satisfied by the initial
double parton distributions $D_{f_{1},f_{2}}(x_{1},x_{2},0)$:

$\sum_{f_{1},f_{2}}\int_{0}^{1}dx_{2}\int_{0}^{1-x_{2}}dx_{1}\,\frac{x_{1}x_{2}}{1-x_{2}}\,D_{f_{1}f_{2}}(x_{1},x_{2},0)=1\;,$
(49)

where the summation over quark/antiquark flavors and gluon is performed. The
departure from 1 of the integral on the lhs is the indication that the
momentum sum rule (7) is violated.

In Fig. 4 we analyze the small $x$ factorization for the double gluon
distribution $x_{1}x_{2}D_{gg}(x_{1},x_{2},t)$ obtained from the full
evolution equations (II) with the initial conditions given the Gaunt-Stirling
(GS) ansatz (dashed lines)

$D_{f_{1}f_{2}}(x_{1},x_{2},0)=D_{f_{1}}(x_{1},0)\,D_{f_{2}}(x_{2},0)\,\frac{(1-x_{1}-x_{2})^{2}}{(1-x_{1})^{2+\alpha(f_{1})}(1-x_{2})^{2+\alpha(f_{2})}}\;,$
(50)

where $\alpha(f)=0$ for a sea quark and gluon and $\alpha(f)=0.5$ for valence
quark, and with the fully factorized (FF) ansatz (dot-dashed lines)

$D_{f_{1}f_{2}}(x_{1},x_{2},0)=D_{f_{1}}(x_{1},0)\,D_{f_{2}}(x_{2},0)\,\theta(1-x_{1}-x_{2})\;.$
(51)

where the single PDFs are given by the MSTW08 parametrization Martin:2009iq .
The main motivation for such initial distributions is the use of single PDFs
which are well determined from global fits to data. This allows to constrain
the DPDFs which presently cannot be experimentally determined due to the
scarcity of data on double parton scattering. In both cases, condition (49) is
only approximately fulfilled with 1.06 and 1.30 on the r.h.s for the GS and FF
prescriptions, respectively, which means that the momentum sum rule for the
double parton distributions is not exactly fulfilled. As a result, the small
$x$ factorization, valid at the initial scale (bottom left plot), is only
approximate at the scale $Q^{2}=100~{}{\rm GeV}^{2}$ (right bottom plot). We
see that the GS initial condition leads to a slightly better factorization.
For a more detailed analysis of this ansatz, see Ref. Gaunt:2009re .

The obvious drawback of the above initial conditions is the sum rules
violation, which in principle makes the QCD evolution inconsistent. In Ref.
Golec-Biernat:2015aza  this problem was addressed in the pure gluonic case by
taking advantage of a particular form of the single PDFs in the MSTW08
parametrization. In principle, the entire construction can be repeated for the
case with quarks, obtaining the double parton distributions satisfying both
the momentum and valence quark number rules. However, such an extension
requires additional constraints imposed on the parameters of the single PDFs,
which unfortunately are not present in realistic single PDFs obtained from
fits to data. Therefore, the GS proposal, which features factorization at
small $x$, is at the moment the best practical way to deal with the initial
conditions for the DPDFs evolution which incorporates the knowledge about the
single PDFs.

It should be emphasized at the end that the small $x$ factorization of the
double parton distributions was discussed here only for the transverse
momentum ${\mathbf{q}}=0$ when the sum rules are valid, see Section II. For
${\mathbf{q}}\neq 0$, there are indications that the factorization does not
hold due to correlations in the ${\mathbf{b}}-$space Blok:2011bu ;
Blok:2013bpa . In particular, after the Fourier transformation to the
${\mathbf{q}}-$space, the following ansatz was proposed in Diehl:2014vaa  for
$x_{1},x_{2}<0.1$,

$D_{f_{1}f_{2}}(x_{1},x_{2},{\mathbf{q}},0)=D_{f_{1}}(x_{1},0)\,D_{f_{2}}(x_{2},0)\exp\\!\left\\{-h_{f_{1}f_{2}}(x_{1},x_{2})\,{\mathbf{q}}^{2}\right\\},$
(52)

where $h_{f_{1}f_{2}}(x_{1},x_{2})$ is some function. Thus, for
${\mathbf{q}}=0$ the factorization (51) holds while for ${\mathbf{q}}\neq 0$
the factorization is always violated due to the exponential factor. This
observation, however, does not invalidate our studies presented in this paper
since the form of the double distrtibutions at ${\mathbf{q}}=0$ is essential
for the generalization to the case with ${\mathbf{q}}\neq 0$. For more details
we refer to Diehl:2014vaa ; Rinaldi:2016jvu ; Diehl:2020xyg .

##  VIII Summary

We analyzed the role of the momentum sum rule for the small $x$ factorization
of the double distributions into a product of single distributions in the pure
gluonic case. In general, the momentum sum rule can be treated as consistency
condition for the QCD evolution equations for single and double parton
distributions. Therefore, it must be also fulfilled by initial conditions for
these equations.

We found that when the momentum sum rule is not fulfilled this can lead to
strong violation of the small $x$ factorization at all scales. On the other
hand, when the momentum sum rule is imposed on the initial conditions, then
small $x$ factorization is restored to a good approximation after the
evolution to high scale. This was observed even when there was no
factorization at the initial scale. We stress that for this to happen, the
DPDF evolution equations must include the splitting term. We also note that
the small $x$ factorization observed at high scales is approximate with
accuracy of a few percent. This indicates residual correlations even at
smallest values of $x$ and large values of $Q^{2}$.

We illustrated the impact of the momentum sum rule on the issue of
factorization on several numerical examples in the pure gluonic case and
provide an analytic understanding of the obtained results using the Mellin
moment space. We also discussed the current status of the initial condition
specifications for the QCD evolution, which are motivated by a good knowledge
of the single PDFs. We conclude that the Gaunt-Stirling prescription for the
initial double parton distributions leads to a phenomenologically acceptable
prescription for the small $x$ behavior of the double parton distributions for
the transverse momentum ${\mathbf{q}}=0$. However, it should be kept in mind
that the small $x$ factorization can be violated for ${\mathbf{q}}\neq 0$ due
to correlations in the ${\mathbf{b}}-$space.

##  IX Acknowledgments

This work was supported by the Polish Narodowe Centrum Nauki Grant No.
2019/33/B/ST2/02588 and by the U.S. Department of Energy Grant No. DE-
SC-0002145. We thank Michal Deak for the participation in the early stages of
this analysis.

##  X Appendix A: Momentum sum rule conservation

We will prove that the momentum sum rule (7) is conserved by the evolution
equation (II) by showing that the single PDFs obtained from Eq. (7),

$D_{f_{2}}(x_{2},t)=\frac{1}{1-x_{2}}\sum_{f_{1}}\int_{0}^{1-x_{2}}dx_{1}\,x_{1}D_{f_{1}f_{2}}(x_{1},x_{2},t)\;,$
(53)

obey the DGLAP equation (13). For this purpose we write the LO DGLAP equation
in the most general form

$\partial_{t}D_{f}(x,t)=\sum_{{f^{\prime}}}\int_{0}^{1}du\,{\cal{P}}_{f{f^{\prime}}}(x,u)\,D_{{f^{\prime}}}(u,t)\;,$
(54)

where $f$ indices denote quark/antiquark flavors and gluon. The splitting
kernel has the following general form

${\cal{P}}_{f{f^{\prime}}}(x,u)=P_{f{f^{\prime}}}^{R}(x,u)-\delta(u-x)\,\delta_{f{f^{\prime}}}\,P_{f}^{V}(x)\;,$
(55)

where $R$ and $V$ denote the real and virtual emission kernels, respectively.
For the real emission kernel we have

$P_{f{f^{\prime}}}^{R}(x,u)=\frac{1}{u}P_{f{f^{\prime}}}\\!\\!\left(\frac{x}{u}\right)\theta(u-x)\;,$
(56)

where $P_{f{f^{\prime}}}$ are the LO real emission Altarelli-Parisi splitting
functions. The virtual kernel is determined by assuming that the momentum sum
rule for single PDFs is valid for any $t$,

$\sum_{f}\int_{0}^{1}dx\,xD_{f}(x,t)={\rm const}\;,$ (57)

which holds true during the DGLAP evolution if

$\sum_{f}\int_{0}^{1}dx\,x\,{\cal{P}}_{f{f^{\prime}}}(x,u)=0\;,$ (58)

for any $f^{\prime}$ and $u\in[0,1]$. Substituting (55) into the above we find

$u\,P_{{f^{\prime}}}^{V}(u)=\sum_{f}\int_{0}^{1}dx\,x\,P_{f{f^{\prime}}}^{R}(x,u)\;.$
(59)

In the introduced notation, the evolution equation (II) for DPDFs is given by

$\displaystyle\partial_{t}D_{f_{1}f_{2}}(x_{1},x_{2},t)$
$\displaystyle=\sum_{{f^{\prime}}}\int_{0}^{1-x_{2}}{du}\,{\cal{P}}_{f_{1}{f^{\prime}}}({x_{1}},{u})\,D_{{f^{\prime}}f_{2}}(u,x_{2},t)$
$\displaystyle+\sum_{{f^{\prime}}}\int_{0}^{1-x_{1}}{du}\,{\cal{P}}_{f_{2}{f^{\prime}}}({x_{2}},{u})\,D_{f_{1}{f^{\prime}}}(x_{1},u,t)$
$\displaystyle+\sum_{{f^{\prime}}}\,P_{{f^{\prime}}\to
f_{1}f_{2}}(x_{1},x_{2})\,D_{{f^{\prime}}}(x_{1}+x_{2},t)\;,$ (60)

where the splitting kernel in the third term reads

$P_{{f^{\prime}}\to
f_{1}f_{2}}(x_{1},x_{2})=P^{R}_{f_{1}{f^{\prime}}}(x_{1},x_{1}+x_{2})=P^{R}_{f_{2}{f^{\prime}}}(x_{2},x_{1}+x_{2})\;.$
(61)

The upper integration limits in Eq. (X) can be set to $1$ from the condition

$D_{f_{1}f_{2}}(x_{1},x_{2},t)=0~{}~{}~{}~{}~{}~{}{\rm
for}~{}~{}~{}~{}~{}~{}x_{1}+x_{2}>1\;.$ (62)

Under this condition the upper limit in Eq. (53) can also be set to $1$. Thus,
after differentiating this equation,

$\partial_{t}D_{f_{2}}(x_{2},t)=\frac{1}{1-x_{2}}\sum_{f_{1}}\int_{0}^{1}dx_{1}\,x_{1}\,\partial_{t}D_{f_{1}f_{2}}(x_{1},x_{2},t)\;,$
(63)

and using (X) on the r.h.s., we obtain

$\displaystyle\partial_{t}D_{f_{2}}(x_{2},t)=\sum_{f_{1}}\int_{0}^{1}dx_{1}\frac{x_{1}}{1-x_{2}}$
$\displaystyle\bigg{\\{}\sum_{{f^{\prime}}}\int_{0}^{1}{du}\,{\cal{P}}_{f_{1}{f^{\prime}}}({x_{1}},{u})\,D_{{f^{\prime}}f_{2}}(u,x_{2},t)$
$\displaystyle+\sum_{{f^{\prime}}}\int_{0}^{1}{du}\,{\cal{P}}_{f_{2}{f^{\prime}}}({x_{2}},{u})\,D_{f_{1}{f^{\prime}}}(x_{1},u,t)$
$\displaystyle+\sum_{{f^{\prime}}}P_{{f^{\prime}}\to
f_{1}f_{2}}(x_{1},x_{2})\,D_{{f^{\prime}}}(x_{1}+x_{2},t)\bigg{\\}}\;.$ (64)

The first term on the rhs vanishes due to condition (58). Thus, we find after
changing the summation and integration order that

$\displaystyle\partial_{t}D_{f_{2}}(x_{2},t)$
$\displaystyle=\sum_{f^{\prime}}\int_{0}^{1}\frac{du}{1-x_{2}}\,{\cal{P}}_{f_{2}{f^{\prime}}}({x_{2}},{u})\bigg{\\{}\sum_{f_{1}}\int_{0}^{1}dx_{1}\,x_{1}D_{f_{1}{f^{\prime}}}(x_{1},u,t)\bigg{\\}}$
$\displaystyle+\sum_{f^{\prime}}\int_{0}^{1}dx_{1}\,\frac{x_{1}}{1-x_{2}}\,P^{R}_{f_{2}{f^{\prime}}}(x_{2},x_{1}+x_{2})\,D_{{f^{\prime}}}(x_{1}+x_{2},t)\;,$
(65)

where we used (61) in the last equation to write the second term in Eq. (X).
The sum over $f_{1}$ in this term disappears since for a given
$({f^{\prime}},f_{2})$ there is only one $f_{1}=f_{1}({f^{\prime}},f_{2})$ in
the sum. Applying (53) in the first term and changing the variable $x_{1}\to
u=x_{1}+x_{2}$ in the second one, we obtain

$\partial_{t}D_{f_{2}}(x_{2},t)=\sum_{f^{\prime}}\int_{0}^{1}du\,\frac{1-u}{1-x_{2}}\,{\cal{P}}_{f_{2}{f^{\prime}}}({x_{2}},{u})\,D_{f^{\prime}}(u,t)+\sum_{f^{\prime}}\int_{x_{2}}^{1+x_{2}}du\,\frac{u-x_{2}}{1-x_{2}}\,P^{R}_{f_{2}{f^{\prime}}}(x_{2},u)\,D_{{f^{\prime}}}(u,t)\;.$
(66)

Because of the property (56) and $D_{f}(x,t)=0$ for $x>1$, the integration
range in the second integral may be shifted to $[0,1]$. Therefore, after
taking into account the form (55) of the kernel
${\cal{P}}_{f_{2}{f^{\prime}}}({x_{2}},{u})$ in the first integral, we find

$\partial_{t}D_{f_{2}}(x_{2},t)=\sum_{f^{\prime}}\int_{0}^{1}du\bigg{\\{}P^{R}_{f_{2}{f^{\prime}}}({x_{2}},{u})-\delta(u-x_{2})\,\delta_{f_{2}{f^{\prime}}}\,P^{V}_{f_{2}{f^{\prime}}}({x_{2}})\bigg{\\}}D_{{f^{\prime}}}(u,t)\;,$
(67)

which is the DGLAP equation (54). Notice the crucial role of the splitting
term in the evolution equation (X) to arrive at this conclusion.

##  XI Appendix B: Asymptotic solution in Mellin moment space

The solution to the evolution equation (III) for the double gluon distribution
in the Mellin space is given by Golec-Biernat:2014nsa

$\displaystyle\tilde{D}_{gg}(n_{1},n_{2},t)$
$\displaystyle={\rm{e}}^{\gamma(n_{1})t}\,{\rm{e}}^{\gamma(n_{2})t}\,\tilde{D}_{gg}(n_{1},n_{2},0)+\int_{0}^{t}dt^{\prime}{\rm{e}}^{\gamma(n_{1})(t-t^{\prime})}\,{\rm{e}}^{\gamma(n_{2})(t-t^{\prime})}\,\tilde{\gamma}(n_{1},n_{2})\,\tilde{D}_{g}(n_{1}+n_{2}-1,t^{\prime})\;.$
(68)

The first term in the sum is the general solution to the homogeneous equation
(III) without the splitting term while the second term is a particular
solution to the nonhomogeneous equation with the splitting term.

The homogeneous solution term in Eq. (68) contains a product of exponentials
which generate two independent DGLAP evolutions in the double gluon
distribution since the solution to the DGLAP equation reads as

$\tilde{D}_{g}(n,t)={\rm{e}}^{\gamma(n)t}\,D_{g}(n,0)\;.$ (69)

In the nonhomogeneous solution term in Eq. (68), the gluon splitting $g\to gg$
at ”time” $t^{\prime}$ is followed by two independent DGLAP evolutions with
the exponential factors up to the final $t$. Since the splitting point
$t^{\prime}$ might occur between $0$ and $t$ we need to integrate over the
whole range of such possibilities.

Let us concentrate on the nonhomogeneous solution in Eq. (68), denoted from
now on by $I$. Using (69), we obtain

$\displaystyle I$
$\displaystyle=\int_{0}^{t}dt^{\prime}{\rm{e}}^{\gamma(n_{1})(t-t^{\prime})+\gamma(n_{2})(t-t^{\prime})+\gamma(n_{1}+n_{2}-1)t^{\prime}}\,\tilde{\gamma}(n_{1},n_{2})\,\tilde{D}_{g}(n_{1}+n_{2}-1,0)=$
$\displaystyle={\rm{e}}^{\gamma(n_{1})t+\gamma(n_{2})t}\,\tilde{\gamma}(n_{1},n_{2})\left(\frac{{\rm{e}}^{(\gamma(n_{1}+n_{2}-1)-\gamma(n_{1})-\gamma(n_{2}))t}-1}{\gamma(n_{1}+n_{2}-1)-\gamma(n_{1})-\gamma(n_{2})}\right)\tilde{D}_{g}(n_{1}+n_{2}-1,0)\;.$
(70)

We will make the approximations that reflect the small $x$ limit, in which the
factorization is observed in the numerical analysis. To this aim, we observe
that the splitting function in this limit becomes: ${{P}}_{gg}(z)\simeq
2C_{A}/z$, which leads to the following approximation to the anomalous
dimension

$\gamma(n)=\int_{0}^{1}dz\,z^{n-1}\,{{P}}_{gg}(z)\simeq\int_{0}^{1}dz\,z^{n-1}\,\frac{2C_{A}}{z}=\frac{2C_{A}}{n-1}\;.$
(71)

Thus, the small $x$ limit corresponds to the limit of $(n-1)\to 0$. For the
double distribution we consider the small $x$ limit studied at the plots,
$x_{1}\to 0$ and $x_{2}={\rm fixed}$ and small, which corresponds in the
Mellin moment space to the limit $(n_{1}-1)\to 0$ and $(n_{2}-1)$ finite. In
this limit, the function $\tilde{\gamma}$ given by Eq. (43) is given by

$\tilde{\gamma}(n_{1},n_{2})\simeq
2C_{A}\int_{0}^{1}dx\,x^{n_{1}-2}(1-x)^{n_{2}-1}=2C_{A}\frac{\Gamma(n_{1}-1)\Gamma(n_{2})}{\Gamma(n_{1}+n_{2}-1)}\simeq\frac{2C_{A}}{n_{1}-1}\;,$
(72)

while for the expression in the denominator in Eq. (XI), we obtain

$\gamma(n_{1}+n_{2}-1)-\gamma(n_{1})-\gamma(n_{2})\;\simeq\;-\frac{2C_{A}}{n_{1}-1}\;.$
(73)

So the nonhomogeneous solution (XI) reads as

$I\,\simeq\,e^{\gamma(n_{1})t+\gamma(n_{2})t}\left[1-\exp\left(-\frac{2C_{A}}{n_{1}-1}t\right)\right]\tilde{D}_{g}(n_{1}+n_{2}-1,0)\;.$
(74)

Here one needs to have ${\rm Re}(n_{1})>1$ which is required for the inverse
Mellin transform to make sense. Thus, for sufficiently large $t$, when the
exponent in the square brackets can be neglected, we find

$I\simeq{\rm{e}}^{\gamma(n_{1})t+\gamma(n_{2})t}\tilde{D}_{g}(n_{1}+n_{2}-1,0)\;.$
(75)

Therefore, the full solution (68) can be approximated by

$\tilde{D}_{gg}(n_{1},n_{2},t)\simeq{\rm{e}}^{\gamma(n_{1})t+\gamma(n_{2})t}\left(\tilde{D}_{gg}(n_{1},n_{2},0)+\tilde{D}_{g}(n_{1}+n_{2}-1,0)\right)\;.$
(76)

This needs to be compared with the empirically found result in the considered
small $x$ limit

$\tilde{D}_{gg}(n_{1},n_{2},t)\simeq{\rm{e}}^{\gamma(n_{1})t+\gamma(n_{2})t}\;\tilde{D}_{g}(n_{1},0)\tilde{D}_{g}(n_{2},0)\;,$
(77)

which gives the following relation for the Mellin moments of the initial
condition to be satisfied for the small $x$ factorization to be valid at
sufficiently large $t$,

$\tilde{D}_{gg}(n_{1},n_{2},0)+\tilde{D}_{g}(n_{1}+n_{2}-1,0)=\tilde{D}_{g}(n_{1},0)\,\tilde{D}_{g}(n_{2},0)\;.$
(78)

## References

  * (1) R. Kirschner,  Phys.Lett. B84, 266 (1979).
  * (2) V. P. Shelest, A. M. Snigirev and G. M. Zinovev,  Phys.Lett. B113, 325 (1982).
  * (3) G. Zinovev, A. Snigirev and V. Shelest,  Theor.Math.Phys. 51, 523 (1982).
  * (4) R. K. Ellis, W. Furmanski and R. Petronzio,  Nucl.Phys. B212, 29 (1983).
  * (5) N. Paver and D. Treleani,  Nuovo Cim. A 70, 215 (1982).
  * (6) N. Paver and D. Treleani,  Phys. Lett. B 146, 252 (1984).
  * (7) N. Paver and D. Treleani,  Z. Phys. C 28, 187 (1985).
  * (8) A. Bukhvostov, G. Frolov, L. Lipatov and E. Kuraev,  Nucl.Phys. B258, 601 (1985).
  * (9) A. Kulesza and W. J. Stirling,  Phys.Lett. B475, 168 (2000), [hep-ph/9912232].
  * (10) A. M. Snigirev,  Phys. Rev. D68, 114012 (2003), [hep-ph/0304172].
  * (11) V. L. Korotkikh and A. M. Snigirev,  Phys. Lett. B594, 171 (2004), [hep-ph/0404155].
  * (12) J. R. Gaunt and W. J. Stirling,  JHEP 03, 005 (2010), [0910.4347].
  * (13) B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman,  Phys.Rev. D83, 071501(R) (2011), [1009.2714].
  * (14) F. A. Ceccopieri,  Phys. Lett. B697, 482 (2011), [1011.6586].
  * (15) M. Diehl and A. Schafer,  Phys. Lett. B698, 389 (2011), [1102.3081].
  * (16) J. R. Gaunt and W. J. Stirling,  JHEP 1106, 048 (2011), [1103.1888].
  * (17) M. G. Ryskin and A. M. Snigirev,  Phys.Rev. D83, 114047 (2011), [1103.3495].
  * (18) B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman,  Eur.Phys.J. C72, 1963 (2012), [1106.5533].
  * (19) C. Kom, A. Kulesza and W. Stirling,  Eur.Phys.J. C71, 1802 (2011), [1109.0309].
  * (20) M. Diehl, D. Ostermeier and A. Schafer,  JHEP 1203, 089 (2012), [1111.0910].
  * (21) M. Luszczak, R. Maciula and A. Szczurek,  Phys. Rev. D85, 094034 (2012), [1111.3255].
  * (22) A. V. Manohar and W. J. Waalewijn,  Phys.Rev. D85, 114009 (2012), [1202.3794].
  * (23) M. G. Ryskin and A. M. Snigirev,  Phys.Rev. D86, 014018 (2012), [1203.2330].
  * (24) J. R. Gaunt,  JHEP 1301, 042 (2013), [1207.0480].
  * (25) B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman,  Eur.Phys.J. C74, 2926 (2014), [1306.3763].
  * (26) W. Broniowski and E. Ruiz Arriola,  Few Body Syst. 55, 381 (2014), [1310.8419].
  * (27) M. Diehl, T. Kasemets and S. Keane,  JHEP 05, 118 (2014), [1401.1233].
  * (28) A. van Hameren, R. Maciula and A. Szczurek,  Phys. Rev. D89, 094019 (2014), [1402.6972].
  * (29) K. Golec-Biernat and E. Lewandowska,  Phys.Rev. D90, 014032 (2014), [1402.4079].
  * (30) F. A. Ceccopieri,  Phys. Lett. B 734, 79 (2014), [1403.2167].
  * (31) R. Maciula and A. Szczurek,  Phys. Rev. D90, 014022 (2014), [1403.2595].
  * (32) A. M. Snigirev, N. A. Snigireva and G. M. Zinovjev,  Phys.Rev. D90, 014015 (2014), [1403.6947].
  * (33) K. Golec-Biernat and E. Lewandowska,  Phys.Rev. D90, 094032 (2014), [1407.4038].
  * (34) J. R. Gaunt, R. Maciula and A. Szczurek,  Phys. Rev. D90, 054017 (2014), [1407.5821].
  * (35) L. A. Harland-Lang, V. A. Khoze and M. G. Ryskin,  J. Phys. G42, 055001 (2015), [1409.4785].
  * (36) B. Blok and M. Strikman,  Eur. Phys. J. C74, 3214 (2014), [1410.5064].
  * (37) R. Maciuła and A. Szczurek,  Phys. Lett. B 749, 57 (2015), [1503.08022].
  * (38) K. Golec-Biernat, E. Lewandowska, M. Serino, Z. Snyder and A. M. Stasto,  Phys. Lett. B 750, 559 (2015), [1507.08583].
  * (39) M. Diehl, J. R. Gaunt, D. Ostermeier, P. Plößl and A. Schäfer,  JHEP 01, 076 (2016), [1510.08696].
  * (40) W. Broniowski, E. Ruiz Arriola and K. Golec-Biernat,  Few Body Syst. 57, 405 (2016), [1602.00254].
  * (41) B. Blok and M. Strikman,  Eur. Phys. J. C 76, 694 (2016), [1608.00014].
  * (42) M. Rinaldi, S. Scopetta, M. C. Traini and V. Vento,  JHEP 10, 063 (2016), [1608.02521].
  * (43) K. Golec-Biernat and A. M. Stasto,  Phys. Rev. D 95, 034033 (2017), [1611.02033].
  * (44) F. A. Ceccopieri, M. Rinaldi and S. Scopetta,  Phys. Rev. D 95, 114030 (2017), [1702.05363].
  * (45) M. Diehl, J. R. Gaunt and K. Schönwald,  JHEP 06, 083 (2017), [1702.06486].
  * (46) M. G. A. Buffing, M. Diehl and T. Kasemets,  JHEP 01, 044 (2018), [1708.03528].
  * (47) B. Blok and M. Strikman,  Adv. Ser. Direct. High Energy Phys. 29, 63 (2018), [1709.00334].
  * (48) M. Diehl and J. R. Gaunt,  Adv. Ser. Direct. High Energy Phys. 29, 7 (2018), [1710.04408].
  * (49) E. Elias, K. Golec-Biernat and A. M. Staśto,  JHEP 01, 141 (2018), [1801.00018].
  * (50) M. Rinaldi and F. A. Ceccopieri,  Phys. Rev. D 97, 071501(R) (2018), [1801.04760].
  * (51) M. Rinaldi and F. A. Ceccopieri,  JHEP 09, 097 (2019), [1812.04286].
  * (52) J. R. Gaunt and T. Kasemets,  Adv. High Energy Phys. 2019, 3797394 (2019), [1812.09099].
  * (53) M. Diehl and R. Nagar,  JHEP 04, 124 (2019), [1812.09509].
  * (54) M. Diehl, J. R. Gaunt, P. Plößl and A. Schäfer,  SciPost Phys. 7, 017 (2019), [1902.08019].
  * (55) W. Broniowski and E. R. Arriola,  Phys. Rev. D 101, 014019 (2020), [1910.03707].
  * (56) M. Diehl, J. R. Gaunt, D. M. Lang, P. Plößl and A. Schäfer,  Eur. Phys. J. C 80, 468 (2020), [2001.10428].
  * (57) G. S. Bali et al.,  JHEP 02, 067 (2021), [2006.14826].
  * (58) B. Cabouat and J. R. Gaunt,  JHEP 10, 012 (2020), [2008.01442].
  * (59) M. Diehl, J. R. Gaunt and P. Ploessl,  JHEP 08, 040 (2021), [2105.08425].
  * (60) M. Diehl, J. R. Gaunt, P. Pichini and P. Plößl,  Eur. Phys. J. C 81, 1033 (2021), [2109.14304].
  * (61) Axial Field Spectrometer Collaboration, T. Akesson et al.,  Z.Phys. C34, 163 (1987).
  * (62) CDF Collaboration, F. Abe et al.,  Phys.Rev.Lett. 79, 584 (1997).
  * (63) CDF Collaboration, F. Abe et al.,  Phys.Rev. D56, 3811 (1997).
  * (64) D0 Collaboration, V. Abazov et al.,  Phys.Rev. D81, 052012 (2010), [0912.5104].
  * (65) ATLAS Collaboration, G. Aad et al.,  New J.Phys. 15, 033038 (2013), [1301.6872].
  * (66) ATLAS Collaboration, G. Aad et al.,  JHEP 1404, 172 (2014), [1401.2831].
  * (67) ATLAS, G. Aad et al.,  Eur. Phys. J. C 75, 229 (2015), [1412.6428].
  * (68) ATLAS, M. Aaboud et al.,  JHEP 11, 110 (2016), [1608.01857].
  * (69) ATLAS, M. Aaboud et al.,  Eur. Phys. J. C 77, 76 (2017), [1612.02950].
  * (70) ATLAS, M. Aaboud et al.,  Phys. Lett. B 790, 595 (2019), [1811.11094].
  * (71) CMS Collaboration, S. Chatrchyan et al.,  JHEP 1403, 032 (2014), [1312.5729].
  * (72) CMS, V. Khachatryan et al.,  JHEP 05, 013 (2017), [1610.07095].
  * (73) CMS, A. M. Sirunyan et al.,  Eur. Phys. J. C 80, 41 (2020), [1909.06265].
  * (74) CMS,  2206.02681.
  * (75) LHCb, R. Aaij et al.,  Phys. Lett. B 707, 52 (2012), [1109.0963].
  * (76) LHCb, R. Aaij et al.,  JHEP 06, 141 (2012), [1205.0975],  [Addendum: JHEP 03, 108 (2014)].
  * (77) LHCb, R. Aaij et al.,  JHEP 07, 052 (2016), [1510.05949].
  * (78) LHCb, R. Aaij et al.,  JHEP 06, 047 (2017), [1612.07451],  [Erratum: JHEP 10, 068 (2017)].
  * (79) J. R. Gaunt,  Double parton scattering in proton-proton collisions,  PhD thesis, University of Cambridge, 2012,  https://www.repository.cam.ac.uk/handle/1810/243945.
  * (80) A. Martin, W. Stirling, R. Thorne and G. Watt,  Eur.Phys.J. C63, 189 (2009), [0901.0002].