# Extraction of the Muon Signals Recorded with the Surface Detector of the
Pierre Auger Observatory Using Recurrent Neural Networks

The Pierre Auger Collaboration

###### Abstract

The Pierre Auger Observatory, at present the largest cosmic-ray observatory
ever built, is instrumented with a ground array of 1600 water-Cherenkov
detectors, known as the Surface Detector (SD). The SD samples the secondary
particle content (mostly photons, electrons, positrons and muons) of extensive
air showers initiated by cosmic rays with energies ranging from $10^{17}~{}$eV
up to more than $10^{20}~{}$eV. Measuring the independent contribution of the
muon component to the total registered signal is crucial to enhance the
capability of the Observatory to estimate the mass of the cosmic rays on an
event-by-event basis. However, with the current design of the SD, it is
difficult to straightforwardly separate the contributions of muons to the SD
time traces from those of photons, electrons and positrons. In this paper, we
present a method aimed at extracting the muon component of the time traces
registered with each individual detector of the SD using Recurrent Neural
Networks. We derive the performances of the method by training the neural
network on simulations, in which the muon and the electromagnetic components
of the traces are known. We conclude this work showing the performance of this
method on experimental data of the Pierre Auger Observatory. We find that our
predictions agree with the parameterizations obtained by the AGASA
collaboration to describe the lateral distributions of the electromagnetic and
muonic components of extensive air showers.

##  1 Introduction

The existence of ultra-high energy cosmic rays (UHECRs) is one of the most
intriguing phenomena in astroparticle physics. UHECRs are atomic nuclei with
energies above $10^{18}$ eV and their flux falls off quickly with energy [1,
2]. The most energetic particles, above $10^{20}$ eV, are detected with a
frequency smaller than one per square kilometer per five centuries. The Pierre
Auger Observatory was designed and built to unveil the origin of these
particles [3]. It is the largest cosmic-ray detector built so far and
comprises an area of 3000 km2. The detection of air showers is performed using
two complementary techniques. On the one hand, the fluorescence light produced
as the shower develops along the atmosphere is measured at four sites of the
Fluorescence Detectors (FD). On the other hand, the particles that reach the
ground are sampled with the SD, an array of 1600 water-Cherenkov detectors
(WCDs), separated 1500 m from each other.

Figure 1:  Left: A schematic view of a surface detector station in the field,
showing its main components. Right: Signal in time measured by a single
station of the SD for a simulated event. The total signal corresponds to the
signal that can be measured by the SD. The signal from the electromagnetic and
muon components are shown independently, as this information is available in
simulations. The total signal is the sum of the muon component and the
electromagnetic component. Our goal is to obtain the muon component using the
total signal as an input.

This work focuses on the data collected with the SD, which operates with a
duty cycle close to 100% and shows full detection efficiency for air showers
above 3.16 EeV ($10^{18.5}$ eV). Charged particles traverse the 12 000 litres
of ultra-pure water contained in each station of the SD and the emitted
Cherenkov photons are collected by three nine-inch photomultiplier tubes
(PMTs); see the left panel of Figure 1. The photomultiplier output is
digitised at 40 MHz (25 ns bins) using 10 bit Flash Analog-Digital Converters
(FADCs). Each signal trace has a length of 768 samples, corresponding to a
time window of 19.2 $\upmu$s. The resulting signal is proportional to the sum
of the Cherenkov radiation in water produced by electromagnetic particles
(photons, electrons and positrons) and muons, see the right panel of Figure 1.
Muons produce some characteristic features in the signal. As highly
penetrating particles, they usually arrive at the stations at earlier times
than the electromagnetic particles that undergo copious multiple scattering,
see Fig. 4 in ref. [4]. The signal that muons produce is spiky, in contrast to
the signal from electromagnetic particles, which is more continuous and spread
in time. Since muons suffer less energy losses in the atmosphere than
electromagnetic particles [5], they dominate the signals of the stations
located far away from the core of the shower.

Based on the Heitler-Matthews extensive air-shower model and assuming the
superposition principle [6], the number of muons $N_{\mu}^{A}$ in an extensive
air shower produced by a nucleus with mass number $A$ can be related to the
number of muons produced in a shower initiated by a proton with the same
energy, $N_{\mu}^{\text{p}}$, through

$N_{\mu}^{A}=N_{\mu}^{\text{p}}A^{1-\beta}.$ (1.1)

where $1-\beta\simeq 0.1$. Therefore, the determination of the muon component
in the air shower is crucial to infer, for each event, the mass of the primary
particle, which is a key ingredient in the searches conducted to pinpoint the
sources of UHECRs. Mass composition inferences are customarily carried out
comparing observables measured in data and simulations.

With the current design of the SD, the separation of the muon and
electromagnetic components can be done for events with large zenith angles [7]
or for distances far off from the shower core [8]. Hence, for a majority of
the recorded data, this separation cannot be performed in a straightforward
and efficient way. In this work, starting from the total time trace, we
estimate the muon signal for each individual WCD as a function of time using
machine learning techniques. The applicability of the method is not restricted
to a limited energy and/or zenith angle range but can be applied to the full
data sample collected so far. However, for practical purposes, in this work we
limit the discussion to the data set referred to as “vertical” events
($\theta<60^{\circ}$, where $\theta$ is the zenith angle) which are those for
which it is harder to estimate the muon component with more straightforward
techniques. As a novelty, we estimate the risetime [9] of the muon signal.

The branch of machine learning research studies algorithms and techniques that
rely on inferring patterns from data. The scientist is only in charge of
developing a suitable architecture, instead of designing every detail of a
model. Numerous techniques were proposed many years ago but several recent
breakthroughs in computational hardware have allowed for these techniques to
be computationally affordable. Nowadays, machine learning is being used
extensively in physics (for a recent review see e.g. ref. [10]). It is
especially well suited for experiments in particle physics [11] and in
astroparticle physics. The large amount of events collected by state-of-the-
art experiments in both fields allows for the exploitation of machine learning
techniques that require vast amounts of data to build and train a model.

An example of a previous work devoted to the extraction of the muon signal
using machine learning can be found in ref. [12]. In this work, we focus our
interest on the study of time series. Hence, from the large set of machine
learning algorithms, we have chosen deep learning since it is the one that has
lately experienced substantial progress when applied to temporal series and
sequences. We go one step beyond to what is discussed in ref. [12], where only
the integral of the muon trace was obtained. For the first time, we estimate
the muon contribution to the signal recorded in each time bin.

This article is structured as follows. In Section 2, we give an overview of
the methodology, the architecture of the neural network and the training
method. In Section 3, we show the results of the neural network. We conclude
by discussing the application of the neural network to data and assess its
performance scrutinizing standard properties of extensive air showers.

##  2 The method

Our approach is based on Recurrent Neural Networks (RNNs). RNNs are specially
well-suited for time series due to their memory mechanism. As the computing
process develops, for each step of the temporal series, the RNN stores
information from a preceding time slot that can be used at a later stage.
Nowadays, RNNs are used successfully in different fields such as natural
language processing or machine translation, see for example [13]. There are
several kinds of RNNs and, among them, we have chosen one of the most common,
known as Long Short Term Memory (LSTM) [14, 15].

###  2.1 The input

The traces recorded by the water-Cherenkov detectors of the SD are the input
to the neural network. The electronics is used to sample the signals in bins
of 25 ns. Once they are calibrated, the signals are expressed in Vertical
Equivalent Muons (VEMs), which correspond to the most-likely signal deposited
by a muon that traverses the center of the detector top-down. The traces used
in this study are obtained after averaging the traces of every functioning PMT
in a WCD.

We use only the first 200 bins (5 $\upmu$s) of each trace. Traces that are
shorter than those 200 bins are padded with zeros at the end. This choice is
made because the most relevant information is encapsulated in the first 200
bins of each trace, while choosing too many bins increases the computational
time and memory needed to train the neural network. From simulations we know
that for $E<10^{19}$ eV, around 90% of the stations have the complete muon
signal in the first 200 bins and the remaining stations have more than 99% of
the muon signal in those 200 bins. For $E>10^{19}$ eV, around 70% of the
stations have the complete muon signal in the first 200 bins and the remaining
stations have around 99% of the muon signal contained in those 200 bins. The
fraction of the muon signal outside the first 200 bins is then negligible and
can be ignored.

Figure 2: Schematic diagram of the geometry of an event, modified from ref.
[16]. The angle $\theta$ is the angle between the shower axis and the zenith,
the distance $r$ is the distance between the shower axis and the perpendicular
projection of the station in the shower plane.

The trace information is not enough to accurately determine the muon component
of extensive air showers. We found that the amount of atmosphere that the
particles traverse plays an important role, since the electromagnetic and muon
components are attenuated differently in the atmosphere. Particles from the
electromagnetic component interact and scatter; the cascade continues down to
energies below 1 MeV until electrons slow down through ionization without
further radiating. In contrast, muons are very penetrating particles that
traverse the atmosphere practically unaffected. This is why they also arrive
earlier than electromagnetic particles. Muons are typically minimum ionizing
particles. Consequently, as the attenuation of electrons, positrons and
photons increases, the traces richer in muons become spikier and shorter in
time.

This is the reason behind our choice of using as input two more variables: the
secant of the reconstructed zenith angle, $\sec\theta$, and the distance to
the core on the plane perpendicular to the shower (the shower plane), $r$, see
Figure 2. $\theta$ is the reconstructed angle between the zenith and the
trajectory of the primary cosmic ray. It results from a geometrical
reconstruction carried out using the arrival time of particles at the stations
and their respective locations [17]. From this reconstruction, the position of
the core of the shower at the ground is estimated. Both variables take into
account the amount of atmosphere crossed by particles. For $\theta\lesssim
80^{\circ}$, the amount of traversed atmosphere is proportional to
$\sec\theta$ and as $r$ increases, the distance travelled by particles through
the atmosphere gets larger.

Other variables, such as the energy of the primary cosmic ray, do not help to
improve the performance of the neural network when they are included. By not
including it, we also avoid the need to match the energy scale of data with
that of simulations, which proves to be a difficult task. The application of
this method to data is then straightforward since the zenith angle and
distance to the core depend only on the geometrical reconstruction of the
event.

###  2.2 Recurrent Neural Network (RNN) architecture

Figure 3: Schematic drawing of the input, architecture and output of the
neural network. See the text for details.

The neural network architecture depicted in Figure 3 is set up as follows. The
input is a vector of 202 components: distance to the core $r$, secant of the
zenith angle $\sec\theta$ and the 200 values of each trace
$S_{1},S_{2},\dots,S_{200}$, where $S_{i}$ is the value of the signal measured
at the time bin $t_{i}$. At the start, this input is split into two sets. One
of them is the set of the time-independent variables, i.e. $r$ and
$\sec\theta$. These variables are fed into a set of fully-connected layers
that will compute the initial hidden state and cell state for the first layer
of LSTMs. The fully-connected layers transform the inital vector of two values
($r$, $\sec\theta$) into a vector of dimension 32 and then to a vector of
dimension 70. The outputs of the fully-connected layers together with the 200
time values of each trace form the input to the LSTMs. Starting from a single
sequence of 200 numbers (the signal $S_{i}$), the LSTMs produce 70, 32 and 32
sequences of 200 numbers, and the last one of these sequences is fed to a
final fully-connected layer. The activation function used for the fully
connected layer is $f(x)=\max(0,x)$ and $\tanh$ for the LSTM layers. The model
has a total of 87 212 trainable or free parameters.

The set of fully-connected layers computes the vector of initial parameters
that encodes information about the amount of atmosphere crossed by the
particles, which in turn changes the shape of the traces. Without this block,
there is insufficient information for the neural network to distinguish
between traces with high or low fractions of muons and will be biased: for
example, it will overestimate the muon component for nearly vertical showers
and underestimate it for more inclined showers. The block of LSTMs is
responsible for computing the traces and using the temporal information of the
input.

###  2.3 Data selection and training

The neural network is trained with a library of simulated showers, with the
same number of showers initiated by protons, helium, oxygen and iron nuclei.
The energy range covered by simulations spans from $10^{18.5}$ eV to
$10^{20.2}$ eV for zenith angles up to $60^{\circ}$. Each shower is simulated
with the CORSIKA software [18] using the EPOS-LHC hadronic interaction model
[19] and then reconstructed using the $\overline{\textrm{Off}}$
$\underline{\textrm{line}}$ software of the Pierre Auger Collaboration [20].
Among those simulated events, we select those that fulfill the quality
condition that the station measuring the highest signal has to be surrounded
by six operating neighbours, thus avoiding events that fall at the edges of
the array. The method is applied to traces whose signals do not show any sign
of saturation caused by an overflow of the read-out electronics or a loss of
the linear behaviour of the PMT. Finally, only traces for which the integral
is more than 5 VEM are included.

A total of around 450 000 events were available. The events were sampled
randomly and assigned to the training, validation or test data sets using a
uniform distribution in energy and $\sec\theta$ for the validation and test
sets; the remaining events were assigned to the training set. The training
data set does not require special care regarding the energy and zenith angle
distributions since the dependence on the energy is mild and the zenith angle
is given as input. The whole event sample was split as follows: around 390 000
in the training set, 22 000 in the validation set and 34 000 in the test set.

Before training, both $r$ and $\sec\theta$ values are scaled to be between 0
and 1 and all the traces are scaled individually to be between 0 and 1. The
true (simulated) muon trace is scaled with the same factor used for the total
trace. The output of the neural network is also between 0 and 1 so as to use
the same factor to rescale back the predicted trace. The loss function that is
minimized in the process of training is the mean squared error, defined for
the trace of a single WCD as

$L=\frac{1}{200}\sum_{i=1}^{200}\left(\widehat{S\rule{0.0pt}{7.0pt}^{\mu}}_{i}-{S\rule{0.0pt}{7.0pt}^{\mu}}_{i}\right)^{2}.$
(2.1)

It corresponds to the average of the squares of the differences between the
true muon trace ${S\rule{0.0pt}{7.0pt}^{\mu}}_{i}$ and the predicted muon
trace111In this work, we use a hat $\ \widehat{}\ $ for all the quantities
that are either predicted by the neural network or computed from predictions
of the neural network. For example, the integral of the muon signal is
$S^{\mu}$, while the integral of the predicted muon signal is
$\widehat{S^{\mu}}$. $\widehat{S\rule{0.0pt}{7.0pt}^{\mu}}_{i}$, for each time
bin $i$. The neural network is trained in batches and for each batch, the
value of $L$ is computed for each trace and averaged over all the samples in
the batch.

The training was done with the optimizer ADAM [21] with a fixed learning rate
of $10^{-4}$ and the default values for the rest of the parameters; see ref.
[21]. Using a batch size of 516 and 150 complete iterations over the whole
training set or epochs on an Nvidia 2080 Ti GPU, the training takes around 8
hours. The loss as a function of the epoch is shown in Figure 4 for both the
training and validation sets. We can see that the curve for the validation set
decreases as the epoch increases. The curve for validation is below the one
for training because, as explained before, a uniform distribution has been
used for the validation data set and there are more events for which the
performance is worse (lower zenith angles) in the training set. To further
assess the goodness of our approach, the difference between the integral of
the predicted and true muon traces is computed after each epoch. In the right
panel of Figure 4, we show the mean and the standard deviation of the
distribution of the difference. We use this measure because it has a
straightforward physical interpretation. The mean controls the bias: i.e., if
the mean is above zero, the neural network is consistently overestimating and
vice versa. The standard deviation tells us how well we are predicting the
muon signal, although this depends on many factors such as $\theta$, as we
will see in Section 3. All the pipeline was implemented in Python 3.8 using
Numpy [22] and Pandas [23] for the data treatment and Pytorch 1.5.0 [24] for
the construction and training of the neural network.

Figure 4: Left: Loss as a function of the epoch, see Eq. 2.1. Right: Mean
value and standard deviation of the difference between the integral of the
true muon signal and the predicted muon signal for the validation set.

Several tests were done to improve the performance of the method. Instead of
predicting the average muon trace, we aimed at predicting the muon signal for
each single PMT. These three independent results were averaged and compared to
the main method (i.e., feeding the network with the average of the signals of
the active PMTs). We found that even though the number of traces nearly
tripled, the performance did not improve, with an increase of 0.2 VEM in the
resolution or, in other words, around 1% with respect to the total signal (see
the bottom right panel of Figure 7). We also tried other changes in the neural
network such as increasing the number of layers or using a decaying learning
rate, with no significant improvements in performance.

##  3 Results

We show examples of the predictions of the neural network in Figure 5. As a
general comment, we observe that, qualitatively, the prediction follows the
shape and peaks of the total signal. The network has learnt to reproduce the
main features of the muon trace: its spiky shape and the fact that most muons
arrive earlier. A thorough discussion of the results and their dependence on
several variables are given in the next sections.

Figure 5: Examples of predicted muon traces for two simulated events with
EPOS-LHC, for a electromagnetic-dominated signal (left) and muon-dominated
signal (right). The prediction (black line) agrees well with the shape of the
simulated muon trace (orange line) for a majority of the time bins. The blue
thicker line corresponds to the total trace, the one measured by a WCD.

###  3.1 Integrals of the trace

Figure 6: Left: Distribution of the difference between the integral of the
predicted muon signal $\widehat{S^{\mu}}$ and the integral of the true muon
signal $S^{\mu}$. Right: Distribution of $\widehat{S^{\mu}}$ and $S^{\mu}$ for
all the stations in the test set.

One way to assess the performance of the method is to compute the integral of
the predicted muon trace,
$\widehat{S^{\mu}}=\sum_{i=1}^{200}\widehat{S\rule{0.0pt}{7.0pt}^{\mu}}_{i}$
and compare it to the integral of the true muon trace
$S^{\mu}=\sum_{i=1}^{200}{S\rule{0.0pt}{7.0pt}^{\mu}}_{i}$. The integral of
the muon trace is an interesting physical observable, since it relates to the
total number of muons that reach the ground. In the left panel of Figure 6, we
show a distribution with the difference between $\widehat{S^{\mu}}$ and
$S^{\mu}$ for a particular bin of energy and zenith angle. The difference is
compatible with zero and does not show a strong dependence with the value of
the true muon signal. In the right panel of Figure 6, we show the distribution
of $\widehat{S^{\mu}}$ and $S^{\mu}$ for all the WCDs in the test set for
showers initiated by a proton primary. The two distributions have similar
shapes.

Figure 7: Left: Mean value (top) and standard deviation (bottom) of the
difference between the predicted and values of the true muon signal as a
function of energy. Right: Relative bias (top) and resolution (bottom) for the
determination of the muon fraction as a function of the energy. Figure 8:
Left: Mean value (top) and standard deviation (bottom) of the difference
between the predicted and values of the true muon signal as a function of
$\sec\theta$. Right: Relative bias (top) and resolution (bottom) for the
determination of the muon fraction as a function of the zenith angle.

In Figure 7 and Figure 8, the performances of the method are depicted as a
function of the energy and $\sec\theta$, respectively. In the left panels, the
mean value and standard deviation of the difference between the predicted and
true values are shown for the different primaries. The mean values are close
to zero and the standard deviation is less than 2 VEMs in most bins. These
performances readily translate into the biases and resolutions for the muon
fraction with respect to the total recorded signal shown in the right panels.
The biases are within 2% for all the energies and angles explored here,
regardless of the primaries. The resolutions are of the order of 11% at
$10^{18.5}~{}$eV, improving with higher energies. The dependence on
$\sec\theta$ reflects the attenuation of the overall signal at large zenith
angles.

Figure 9: Integral of the predicted muon trace as a function of the integral
of the true muon trace. The black line corresponds to a perfect prediction and
$\rho$ is the Pearson correlation coefficient.

The performance of the method can be also evaluated as a function of the
simulated muon signal ($S^{\mu}$). In Figure 9, we show the density plot of
the predictions and the true values. The predictions are highly correlated
with a Pearson correlation coefficient that is 0.99 for all the primaries. The
mean and standard deviation of the distribution of $\widehat{S^{\mu}}-S^{\mu}$
are shown in Figure 10 as a function of $S^{\mu}$. The mean is close to zero
and rarely exceeds 2 VEM. These results are observed to be zenith-angle
dependent. The standard deviation increases as $S^{\mu}$ increases for more
vertical events (left panel), while it remains constant at large angles (right
panel) as a result of the preponderance of the muon signal.

Figure 10: Mean and standard deviation of the difference between the integral
of the predicted muon signal $\widehat{S^{\mu}}$ and the integral of the true
muon signal $S^{\mu}$ for all the stations from events with the energies and
zenith angles specified in the boxes.

For $S^{\mu}>5$ VEM, the resolution in the determination of the muon signal
for each individual WCD goes from about 20% for the case of vertical events to
10% for events with $\sec\theta\simeq 2.0$

###  3.2 Risetime of the muon signal

The risetime of the recorded signals is another observable worth analysing,
since it has been used successfully by the Pierre Auger Collaboration to
extract valuable information about extensive air showers [16, 9]. The risetime
$t_{1/2}$ is defined as the difference between two times, $t_{10}$ and
$t_{50}$. $t_{10}$ is the time corresponding to the first time bin where the
integral of the signal reaches 10 % of its total value, while $t_{50}$ is
obtained when the integral reaches 50 % of the total. The risetime gives
information about the shape of the trace: traces where the signal is
concentrated in a few bins will have a shorter risetime, while traces where
the signal is spread over time will have larger risetimes. In particular, the
muon component has a smaller risetime than the electromagnetic component,
since muons arrive earlier and in a shorter window of time [4].

In Figure 11, we compare the risetime of the predicted muon trace
$\widehat{t\rule{0.0pt}{7.0pt}^{\mu}_{1/2}}$ with the risetime of the
simulated muon trace $t\rule{0.0pt}{7.0pt}^{\mu}_{1/2}$. The standard
deviation is less than 100 ns (4 time bins) for most values of the true
risetime. This is a small time compared to the risetime of the total signal,
which has a mean of 500 ns and 150 ns for the left and rights bins of Figure
11, respectively. Values with $t\rule{0.0pt}{7.0pt}^{\mu}_{1/2}<100$ ns
correspond to around 1.5% of the samples in the left panel of Figure 11, since
for vertical showers it is very rare to have the muon signal concentrated in
very few bins, while for inclined showers it is rare that the muon signal is
spread over a large time. The mean value approaches zero and the performance
improves as the zenith angle increases. This means that the network
successfully predicts not only the integral of the muon trace but also the
shape of the muon trace.

Figure 11: Mean and standard deviation of the difference between the risetime
of the predicted muon signal $\widehat{t\rule{0.0pt}{7.0pt}^{\mu}_{1/2}}$ and
the risetime of the true muon signal $t\rule{0.0pt}{7.0pt}^{\mu}_{1/2}$ for
all the stations from events with the energies and zenith angles specified in
the boxes.

###  3.3 Hadronic interaction model

We have chosen to train our neural network on simulations with the EPOS-LHC
generator of hadronic interactions and presented some tests in the previous
sections using the same event generator. We now test our method using
simulations done with QGSJetII-04 [25] and Sibyll 2.3 [26] event generators.
In other words, we test the predictions for simulations that are not only
unknown to the neural network but that also have been generated using a
different model of hadronic interactions.

#### QGSJetII-04.

In Figure 12, an example of a trace obtained with simulations produced using
QGSJetII-04 is shown. The prediction follows the shape of the peaks,
predicting quite accurately the muon signal. The difference between true and
predicted muon signals does not show a strong deviance from zero.

Figure 12: Left: Example of a predicted trace for a simulation of a proton
generated air shower done with QGSJetII-04. Right: Distribution of
$\widehat{S^{\mu}}-S^{\mu}$ for all the stations in the bin specified for
simulations using proton and iron nuclei.

#### Sibyll 2.3.

In Figure 13, we show the result of predicting the muon traces for simulations
performed with Sibyll 2.3. By comparing the values of the bias and of the
resolution, we can see that the performance is similar to the case where the
predictions were based on simulations done with QGSJetII-04.

Figure 13: Left: Example of a predicted trace for a simulation done with
Sibyll 2.3 with a proton as primary cosmic ray. Right: Distribution of
$\widehat{S^{\mu}}-S^{\mu}$ for all the stations in the bin specified for
simulations using proton and iron nuclei.

With these results we show that the neural network predictions rely
essentially on the response of the detectors to muons and are thus independent
of the hadronic interaction model used to simulate extensive air showers. For
completeness, we have also carried out the opposite exercise: train the neural
network with QGSJetII-04 and Sibyll 2.3 and predict for the other two models,
respectively. The outcome is in a good agreement with the case where the
neural network learns from events simulated with EPOS-LHC.

Figure 14: Examples of the predicted muon traces for two WCDs that belong to
two different events recorded by the SD.

###  3.4 Comparison to data

In this section we assess how the neural network performs when applied to
experimental data. The sample of events selected have the same cuts as the
simulations used for training. In total, there are 177 000 events registered
from 2004 to 2017. In Figure 14, we show examples of the muon trace predicted
for two typical signals recorded with two independent WCDs. We observe that
the output of the neural network produces features similar to those shown by
the simulated traces of Figure 5, namely: the predicted muon signal is larger
at earlier times since muons arrive earlier; in addition, the predicted muon
trace exhibits a spiky structure.

Figure 15: Left: Lateral distribution of the predicted muon and
electromagnetic components for data. The solid and dashed lines are the
outcome of two independent fits (see text for details). Right: Measured
risetime of the total signal and risetime of the predicted muon signal as a
function of the distance for two different zenith angle ranges.

A straightforward way to verify the robustness of the muon signal measurement
is to study the behaviour of the total muon signal with the distance to the
air-shower axis, i.e. the muon lateral distribution function (LDF). The muon
LDF is illustrated in Figure 15, together with the electromagnetic one. For
each WCD, we can obtain the electromagnetic component simply by subtracting
the predicted muon signal from the total recorded signal. This endows us with
the additional possibility to study the behaviour of the LDF of
electromagnetic signals.

In addition, we compare our data to the parameterizations of the LDF of muon
and electromagnetic particles that best fit the Akeno measurements. The
authors of ref. [27] studied the properties of muons above 1 GeV at different
distances to the shower core for events with $\sec\theta<1.2$. The behaviour
of the Akeno data is well reproduced by the function proposed by Greisen in
ref. [28], Eqs. (4) and (6) in ref. [27]. We fit our data with those
expressions and the parameters that best reproduce Akeno data, just leaving
the overall normalization free. The outcome of the fit is shown as a solid
line in Figure 15 for the energy range $10^{18.6}$ eV to $10^{18.7}$ eV. The
level of agreement between our measurements and the muon lateral distribution
function extracted from Akeno data is remarkable.

In ref. [29], the Akeno data are used to extract a parameterization of the
density of the electromagnetic signal as a function of the distance to the
shower core (Eq. (2.1) in ref. [29]). We use that parameterization to fit our
predicted electromagnetic signal, leaving the overall normalization free and
the parameter referred to as $\eta$ in the original publication. $\eta$ is
left as a free parameter since it is a function of the measured
electromagnetic signal (see Eq. (2.2) in ref. [29]). As for the case of muons,
the agreement with the Akeno parameterization of the electromagnetic signals
is very good (dashed line in Figure 15).

The data sample was also used to carry out a cross-check involving the arrival
times of the particles: we compared the behaviour with the distance to the
core and with the zenith angle of the risetimes measured for the total signal
and the predicted muon signal, see the right panel of Figure 15. As expected,
the muon risetime is consistently smaller since muons arrive earlier. As the
zenith angle increases, the fraction of the muon signal grows, thus the
risetime of the total signal and the muon risetime become more similar. This
behaviour is correctly reproduced by the predictions of the RNN.

##  4 Conclusions

We have shown that over a broad range of energies and zenith angles, a
Recurrent Neural Network is a suitable tool to predict the muon signals that
are part of the time traces measured with each individual WCD of the SD of the
Pierre Auger Observatory. The neural network is trained on simulations, but
the predictions are independent of the model used to simulate hadronic
interactions. The muon signal for each WCD can be measured with a resolution
that decreases from 20% to 10% as the zenith angle increases. Likewise, the
muon fraction with respect to the total signal is estimated with biases that
are within a 2% and a resolution that is lower than 11%.

When applied to data, the behaviour of the extracted muon and electromagnetic
signals agrees well with relevant measurements found in the literature while
other observables, such as the risetime, follow the behaviour dictated by
basic physics principles.

The combination of this algorithm with the future data collected with the
upgraded Observatory [30] will represent a major step forward since arguably,
we will achieve an unprecedented resolution in our capability to make mass
estimates on an event-by-event basis.

## Acknowledgments

The successful installation, commissioning, and operation of the Pierre Auger
Observatory would not have been possible without the strong commitment and
effort from the technical and administrative staff in Malargüe. We are very
grateful to the following agencies and organizations for financial support:

Argentina – Comisión Nacional de Energía Atómica; Agencia Nacional de
Promoción Científica y Tecnológica (ANPCyT); Consejo Nacional de
Investigaciones Científicas y Técnicas (CONICET); Gobierno de la Provincia de
Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in
gratitude for their continuing cooperation over land access; Australia – the
Australian Research Council; Brazil – Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq); Financiadora de Estudos e Projetos (FINEP);
Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ); São Paulo
Research Foundation (FAPESP) Grants No. 2019/10151-2, No. 2010/07359-6 and No.
1999/05404-3; Ministério da Ciência, Tecnologia, Inovações e Comunicações
(MCTIC); Czech Republic – Grant No. MSMT CR LTT18004, LM2015038, LM2018102,
CZ.02.1.01/0.0/0.0/16_013/0001402, CZ.02.1.01/0.0/0.0/18_046/0016010 and
CZ.02.1.01/0.0/0.0/17_049/0008422; France – Centre de Calcul IN2P3/CNRS;
Centre National de la Recherche Scientifique (CNRS); Conseil Régional Ile-de-
France; Département Physique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS);
Département Sciences de l’Univers (SDU-INSU/CNRS); Institut Lagrange de Paris
(ILP) Grant No. LABEX ANR-10-LABX-63 within the Investissements d’Avenir
Programme Grant No. ANR-11-IDEX-0004-02; Germany – Bundesministerium für
Bildung und Forschung (BMBF); Deutsche Forschungsgemeinschaft (DFG);
Finanzministerium Baden-Württemberg; Helmholtz Alliance for Astroparticle
Physics (HAP); Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF);
Ministerium für Innovation, Wissenschaft und Forschung des Landes Nordrhein-
Westfalen; Ministerium für Wissenschaft, Forschung und Kunst des Landes Baden-
Württemberg; Italy – Istituto Nazionale di Fisica Nucleare (INFN); Istituto
Nazionale di Astrofisica (INAF); Ministero dell’Istruzione, dell’Universitá e
della Ricerca (MIUR); CETEMPS Center of Excellence; Ministero degli Affari
Esteri (MAE); México – Consejo Nacional de Ciencia y Tecnología (CONACYT) No.
167733; Universidad Nacional Autónoma de México (UNAM); PAPIIT DGAPA-UNAM; The
Netherlands – Ministry of Education, Culture and Science; Netherlands
Organisation for Scientific Research (NWO); Dutch national e-infrastructure
with the support of SURF Cooperative; Poland -Ministry of Science and Higher
Education, grant No. DIR/WK/2018/11; National Science Centre, Grants No.
2013/08/M/ST9/00322, No. 2016/23/B/ST9/01635 and No. HARMONIA
5–2013/10/M/ST9/00062, UMO-2016/22/M/ST9/00198; Portugal – Portuguese national
funds and FEDER funds within Programa Operacional Factores de Competitividade
through Fundação para a Ciência e a Tecnologia (COMPETE); Romania – Romanian
Ministry of Education and Research, the Program Nucleu within MCI
(PN19150201/16N/2019 and PN19060102) and project PN-
III-P1-1.2-PCCDI-2017-0839/19PCCDI/2018 within PNCDI III; Slovenia – Slovenian
Research Agency, grants P1-0031, P1-0385, I0-0033, N1-0111; Spain – Ministerio
de Economía, Industria y Competitividad (FPA2017-85114-P and
PID2019-104676GB-C32, Xunta de Galicia (ED431C 2017/07), Junta de Andalucía
(SOMM17/6104/UGR, P18-FR-4314) Feder Funds, RENATA Red Nacional Temática de
Astropartículas (FPA2015-68783-REDT) and María de Maeztu Unit of Excellence
(MDM-2016-0692); USA – Department of Energy, Contracts No. DE-AC02-07CH11359,
No. DE-FR02-04ER41300, No. DE-FG02-99ER41107 and No. DE-SC0011689; National
Science Foundation, Grant No. 0450696; The Grainger Foundation; Marie Curie-
IRSES/EPLANET; European Particle Physics Latin American Network; and UNESCO.

The authors gratefully acknowledge the support of NVIDIA Corporation with the
donation of GPU hardware used for this research.

## References

  * [1] Pierre Auger collaboration, _Features of the energy spectrum of cosmic rays above $2.5{\times}10^{18}$ eV using the Pierre Auger Observatory_, _Phys. Rev. Lett._ 125 (2020) 121106 [2008.06488].
  * [2] Pierre Auger collaboration, _A measurement of the cosmic-ray energy spectrum above $2.5{\times}10^{18}$ eV using the Pierre Auger Observatory_, 2008.06486.
  * [3] Pierre Auger collaboration, _The Pierre Auger Cosmic Ray Observatory_ , _Nucl. Instrum. Meth. A_ 798 (2015) 172 [1502.01323].
  * [4] J. Linsley and L. Scarsi, _Arrival Times of Air Shower Particles at Large Distances from the Axis_ , _Physical Review_ 128 (1962) 2384.
  * [5] T. K. Gaisser, R. Engel and E. Resconi, _Cosmic Rays and Particle Physics: 2nd Edition_. Cambridge University Press, 6, 2016.
  * [6] J. Matthews, _A Heitler model of extensive air showers_ , _Astropart. Phys._ 22 (2005) 387.
  * [7] Pierre Auger collaboration, _Muons in Air Showers at the Pierre Auger Observatory: Mean Number in Highly Inclined Events_ , _Phys. Rev. D_ 91 (2015) 032003 [1408.1421].
  * [8] Pierre Auger collaboration, _Muons in Air Showers at the Pierre Auger Observatory: Measurement of Atmospheric Production Depth_ , _Phys. Rev. D_ 90 (2014) 012012 [1407.5919].
  * [9] Pierre Auger collaboration, _Inferences on mass composition and tests of hadronic interactions from 0.3 to 100 EeV using the water-Cherenkov detectors of the Pierre Auger Observatory_ , _Phys. Rev. D_ 96 (2017) 122003 [1710.07249].
  * [10] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby et al., _Machine learning and the physical sciences_ , _Rev. Mod. Phys._ 91 (2019) 045002 [1903.10563].
  * [11] _A Living Review of Machine Learning for Particle Physics_ , https://iml-wg.github.io/HEPML-LivingReview/.
  * [12] A. Guillen, A. Bueno, J. Carceller, J. Martinez-Velazquez, G. Rubio, C. Todero Peixoto et al., _Deep learning techniques applied to the physics of extensive air showers_ , _Astropart. Phys._ 111 (2019) 12 [1807.09024].
  * [13] A. Karpathy, _The Unreasonable Effectiveness of Recurrent Neural Networks_ , (2015), https://karpathy.github.io/2015/05/21/rnn-effectiveness/.
  * [14] S. Hochreiter and J. Schmidhuber, _Long Short-Term Memory_ , _Neural Comput._ 9 (1997) 1735.
  * [15] C. Olah, _Understanding LSTM Networks_ , (2015), https://colah.github.io/posts/2015-08-Understanding-LSTMs/.
  * [16] Pierre Auger collaboration, _Azimuthal Asymmetry in the Risetime of the Surface Detector Signals of the Pierre Auger Observatory_ , _Phys. Rev. D_ 93 (2016) 072006 [1604.00978].
  * [17] Pierre Auger collaboration, _Reconstruction of Events Recorded with the Surface Detector of the Pierre Auger Observatory_ , _Journal of Instrumentation_ 15 (2020) P10021 [2007.09035].
  * [18] D. Heck, J. Knapp, J. Capdevielle, G. Schatz and T. Thouw, _CORSIKA: A Monte Carlo code to simulate extensive air showers_ , .
  * [19] T. Pierog, I. Karpenko, J. Katzy, E. Yatsenko and K. Werner, _EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider_ , _Phys. Rev. C_ 92 (2015) 034906 [1306.0121].
  * [20] Pierre Auger collaboration, _The offline software framework of the Pierre Auger Observatory_ , in _2005 IEEE Nuclear Science Symposium and Medical Imaging Conference_ , no. 2, pp. 1000–1003, 2005, DOI [astro-ph/0601016].
  * [21] D. P. Kingma and J. Ba, _Adam: A Method for Stochastic Optimization_ , 1412.6980.
  * [22] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau et al., _Array programming with NumPy_ , _Nature_ 585 (2020) 357–362.
  * [23] Wes McKinney, _Data Structures for Statistical Computing in Python_ , _Proceedings of the 9th Python in Science Conference_ (2010) 56.
  * [24] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan et al., _Pytorch: An imperative style, high-performance deep learning library_ , in _Advances in Neural Information Processing Systems 32_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox and R. Garnett, eds., pp. 8024–8035, Curran Associates, Inc., (2019), http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * [25] S. Ostapchenko, _Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model_ , _Phys. Rev. D_ 83 (2011) 014018 [1010.1869].
  * [26] E.-J. Ahn, R. Engel, T. K. Gaisser, P. Lipari and T. Stanev, _Cosmic ray interaction event generator SIBYLL 2.1_ , _Phys. Rev. D_ 80 (2009) 094003 [0906.4113].
  * [27] AGASA collaboration, _Muons ( >= 1-GeV) in large extensive air showers of energies between $10^{16.5}$ eV and $10^{19.5}$ eV observed at Akeno_, _J. Phys. G_ 21 (1995) 1101.
  * [28] K. Greisen, _Cosmic Ray showers_ , _Annual Review of Nuclear Science_ 10 (1960) 63.
  * [29] M. Nagano, M. Teshima, Y. Matsubara, H. Dai, T. Hara, N. Hayashida et al., _Energy spectrum of primary cosmic rays above $10^{17}$ eV determined from the extensive air shower experiment at Akeno_, _J. Phys. G_ 18 (1992) 423.
  * [30] Pierre Auger collaboration, _The Pierre Auger Observatory Upgrade - Preliminary Design Report_ , 1604.03637.

$\bullet$

## The Pierre Auger Collaboration

A. Aab80, P. Abreu72, M. Aglietta52,50, J.M. Albury12, I. Allekotte1, A.
Almela8,11, J. Alvarez-Muñiz79, R. Alves Batista80, G.A. Anastasi61,50, L.
Anchordoqui87, B. Andrada8, S. Andringa72, C. Aramo48, P.R. Araújo Ferreira40,
J. C. Arteaga Velázquez66, H. Asorey8, P. Assis72, G. Avila10, A.M. Badescu75,
A. Bakalova30, A. Balaceanu73, F. Barbato43,44, R.J. Barreira Luz72, K.H.
Becker36, J.A. Bellido12,68, C. Berat34, M.E. Bertaina61,50, X. Bertou1, P.L.
Biermannb, T. Bister40, J. Biteau35, J. Blazek30, C. Bleve34, M. Boháčová30,
D. Boncioli55,44, C. Bonifazi24, L. Bonneau Arbeletche19, N. Borodai69, A.M.
Botti8, J. Brackd, T. Bretz40, P.G. Brichetto Orchera8, F.L. Briechle40, P.
Buchholz42, A. Bueno78, S. Buitink14, M. Buscemi45, K.S. Caballero-Mora65, L.
Caccianiga57,47, F. Canfora80,82, I. Caracas36, J.M. Carceller78, R.
Caruso56,45, A. Castellina52,50, F. Catalani17, G. Cataldi46, L. Cazon72, M.
Cerda9, J.A. Chinellato20, K. Choi13, J. Chudoba30, L. Chytka31, R.W. Clay12,
A.C. Cobos Cerutti7, R. Colalillo58,48, A. Coleman93, M.R. Coluccia46, R.
Conceição72, A. Condorelli43,44, G. Consolati47,53, F. Contreras10, F.
Convenga54,46, D. Correia dos Santos26, C.E. Covault85, S. Dasso5,3, K.
Daumiller39, B.R. Dawson12, J.A. Day12, R.M. de Almeida26, J. de Jesús8,39,
S.J. de Jong80,82, G. De Mauro80,82, J.R.T. de Mello Neto24,25, I. De
Mitri43,44, J. de Oliveira26, D. de Oliveira Franco20, F. de Palma54,46, V. de
Souza18, E. De Vito54,46, M. del Río10, O. Deligny32, A. Di Matteo50, C.
Dobrigkeit20, J.C. D’Olivo67, R.C. dos Anjos23, M.T. Dova4, J. Ebr30, R.
Engel37,39, I. Epicoco54,46, M. Erdmann40, C.O. Escobara, A. Etchegoyen8,11,
H. Falcke80,83,82, J. Farmer92, G. Farrar90, A.C. Fauth20, N. Fazzinia, F.
Feldbusch38, F. Fenu52,50, B. Fick89, J.M. Figueira8, A. Filipčič77,76, T.
Fodran80, M.M. Freire6, T. Fujii92,e, A. Fuster8,11, C. Galea80, C.
Galelli57,47, B. García7, A.L. Garcia Vegas40, H. Gemmeke38, F. Gesualdi8,39,
A. Gherghel-Lascu73, P.L. Ghia32, U. Giaccari80, M. Giammarchi47, M. Giller70,
J. Glombitza40, F. Gobbi9, F. Gollan8, G. Golup1, M. Gómez Berisso1, P.F.
Gómez Vitale10, J.P. Gongora10, J.M. González1, N. González13, I. Goos1,39, D.
Góra69, A. Gorgi52,50, M. Gottowik36, T.D. Grubb12, F. Guarino58,48, G.P.
Guedes21, E. Guido50,61, S. Hahn39,8, P. Hamal30, M.R. Hampel8, P. Hansen4, D.
Harari1, V.M. Harvey12, A. Haungs39, T. Hebbeker40, D. Heck39, G.C. Hill12, C.
Hojvata, J.R. Hörandel80,82, P. Horvath31, M. Hrabovský31, T. Huege39,14, J.
Hulsman8,39, A. Insolia56,45, P.G. Isar74, P. Janecek30, J.A. Johnsen86, J.
Jurysek30, A. Kääpä36, K.H. Kampert36, B. Keilhauer39, J. Kemp40, H.O.
Klages39, M. Kleifges38, J. Kleinfeller9, M. Köpke37, N. Kunka38, B.L. Lago16,
R.G. Lang18, N. Langner40, M.A. Leigui de Oliveira22, V. Lenok39, A.
Letessier-Selvon33, I. Lhenry-Yvon32, D. Lo Presti56,45, L. Lopes72, R.
López62, L. Lu94, Q. Luce37, A. Lucero8, J.P. Lundquist76, A. Machado
Payeras20, G. Mancarella54,46, D. Mandat30, B.C. Manning12, J. Manshanden41,
P. Mantscha, S. Marafico32, A.G. Mariazzi4, I.C. Mariş13, G. Marsella59,45, D.
Martello54,46, H. Martinez18, O. Martínez Bravo62, M. Mastrodicasa55,44, H.J.
Mathes39, J. Matthews88, G. Matthiae60,49, E. Mayotte36, P.O. Mazura, G.
Medina-Tanco67, D. Melo8, A. Menshikov38, K.-D. Merenda86, S. Michal31, M.I.
Micheletti6, L. Miramonti57,47, S. Mollerach1, F. Montanet34, C. Morello52,50,
M. Mostafá91, A.L. Müller8, M.A. Muller20, K. Mulrey14, R. Mussa50, M.
Muzio90, W.M. Namasaka36, A. Nasr-Esfahani36, L. Nellen67, M. Niculescu-
Oglinzanu73, M. Niechciol42, D. Nitz89, D. Nosek29, V. Novotny29, L. Nožka31,
A Nucita54,46, L.A. Núñez28, M. Palatka30, J. Pallotta2, P. Papenbreer36, G.
Parente79, A. Parra62, M. Pech30, F. Pedreira79, J. Pȩkala69, R. Pelayo64, J.
Peña-Rodriguez28, E.E. Pereira Martins37,8, J. Perez Armand19, C. Pérez
Bertolli8,39, M. Perlin8,39, L. Perrone54,46, S. Petrera43,44, T. Pierog39, M.
Pimenta72, V. Pirronello56,45, M. Platino8, B. Pont80, M. Pothast82,80, P.
Privitera92, M. Prouza30, A. Puyleart89, S. Querchfeld36, J. Rautenberg36, D.
Ravignani8, M. Reininghaus39,8, J. Ridky30, F. Riehn72, M. Risse42, V.
Rizi55,44, W. Rodrigues de Carvalho19, J. Rodriguez Rojo10, M.J. Roncoroni8,
M. Roth39, E. Roulet1, A.C. Rovero5, P. Ruehl42, S.J. Saffi12, A. Saftoiu73,
F. Salamida55,44, H. Salazar62, G. Salina49, J.D. Sanabria Gomez28, F.
Sánchez8, E.M. Santos19, E. Santos30, F. Sarazin86, R. Sarmento72, C.
Sarmiento-Cano8, R. Sato10, P. Savina54,46,32, C.M. Schäfer39, V. Scherini46,
H. Schieler39, M. Schimassek37,8, M. Schimp36, F. Schlüter39,8, D. Schmidt37,
O. Scholten81,14, P. Schovánek30, F.G. Schröder93,39, S. Schröder36, J.
Schulte40, S.J. Sciutto4, M. Scornavacche8,39, A. Segreto51,45, S. Sehgal36,
R.C. Shellard15, G. Sigl41, G. Silli8,39, O. Sima73,f, R. Šmída92, P.
Sommers91, J.F. Soriano87, J. Souchard34, R. Squartini9, M. Stadelmaier39,8,
D. Stanca73, S. Stanič76, J. Stasielak69, P. Stassi34, A. Streich37,8, M.
Suárez-Durán28, T. Sudholz12, T. Suomijärvi35, A.D. Supanitsky8, J. Šupík31,
Z. Szadkowski71, A. Taboada37, A. Tapia27, C. Taricco61,50, C.
Timmermans82,80, O. Tkachenko39, P. Tobiska30, C.J. Todero Peixoto17, B.
Tomé72, A. Travaini9, P. Travnicek30, C. Trimarelli55,44, M. Trini76, M.
Tueros4, R. Ulrich39, M. Unger39, L. Vaclavek31, M. Vacula31, J.F. Valdés
Galicia67, L. Valore58,48, E. Varela62, V. Varma K.C.8,39, A. Vásquez-
Ramírez28, D. Veberič39, C. Ventura25, I.D. Vergara Quispe4, V. Verzi49, J.
Vicha30, J. Vink84, S. Vorobiov76, H. Wahlberg4, C. Watanabe24, A.A. Watsonc,
M. Weber38, A. Weindl39, L. Wiencke86, H. Wilczyński69, T. Winchen14, M.
Wirtz40, D. Wittkowski36, B. Wundheiler8, A. Yushkov30, O. Zapparrata13, E.
Zas79, D. Zavrtanik76,77, M. Zavrtanik77,76, L. Zehrer76, A. Zepeda63

1

Centro Atómico Bariloche and Instituto Balseiro (CNEA-UNCuyo-CONICET), San
Carlos de Bariloche, Argentina

2

Centro de Investigaciones en Láseres y Aplicaciones, CITEDEF and CONICET,
Villa Martelli, Argentina

3

Departamento de Física and Departamento de Ciencias de la Atmósfera y los
Océanos, FCEyN, Universidad de Buenos Aires and CONICET, Buenos Aires,
Argentina

4

IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina

5

Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA), Buenos
Aires, Argentina

6

Instituto de Física de Rosario (IFIR) – CONICET/U.N.R. and Facultad de
Ciencias Bioquímicas y Farmacéuticas U.N.R., Rosario, Argentina

7

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET,
UNSAM), and Universidad Tecnológica Nacional – Facultad Regional Mendoza
(CONICET/CNEA), Mendoza, Argentina

8

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET,
UNSAM), Buenos Aires, Argentina

9

Observatorio Pierre Auger, Malargüe, Argentina

10

Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe,
Argentina

11

Universidad Tecnológica Nacional – Facultad Regional Buenos Aires, Buenos
Aires, Argentina

12

University of Adelaide, Adelaide, S.A., Australia

13

Université Libre de Bruxelles (ULB), Brussels, Belgium

14

Vrije Universiteit Brussels, Brussels, Belgium

15

Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil

16

Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Nova Friburgo,
Brazil

17

Universidade de São Paulo, Escola de Engenharia de Lorena, Lorena, SP, Brazil

18

Universidade de São Paulo, Instituto de Física de São Carlos, São Carlos, SP,
Brazil

19

Universidade de São Paulo, Instituto de Física, São Paulo, SP, Brazil

20

Universidade Estadual de Campinas, IFGW, Campinas, SP, Brazil

21

Universidade Estadual de Feira de Santana, Feira de Santana, Brazil

22

Universidade Federal do ABC, Santo André, SP, Brazil

23

Universidade Federal do Paraná, Setor Palotina, Palotina, Brazil

24

Universidade Federal do Rio de Janeiro, Instituto de Física, Rio de Janeiro,
RJ, Brazil

25

Universidade Federal do Rio de Janeiro (UFRJ), Observatório do Valongo, Rio de
Janeiro, RJ, Brazil

26

Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil

27

Universidad de Medellín, Medellín, Colombia

28

Universidad Industrial de Santander, Bucaramanga, Colombia

29

Charles University, Faculty of Mathematics and Physics, Institute of Particle
and Nuclear Physics, Prague, Czech Republic

30

Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic

31

Palacky University, RCPTM, Olomouc, Czech Republic

32

CNRS/IN2P3, IJCLab, Université Paris-Saclay, Orsay, France

33

Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Sorbonne
Université, Université de Paris, CNRS-IN2P3, Paris, France

34

Univ. Grenoble Alpes, CNRS, Grenoble Institute of Engineering Univ. Grenoble
Alpes, LPSC-IN2P3, 38000 Grenoble, France

35

Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France

36

Bergische Universität Wuppertal, Department of Physics, Wuppertal, Germany

37

Karlsruhe Institute of Technology (KIT), Institute for Experimental Particle
Physics, Karlsruhe, Germany

38

Karlsruhe Institute of Technology (KIT), Institut für Prozessdatenverarbeitung
und Elektronik, Karlsruhe, Germany

39

Karlsruhe Institute of Technology (KIT), Institute for Astroparticle Physics,
Karlsruhe, Germany

40

RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany

41

Universität Hamburg, II. Institut für Theoretische Physik, Hamburg, Germany

42

Universität Siegen, Department Physik – Experimentelle Teilchenphysik, Siegen,
Germany

43

Gran Sasso Science Institute, L’Aquila, Italy

44

INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy

45

INFN, Sezione di Catania, Catania, Italy

46

INFN, Sezione di Lecce, Lecce, Italy

47

INFN, Sezione di Milano, Milano, Italy

48

INFN, Sezione di Napoli, Napoli, Italy

49

INFN, Sezione di Roma “Tor Vergata”, Roma, Italy

50

INFN, Sezione di Torino, Torino, Italy

51

Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo (INAF), Palermo,
Italy

52

Osservatorio Astrofisico di Torino (INAF), Torino, Italy

53

Politecnico di Milano, Dipartimento di Scienze e Tecnologie Aerospaziali ,
Milano, Italy

54

Università del Salento, Dipartimento di Matematica e Fisica “E. De Giorgi”,
Lecce, Italy

55

Università dell’Aquila, Dipartimento di Scienze Fisiche e Chimiche, L’Aquila,
Italy

56

Università di Catania, Dipartimento di Fisica e Astronomia, Catania, Italy

57

Università di Milano, Dipartimento di Fisica, Milano, Italy

58

Università di Napoli “Federico II”, Dipartimento di Fisica “Ettore Pancini”,
Napoli, Italy

59

Università di Palermo, Dipartimento di Fisica e Chimica ”E. Segrè”, Palermo,
Italy

60

Università di Roma “Tor Vergata”, Dipartimento di Fisica, Roma, Italy

61

Università Torino, Dipartimento di Fisica, Torino, Italy

62

Benemérita Universidad Autónoma de Puebla, Puebla, México

63

Centro de Investigación y de Estudios Avanzados del IPN (CINVESTAV), México,
D.F., México

64

Unidad Profesional Interdisciplinaria en Ingeniería y Tecnologías Avanzadas
del Instituto Politécnico Nacional (UPIITA-IPN), México, D.F., México

65

Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, México

66

Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Michoacán, México

67

Universidad Nacional Autónoma de México, México, D.F., México

68

Universidad Nacional de San Agustin de Arequipa, Facultad de Ciencias
Naturales y Formales, Arequipa, Peru

69

Institute of Nuclear Physics PAN, Krakow, Poland

70

University of Łódź, Faculty of Astrophysics, Łódź, Poland

71

University of Łódź, Faculty of High-Energy Astrophysics,Łódź, Poland

72

Laboratório de Instrumentação e Física Experimental de Partículas – LIP and
Instituto Superior Técnico – IST, Universidade de Lisboa – UL, Lisboa,
Portugal

73

“Horia Hulubei” National Institute for Physics and Nuclear Engineering,
Bucharest-Magurele, Romania

74

Institute of Space Science, Bucharest-Magurele, Romania

75

University Politehnica of Bucharest, Bucharest, Romania

76

Center for Astrophysics and Cosmology (CAC), University of Nova Gorica, Nova
Gorica, Slovenia

77

Experimental Particle Physics Department, J. Stefan Institute, Ljubljana,
Slovenia

78

Universidad de Granada and C.A.F.P.E., Granada, Spain

79

Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago
de Compostela, Santiago de Compostela, Spain

80

IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands

81

KVI – Center for Advanced Radiation Technology, University of Groningen,
Groningen, The Netherlands

82

Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science
Park, Amsterdam, The Netherlands

83

Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The
Netherlands

84

Universiteit van Amsterdam, Faculty of Science, Amsterdam, The Netherlands

85

Case Western Reserve University, Cleveland, OH, USA

86

Colorado School of Mines, Golden, CO, USA

87

Department of Physics and Astronomy, Lehman College, City University of New
York, Bronx, NY, USA

88

Louisiana State University, Baton Rouge, LA, USA

89

Michigan Technological University, Houghton, MI, USA

90

New York University, New York, NY, USA

91

Pennsylvania State University, University Park, PA, USA

92

University of Chicago, Enrico Fermi Institute, Chicago, IL, USA

93

University of Delaware, Department of Physics and Astronomy, Bartol Research
Institute, Newark, DE, USA

94

University of Wisconsin-Madison, Department of Physics and WIPAC, Madison, WI,
USA

—–

a

Fermi National Accelerator Laboratory, Fermilab, Batavia, IL, USA

b

Max-Planck-Institut für Radioastronomie, Bonn, Germany

c

School of Physics and Astronomy, University of Leeds, Leeds, United Kingdom

d

Colorado State University, Fort Collins, CO, USA

e

now at Hakubi Center for Advanced Research and Graduate School of Science,
Kyoto University, Kyoto, Japan

f

also at University of Bucharest, Physics Department, Bucharest, Romania


# Sector decomposition via computational geometry

Toshiaki Kaneko
High Energy Accelerator Research Organization (KEK)
1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan
Takahiro Ueda
Graduate School of Pure and Applied Sciences, University of Tsukuba,
1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan

###### Abstract

A non-iterative method is presented for the factorization step of sector
decomposition method, which separates infrared divergent part from loop
integration. This method is based on a classification of asymptotic behavior
of polynomials. The problem is converted to ones for convex body in Euclidean
space. They are solved with algorithms developed in computational geometry. A
test implementation shows that this method produces less number of decomposed
sectors than usual iterative sector decompositions.

Talk given at 13th International Workshop on Advanced Computing and Analysis
Techniques in Physics Research, February 26, 2010, Jaipur, India

##  1 Introduction

In the calculation of Feynman amplitude with massless particles, one has to
regulate infrared divergences (IR), which cancel out among loop corrections
and real emission processes. In perturbative QCD, this divergences are
regulated by $D$-dimensional method. Divergent part is expressed as poles in
terms of ${\epsilon}=(4-D)/2$. Sector decomposition method developed in Refs.
[1, 2, 3, 4] are widely used for separating IR divergences.

As a simple example of separation of infrared divergences, let us consider the
following one-dimensional integral for ${\epsilon}<0$ with a regular function
$f$ such that $f(0)$ is a non-zero finite value:

$\begin{split}I&=\int_{0}^{1}dx\;x^{-1-{\epsilon}}f(x)\\\
&=\int_{0}^{1}dx\;x^{-1-{\epsilon}}f(0)+\int_{0}^{1}dx\;x^{-1-{\epsilon}}\left(f(x)-f(0)\right)\\\
&=-\frac{f(0)}{{\epsilon}}+\int_{0}^{1}dx\;\left(x^{-{\epsilon}}f^{\prime}(0)+\frac{1}{2}x^{1-{\epsilon}}f^{\prime\prime}(0)+...\right).\end{split}$
(1)

When ${\epsilon}=0$, the factor $x^{-1-{\epsilon}}$ produces a logarithmic
divergence by the integration around $x\sim 0$. While for ${\epsilon}<0$, the
divergence is regularized and converted to a pole of ${\epsilon}$. The
coefficient of the pole is $f(0)$, which is the first term of the Taylor
expansion of $f(x)$ around $x=0$. The rest of the integration becomes finite
for ${\epsilon}\rightarrow 0$. This example shows that, when the singular part
of the integrand is factored out, one can separate divergent part as poles in
terms of ${\epsilon}$.

For the general loop integrations, it will also be possible to separate IR
divergence, if the singular part is factored out. Although the factorization
of a multi-variate polynomial is not a trivial problem, it is solved by the
method of sector decomposition. This method shows that when the integration
domain is properly decomposed into sectors and selecting good variables in
each sector, the divergent part is factored out. This method provides at the
same time a practical procedure to find such sectors and variables. Then we
can obtain the coefficients of the poles of ${\epsilon}$ and the finite part
as a Laurent series in terms of ${\epsilon}$.

Usual sector decomposition method divides the integration domain and find
appropriate new variables by an iterative way. It was found that a simple
iterative method may fall into an infinite loop and strategies are proposed to
avoid this problem by Refs. [5, 6]. Other strategies have been proposed and a
practical system has been constructed by Refs. [7, 8, 9]. Since one has to
calculate integration for each sector, a method is preferable when it produces
less number of sectors.

We propose another method based on the classification of asymptotic behavior
of polynomials around the origin. The problem is converted to one for convex
bodies in Euclidean space and solved with algorithms developed in
computational geometry. These are deterministic algorithms without iterations.
This talk presents, based on Ref. [10], our basic ideas and show the results
of this method.

##  2 Sector decomposition

The procedure of sector decomposition consists of the following steps:

  1. 1.

Start from Feynman parameter representation of a loop integration.

  2. 2.

Primary sector decomposition:

The $\delta$-function is integrated and the integration domain is decomposed
such that the singularities appear in the specific position in each sector.

  3. 3.

Factorization of the integrand:

The sectors obtained in the previous step are decomposed into finer ones. New
variables are found such that singular parts are factored out.

  4. 4.

Separation of poles in terms of ${\epsilon}$:

The regular part of the integrand is expanded in terms of ${\epsilon}$. The
coefficients of Laurent series in terms of ${\epsilon}$ are obtained.

  5. 5.

Integration of coefficients:

The coefficients are now free from IR divergences. However, they will still
include multi-dimensional integration. They will be integrated out by analytic
or numerical methods.

Figure 1: An example of primary sector decomposition in 3-dimensional case.
Integration plane is divided into 3 sectors and sub-domain $abcd$ in (a) is
mapped to a square shown by (b).

Feynman parameter representation of a loop integration includes
$\delta$-function which defines the $(N-1)$-dimensional hyperplane. The
primary sector decomposition integrates over one of integration variables and
decompose the integration domain into $N$ sectors. An example is shown by Fig.
1 for the 3-dimensional case. By the decomposition of the integration domain
and selecting new variables, sectors become $(N-1)$-dimensional cubes. The
original boundary of the integration domain is mapped to a part of the new
boundary specified by conditions that some of new variables are 0. The
resulting expression of the integration for $l$-th sector becomes:

$G_{l}=\int_{0}^{1}d^{N-1}t\;t^{\nu-\mathbb{I}}\;{\cal
U}_{l}^{\gamma}(t)\,{\cal F}_{l}^{\beta}(t),\quad
t^{\nu-\mathbb{I}}=\prod_{j=0}^{N-1}t_{j}^{\nu_{j}-1},$ (2)

where ${\cal F}_{l}(t)$ and ${\cal U}_{l}(t)$ are polynomials of variables
$t$. Parameter $\nu_{j}$ is integer determined by the power of propagators,
and $\mathbb{I}$ is an $(N-1)$-dimensional vector whose elements are all 1.
Powers $\beta$ and $\gamma$ are functions of ${\epsilon}$. If ${\cal
F}_{l}(t)$ or ${\cal U}_{l}(t)$ become 0 at the boundary, and corresponding
power $\beta$ or $\gamma$ becomes negative integer for ${\epsilon}=0$, the
integration diverges.

##  3 Geometric method

We want to express polynomial ${\cal F}_{l}(t)$ in the form:

${\cal F}_{l}(t(z))=C_{a}z^{b_{a}}\left(1+H_{a}(z)\right),\quad
H_{a}(0)=0,\quad z^{b_{a}}=\prod_{j}z_{j}^{(b_{a})_{j}},$ (3)

where $z$ is a set of new variable, $H_{a}(z)$ is a polynomial of $z$, $b_{a}$
is an integer vector and $C_{a}$ is a constant. In order to seek such an
expression, we consider the following example taken from one-loop box
integration:

$\begin{split}{\cal
F}_{l}(t)=-s_{23}t_{2}t_{3}-s_{12}t_{1}-s_{4}t_{1}t_{3}=t_{2}t_{3}\left(-s_{23}-s_{12}\frac{t_{1}}{t_{2}t_{3}}-s_{4}\frac{t_{1}}{t_{2}}\right).\end{split}$
(4)

If terms $t_{1}/t_{2}t_{3}$ and $t_{1}/t_{2}\rightarrow 0$ when
$t_{2}t_{3}\rightarrow 0$, the asymptotic behavior of ${\cal F}_{l}$ is
determined by the term $-s_{23}t_{2}t_{3}$. The singular behavior of the
integrand around the origin is determined by this term. This condition is
satisfied by taking new variables $z$ defined by
$t_{2}=z_{2},\;t_{3}=z_{3},\;t_{1}=z_{1}z_{2}z_{3}$. With the range of $z_{j}$
being limited to $(0,1)$, a sub-domain of the integration domain in $t$-space
is obtained. This example shows that a term of a polynomial in some sub-domain
will determine asymptotic behavior of the polynomial around the origin. We
call this term dominant in this sub-domain.

In order to see dominant terms more closely around the origin, let us change
variable to $y_{j}=-\log t_{j}$ or $t_{j}=e^{-y_{j}}$. Monomial
$t^{c}=t_{1}^{c_{1}}t_{2}^{c_{2}}\cdots t_{N-1}^{c_{N-1}}$ becomes
$e^{-(c,y)}$, where $(c,y)=\sum{c_{j}y_{j}}$ is the inner product defined in
$(N-1)$-dimensional Euclidean space. This shows that a monomial corresponds
one-to-one to a integer vector $c$ in this Euclidean space (let us call this
power vector). Let $Z^{{\cal F}_{l}}$ is the set of power vectors
corresponding to the terms included in a polynomial ${\cal F}_{l}$. Polynomial
${\cal F}_{l}$ is expressed by:

${\cal F}_{l}(t)=\sum_{c\in Z^{{\cal F}_{l}}}a_{c}t^{c}=\sum_{c\in Z^{{\cal
F}_{l}}}a_{c}e^{-(c,y)}.$ (5)

In order to find the dominant term, let us consider a limit
$\lambda\rightarrow+\infty$ for $y_{j}=\lambda u_{j}\rightarrow+\infty$ with a
fixed non-negative real vector $u$. A term $t^{b}=e^{-(b,y)}$ is dominant in
this limit when $(b,y)\leq(c,y)$ for all $c\in Z^{{\cal F}_{l}}$. Conversely,
let us fix vector $b$ and vary $u$ or $y$. The term with this power vector $b$
is dominant in the sub-domain defined by:

$\Delta_{b}^{{\cal F}_{l}}:=\\{y\in\mathbb{R}_{\geq 0}^{N-1}|(c-b,y)>0,\forall
c\in Z^{{\cal F}_{l}}\\}.$ (6)

This sub-domain forms a convex polyhedral cone. The polynomial ${\cal F}_{l}$
is expressed as the following:

${\cal F}_{l}(t(y))=\sum_{b\in Z^{{\cal F}_{l}}}\theta(y\in\Delta_{b}^{{\cal
F}_{l}})e^{-(b,y)}\Bigl{[}a_{b}+\sum_{c\in Z^{{\cal
F}_{l}}-\\{b\\}}a_{c}e^{-(c-b,y)}\Bigr{]}.$ (7)

For a vector $y\in\Delta_{b}^{{\cal F}_{l}}$, factor $e^{-(c-b,y)}\rightarrow
0$ for $|y|\rightarrow\infty$, since $(c-b,y)>0$ holds for all $c$. This
implies that the term $t^{b}$ is dominant in the sub-domain $\Delta_{b}^{{\cal
F}_{l}}$. Thus, asymptotic behavior of ${\cal F}_{l}$ is classified in terms
of these sub-domains $\\{\Delta_{b}^{{\cal F}_{l}}\\}$.

Figure 2: An example of geometric sector decomposition

When this method is applied also to ${\cal U}_{l}$, the singular parts are
factored out from the integration. Let us return to our example Eq. (4) of
one-loop box. In this case ${\cal U}_{l}$ is $1+t_{1}+t_{2}+t_{3}$ and this
expression shows that ${\cal U}_{l}$ does not produce IR singularity (${\cal
U}_{l}\rightarrow 1$ for $t_{i}\rightarrow 0$). With a simple calculation, we
obtain from Eq. (4):

$\begin{split}t^{b_{0}}&=t_{1},\quad t^{b_{1}}=t_{2}t_{3},\quad
t^{b_{2}}=t_{1}t_{3},\\\
Z^{\mathcal{F}_{l}}&=\\{b_{0}=(1,0,0),\;b_{1}=(0,1,1),\;b_{2}=(1,0,1)\\},\\\
\Delta_{b_{0}}&=\\{x_{1}(1,1,0)+x_{2}(1,0,1)+x_{3}(0,1,0)+x_{4}(0,0,1)\;|\;x_{1},x_{2},x_{3},x_{4}\geq
0\\},\\\
\Delta_{b_{1}}&=\\{x_{1}(1,0,0)+x_{2}(1,1,0)+x_{3}(1,0,1)\;|\;x_{1},x_{2},x_{3}\geq
0\\},\\\ \Delta_{b_{2}}&=\emptyset.\end{split}$ (8)

The sub-domains are shown in Fig. 2. The fact that the last sub-domain is
empty means that term $t_{1}t_{3}$ never becomes dominant. The second domain
is a triangular cone for which variables $x$ take simple range of values.
However, the first sub-domain includes four parameters in 3-dimensional space.
When this sub-domain is divided further into two triangular cones, this
redundancy disappears:

$\displaystyle\Delta_{b_{0}}$
$\displaystyle=\Delta_{b_{0}}^{(1)}\cup\Delta_{b_{0}}^{(2)},$
$\displaystyle\Delta_{b_{0}}^{(1)}$
$\displaystyle=\\{x_{1}(1,1,0)+x_{2}(1,0,1)+x_{4}(0,0,1)\;|\;x_{1},x_{2},x_{4}\geq
0\\},$ $\displaystyle\Delta_{b_{0}}^{(2)}$
$\displaystyle=\\{x_{1}(1,1,0)+x_{3}(0,1,0)+x_{4}(0,0,1)\;|\;x_{1},x_{3},x_{4}\geq
0\\}.$

Now we subsequently change variable from $y_{j}$ to $x_{i}$ and then to
$z_{i}=e^{-x_{i}}$. The integration domain for $z$ becomes 3-dimensional unit
cube. Finally we obtain a sector decomposition as the following:

$\displaystyle\begin{array}[]{lll}\hline\cr\text{sector}&~{}\quad\text{variables}&~{}\quad
Z^{\mathcal{F}_{l}}\\\
\hline\cr\Delta_{b_{0}}^{(1)}&t_{1}=z_{1}z_{2},\;t_{2}=z_{1},\;t_{3}=z_{2}z_{4}&z_{1}z_{2}(-s_{12}-s_{23}z_{4}-s_{4}z_{2}z_{4})\\\
\Delta_{b_{0}}^{(2)}&t_{1}=z_{1},\;t_{2}=z_{1}z_{3},\;t_{3}=z_{4}&z_{1}(-s_{12}-s_{23}z_{3}z_{4}-s_{4}z_{4})\\\
\Delta_{b_{1}}&t_{1}=z_{1}z_{2}z_{3},\;t_{2}=z_{2},\;t_{3}=z_{3}&z_{2}z_{3}(-s_{12}z_{1}-s_{23}-s_{4}z_{1}z_{3})\\\
\hline\cr\end{array}$

It is easy to calculate Jacobian, which is found to be a monomial of $z$.

For the general case, the decomposition of the integration domain in
accordance with the classification of the asymptotic behavior, corresponding
to Eqs. (6) and (8), becomes a problem in Euclidean geometry. It is not so
hard to solve this when one uses convex hull algorithms developed in
computational geometry. This decomposition is uniquely determined once a
polynomial is given. It is independent of the choice of an algorithm. In order
to obtain final representation of sector decomposition, corresponding to the
previous table, it is necessary to triangulate convex polyhedral cones. This
is performed with triangulation algorithms. The triangulation of convex
polyhedral cones is not determined uniquely and the number of sectors will
depend on the algorithm.

##  4 Test implementation and conclusion

We have made a test implementation of this method. Our convex hull algorithm
is a modified one for convex polyhedral cones based on the algorithm described
in Ref. [11] for polytopes. For triangulation of convex polyhedral cones, we
have developed our own algorithm. The input to our program is given by the
package described in Ref. [12] and output is passed to the same package for
the separation of poles and integrations of coefficients. We have checked our
program by comparison with another package qhull (Ref. [13]) for the convex
hull algorithm, calculating integration volume for the triangulation, and
comparison of integrated values of several diagrams with references. The
number of decomposed sectors are shown by Table 1.

Table 1:  Comparison of number of sectors among different methods. Numbers in column “H” are cited from Ref. [4]. Columns “A”, “B”, “C”, “S” and “X” indicate corresponding strategy described in [5] and [7]. As shown in [7], “F” means that the sector decomposition fails and “M” means that the memory overflow happened during the sector decomposition on a 8Gb machine. The numbers with “*” are given by [14]. “This method” indicates the number of sectors obtained by our method. “Exponential S.D.” indicates the number of sectors before the triangulation.  Diagram | A | B | C | S | X | H | This | Exponential
---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  | method | S.D.
Bubble | 2 | 2 | 2 | 2* | 2 |  | 2 | 2
Triangle | 3 | 3 | 3 | 3* | 3 |  | 3 | 3
Box | 12 | 12 | 12 | 12 | 12 |  | 12 | 8
Tbubble | 58 | 48 | 48 | 48* | 48 |  | 48 | 36
Double box, $p_{i}^{2}=0$ | 775 | 586 | 586 | 362 | 293 | 282 | 266 | 106
Double box, $p_{4}^{2}\neq 0$ | 543* | 245* | 245* | 230* | 192* | 197 | 186 | 100
Double box, $p_{i}^{2}=0$ | 1138 | 698 | 698 | 441* | 395 |  | 360 | 120
nonplanar |  |  |  |  |  |  |  |
D420 | 8898 | 564 | 564 | 180 | F |  | 168 | 100
3 loop vertex (A8) | 4617* | 1196* | 1196* | 871* | 750* | 684 | 684 | 240
Triple box | M | 114256 | 114256 | 22657 | 10155 |  | 6568 | 856

We have proposed a factorization algorithm in sector decomposition. Our method
is based on a classification of asymptotic behavior of polynomials. The
problem is converted to one for convex bodies in Euclidean space. In order to
find sector decomposition for a given integrand, we employed algorithms
developed in computational geometry. This method is deterministic and never
falls into an infinite loop. A test implementation shows that the number of
decomposed sectors is less than iterated sector decomposition combined with
several strategies.

This work is supported in part by Ministry of Education, Science, and Culture,
Japan under Grant-in-Aid Nos. 20340063 and 21540286.

## References

  * [1] T. Binoth and G. Heinrich, Nucl. Phys. B 585 (2000) 741 [arXiv:hep-ph/0004013].
  * [2] T. Binoth and G. Heinrich, Nucl. Phys. B 680 (2004) 375 [arXiv:hep-ph/0305234].
  * [3] T. Binoth and G. Heinrich, Nucl. Phys. B 693 (2004) 134 [arXiv:hep-ph/0402265].
  * [4] G. Heinrich, Int. J. Mod. Phys. A 23 (2008) 1457 [arXiv:0803.4177 [hep-ph]].
  * [5] C. Bogner and S. Weinzierl, Comput. Phys. Commun. 178 (2008) 596 [arXiv:0709.4092 [hep-ph]].
  * [6] C. Bogner and S. Weinzierl, Nucl. Phys. Proc. Suppl. 183 (2008) 256 [arXiv:0806.4307 [hep-ph]].
  * [7] A. V. Smirnov and M. N. Tentyukov, Comput. Phys. Commun. 180 (2009) 735 [arXiv:0807.4129 [hep-ph]].
  * [8] A. V. Smirnov and V. A. Smirnov, JHEP 0905 (2009) 004 [arXiv:0812.4700 [hep-ph]].
  * [9] A. V. Smirnov, V. A. Smirnov and M. Tentyukov, [arXiv:0912.0158 [hep-ph]].
  * [10] T. Kaneko and T. Ueda, [arXiv:0908.2897 [hep-ph]], to appear in Comput. Phys. Commun.
  * [11] H. Edelsbrunner. “Algorithms in Combinatorial Geometry,” Springer, 1987.
  * [12] T. Ueda and J. Fujimoto, PoS (ACAT08) (2008) 120 [arXiv:0902.2656 [hep-ph]].
  * [13] Program package qhull, `http://www.qhull.org/`.
  * [14] Private communication with A. V. Smirnov.