text:
\section{Introduction}
\label{se:Introduction}

\textsc{MaGe}\ (MAjorana-GErda) is a \textsc{Geant4}-based~\cite{Agostinelli:2002hh,Allison:2006} physics simulation
software
framework jointly developed by the \textsc{Majorana}\ and
\textsc{Gerda}\ collaborations~\cite{Aalseth:2004yt, Schoenert:2005}. Both experiments will
search for the neutrinoless double-beta decay ($0\nu\beta\beta$\ decay)
of the $^{76}\mathrm{Ge}$\ isotope using arrays of isotopically enriched High-Purity Germanium (HPGe) detectors.
The discovery of $0\nu\beta\beta$\ decay is the only
practical way to determine if the neutrino is a Majorana particle. For
further details on the physics motivation, see the review article
in~\cite{Avignone:2007fu}. The current lower limit on the $0\nu\beta\beta$\ decay
half-life of $^{76}\mathrm{Ge}$\ is
$1.9\times10^{25}\,\mathrm{years}$~\cite{Klap:2001}, making this decay extremely rare if it exists.
This requires great care to reduce experimental backgrounds from naturally occurring radioactivity
and cosmic rays. This is achieved via careful material selection and assay, a deep underground
location, passive and active shielding, and analysis cuts.
The purpose of \textsc{MaGe}\ is to simulate the response of the \textsc{Majorana}\ and \textsc{Gerda}\
experiments to ionizing radiation from backgrounds, calibration sources, and $0\nu\beta\beta$\ decays.
\textsc{MaGe}\ is also used to simulate the response of prototype detectors, test-stands, and low-background
assay systems.
In the prototyping phase, the simulation
is used as a virtual test stand to guide detector design, to estimate the
effectiveness of proposed background reduction techniques, and to estimate
the experimental sensitivity. During experimental operation, \textsc{MaGe}\ will be used to simulate and characterize unexpected
backgrounds and determine the ultimate sensitivity of the experiments. It will also provide
probability distributions for signal extraction analyses.
The combination of the two collaboration's simulation package into one framework
reduces duplication of efforts, eases comparison between simulated data and experiment, and simplifies the addition of new simulated detector geometries. \textsc{Majorana}\ and \textsc{Gerda}\ are currently constructing detectors with tens of kilograms of enriched isotopes, but it is the goal of parts of the two collaborations to merge and pursue a joint effort towards a tonne-scale germanium experiment. Having a joint simulation package during the very early phases will ease this future integration. This paper focusses on the software framework, custom event generators, and physics lists of \textsc{MaGe} .

The code and physics requirements are given in section~\ref{section:Requirements}. The code structure of \textsc{MaGe}\ is discussed in section~\ref{section:structure}.
Section~\ref{section:physics} describes the implemented \textsc{Geant4}\ physics lists that are optimized for low-energy (sub-keV to few MeV), low background applications. Validation of \textsc{MaGe}\
simulations against experimental data is discussed in Sect.~\ref{section:validation}.
Conclusions are provided in the last section.

\section{Requirements}
\label{section:Requirements}

The requirements for the \textsc{MaGe}\ framework can be subdivided into physics requirements and software requirements.
The physics requirements define the physical processes that have to be simulated to find the response of the detectors.
Software requirements are driven by the use of the \textsc{Geant4}\ toolkit as basis for \textsc{MaGe}\ and the anticipated end-users.

\subsection{Physics Requirements}
\label{se:physics_requirements}

\textsc{Geant4}\ is a simulation toolkit that uses Monte Carlo techniques to simulate the propagation of particles and nuclei through matter. It has extensive capabilities to simulate different experimental geometries, propagating particles, and particle interactions, and has the foundation of the physics requirements for \textsc{MaGe}.
The choice of \textsc{Geant4}\ over other packages as the basis for \textsc{MaGe}\ was motivated
by its flexibility and active development within the particle and
medical physics communities, as well as its C++ and object-oriented structure. \textsc{Geant4}\ is open-source and allows collaborations with members from multiple countries to use it. It is a standard simulation tool for LHC experiments and will likely be supported for at least another decade.

The most stringent
requirements for \textsc{MaGe}\ are the proper simulation of the relevant
background sources for $0\nu\beta\beta$\ decay
experiments. \textsc{Geant4}\ fulfills some of these requirements, specifically it simulates:
\begin{enumerate}
\item Electromagnetic interactions of electrons and $\gamma$-rays at MeV and keV energies.
\item Radioactive isotope decay chains and nuclear de-excitations.
\item Interactions of thermal and fast neutrons.
\item Electromagnetic and hadronic showers
initiated by cosmic-ray muons.
\item Penetration depths and ionization energy loss profiles of
$\alpha$-particles.
\end{enumerate}
The list of implemented physics models in \textsc{MaGe}\ was optimized for
low-background, underground physics
applications~\cite{Bauer:2004an}, with an emphasis on low-energy
interactions and hadronic interactions resulting from cosmic-ray
spallation.
\textsc{MaGe}\ has implemented different selections of models in its \textsc{Geant4}\ physics list. These are tailored to fit specific physics applications, and the required processes can be selected by the user to optimize computation time and data storage requirements. For example, the simulation of muon propagation through rocks requires different models than that of background from radioactive decays in detector components. Many
important tuning parameters, such as the	 production cuts for $\delta$-rays and soft bremsstrahlung photons, may also be set by the user.
Specific details of the physics lists implemented in \textsc{MaGe}\ are discussed in Section~\ref{section:physics}. The \textsc{MaGe}\ code is regularly updated and ported in order to make
it compatible with the most recent \textsc{Geant4}\ releases.\\

Other physics processes that are not part of \textsc{Geant4}\ had to be simulated. These include:
\begin{enumerate}
\item Electric-field solvers are required to simulate the trajectories of charge carriers inside the HPGe crystals under the influence of the biasing field.
\item Generators of electronic waveforms by charge carriers in the HPGe detectors as they drift inside the crystal towards the collection electrodes.
\item Electronic transfer functions of generated pulses into simulated detector pulses.
\item A generic surface sampler that uniformly and randomly samples points on arbitrary surfaces had to be implemented to simulate surface alpha contaminations. The algorithm of this sampler is described in Ref.~\cite{surfacesampler}.
\item Event generators for physics processes other than normal nuclear decay, such as two neutrino double-beta decay for different models, $0\nu\beta\beta$\ decay, and cosmic-ray muons at depth.
\end{enumerate}
Some of the existing aspects of \textsc{Geant4}\ had to be extended or improved to fulfill the requirements of \textsc{MaGe} . These are described later in this paper.

\subsection{Software Requirements}
\label{se:software_requirement}

The use of \textsc{Geant4}\ as the basis of \textsc{MaGe}\ made C++ the natural choice for \textsc{MaGe} , and \textsc{MaGe}\ makes full use of the object-oriented nature of C++ and \textsc{Geant4}{}. The software framework is required to:
\begin{enumerate}
\item Allow run-time selection of detector configuration, event generators and output format.
\item Allow parallel and independent
development of different branches of the code. For example, the development of the geometric descriptions of the \textsc{Gerda}\ and \textsc{Majorana}\ detector geometries should proceed independently and with minimal interference.
\item Ease-of-maintenance over the full life time of the
experiments.
\item Perform simulations with different
configurations by physicists that are inexperienced in programming.
\end{enumerate}
The collaboration does not maintain the \textsc{MaGe}\ source code for public
release, since most of components are experiment-specific and not useful beyond the two collaborations. However, some select components, such as physics lists and event generators, are provided to
interested parties on a case-by-case basis. The object-oriented nature of the \textsc{MaGe}\ framework makes the transfer of
such code to other users straightforward.

\section{Structure of \textsc{MaGe}\ }
\label{section:structure}

To realize the requirements outlined in the previous section we subdivided a simulation task into different  components. These components are geometries, generators, output formats, and physics lists. The user instantiates one instance of a class corresponding to each component at run-time, typically via a \textsc{Geant4}\ messenger in a macro file. A component may then also instantiate other components, and have its own macro commands that allow the user to further refine the simulation parameters. This design allows the simulation of many different detectors, prototypes and
validations to be performed within the same executable using
the same physics processes, geometries and other codes. This greatly eases cross-comparisons and reduces coding and debugging effort. The components are described in this section.

\subsection{Geometries}
\label{se:geometries}
These are the physical geometries of the detectors or experiments that are being simulated. \textsc{MaGe}\
supports geometric description via the \textsc{Geant4}\ geometry description classes or through an
interface with the Geometry Description Mark-up Language~\cite{gdml} (GDML).
\textsc{MaGe}\ currently has about 30 user-selectable geometries. The geometries
range from simple cylindrical crystals to a full detector array with 60 crystals, mounting
components, shield and surrounding room. Each geometry is encoded in a class that derives
from a geometry base class that contains the following basic components of a
geometry:
\begin{itemize}
\item A unique identifying serial number string.
\item A detector name string.
\item A \emph{ConstructDetector()} method that is invoked by \textsc{Geant4}\ to construct the detector geometry during run-time.
\item An associated \emph{G4LogicalVolume}.
\item A setting of the importance value of the region, used by \textsc{Geant4}\
when performing simulations that require importance sampling of geometries to optimize performance.
\end{itemize}
This design also allows the reuse of existing geometry classes, since
the classes describing a geometry can be instantiated within a class
that requires that component. For example, a detailed germanium
crystal has been coded that is used many times in other simulated
detector geometries. This crystal can be simulated on its own, or be
instantiated many times in a complex detector array. Shown in
Fig.~\ref{fig:MJDEMOGeometry} is an example of how multiple
stand-alone geometry classes are combined to create a complex detector. A rendering of the detector is shown in Fig.~\ref{fig:MJDEMOrendering}. The \textsc{Gerda}\ detector array rendering shown
in Fig.~\ref{fig:GerdaArray} is constructed in a similar way. A bonus
of this approach is that the \textsc{Gerda}\ and \textsc{Majorana}\ collaborations can share
the same basic geometries, such as crystals. This reduces redundant
code and increases code scrutiny.

\begin{figure}[tbh]
\flushleft
\includegraphics[width=0.4\textwidth]{plots/MJDEMOdiagram.pdf}
\caption{Classes used to create the geometry of the \textsc{Majorana}\ {\sc Demonstrator}. This detector consists of 35 crystals arranged in 7 columns of 5 crystals each. The crystals are placed inside a cryostat that in turn is placed inside a layered shield. This diagram illustrates how this design is realized in the simulation. Each of these classes can be instantiated on its own, or as part of a larger geometry, such as this. A rendering of this geometry is shown in Fig.~\ref{fig:MJDEMOrendering}.}\label{fig:MJDEMOGeometry}
\end{figure}

\begin{figure}[tbh]
\centering
\includegraphics[width=0.45\textwidth]{plots/MJDEMOrendering.pdf}
\caption{Rendering of \textsc{Majorana}\ {\sc Demonstrator} with cryostat and shields removed that was created using the class structure in Fig.~\ref{fig:MJDEMOGeometry}. Each of the 7 crystal columns with associated support structures is a stand-alone class, as is each of the 5 crystals in a column. }
\label{fig:MJDEMOrendering}
\end{figure}

\begin{figure}[tbh]
\centering
\includegraphics[width=0.45\textwidth]{plots/GerdaArray.pdf}
\caption{Rendering of \textsc{Gerda}\ detector array with cryostat and shields
removed. On the right hand side are 18-fold segmented HPGe detectors,
while left are unsegmented ones. Detectors of the same type are only
modeled once. The codes are reused to construct the whole array.}
\label{fig:GerdaArray}
\end{figure}

\subsection{Event Generators}
\label{se:event_generators}
These classes generate the initial conditions of each event to be simulated. \textsc{MaGe}\ has several event generators that all inherit from an \emph{MGVGenerator} base class. One such generator is instantiated by the user at run-time via a \textsc{Geant4}\ messenger. The \emph{MGVGenerator} base class contains the following relevant virtual methods that are defined in the daughter class:
\begin{itemize}
\item A \emph{BeginOfRunAction()} method that is executed at the beginning of a \textsc{Geant4}\ run. It is primarily used to compute lookup tables from analytic expressions for distributions to be sampled.
\item A \emph{EndOfRunAction()} method that is executed at the end of a run. It is primarily used to deallocate memory used by lookup tables.
\item A \emph{GeneratePrimaryVertex(G4Event *event)} method that passes the initial event particle type, position and
momentum to the \emph{G4Event} object.
\item A \emph{SetParticlePosition(G4ThreeVector vec)} method that can be used by another generator, such as the position sampler described below, to set the position for the initial vertex.
\end{itemize}
The generator may also have a \textsc{Geant4}\ messenger associated with it to allow the user to set parameters at run-time. Using this base class, we have created the following specific generators for \textsc{MaGe}:
\begin{itemize}
\item A wrapper for the Radioactive Decay Module~\cite{G4RDM}
(RDM) in \textsc{Geant4}\ that generates radioactive decays from unstable nuclei at a point.
\item ``Hybrid" generators that use the RDM generator to generate the initial decay of a nucleus, but then uses a customized surface or volume sampler to place it at a specific location in our simulated geometry. We developed a volume sampler that can generate points uniformly distributed in any {\sc
Geant4} solid. This is required to
simulate radioactive contamination embedded in detector components. A
surface sampler was also developed that creates points uniformly
distributed on the surface of any {\sc Geant4} solid.
This surface sampler is complex and is described further in~\cite{surfacesampler}.
A surface sampler is required to simulate surface contaminations, in particular $\alpha$-emitters and $\beta$-emitters on detector surfaces inside the cryostat that are within line-of-sight of the HPGe crystals.
\item Generators to simulate neutron and
muon backgrounds in underground laboratories, using either theoretical
models or data-driven approaches. Interfaces are available which read
initial conditions for an event from other codes, such as
\textsc{Sources4A}~\cite{sources4a} for neutron flux and
\textsc{Musun}~\cite{musun} for muon flux, as well as other selected analytical models, such as that by Wang et al.~\cite{wang}.
\item A wrapper for the FORTRAN-based \textsc{Decay0}~\cite{decay0} generator. It is primarily used to generate
different types of double-beta decays, but can also simulate normal nuclear decays. \textsc{Decay0} is able to
simulate properly the angular correlation between $\gamma$-rays in nuclear de-excitation cascades, that is not accounted for
in the \textsc{Geant4}\ Radioactive Decay Module.
\item A customized generator to simulate the neutron and gamma flux from an AmBe source used for studying neutron interactions.
\end{itemize}

\subsection{Physics Lists}
\label{se:physics_lists}
A collection of {\sc Geant4} physics
processes is called a physics list. They define the particles that are
included in the simulation and the decays and interactions they can
undergo.  There are several physics lists implemented in {\sc MaGe},
each optimized for the particular problem being simulated.  One list
is optimized at higher energies for simulating cosmic-ray muon
interactions, while others are optimized for standard electromagnetic
interactions at lower energies, i.e. below 10 MeV. These are used to
simulate the response of detectors to the decay of radioactive
isotopes. Each physics list is contained in its own class that is
instantiated at runtime. The physics lists and requirements are discussed in more detail in section~\ref{section:physics}.

\subsection{Output Format}
\label{se:output_format}
During the simulation \textsc{Geant4}\ generates complete information about the trajectory and interactions of particles as they propagate through the detector. Although all of this information is available to the user, it is typically processed, parsed and saved to an output file for further analysis after the simulation run is complete. Different detector or even different simulations of the same detectors have specific output requirements. \textsc{MaGe}\ defines a generic base output class and
does not have a built-in output format, but we have implemented interfaces to the AIDA-compliant~\cite{aida}
and ROOT~\cite{root} analysis tools, as well as simple text-based output.
At one extreme \textsc{MaGe}\ has implemented a class that saves all of the \textsc{Geant4}\ information for each step, at the other end is an output class that only generates a histogram of energy deposited in a single HPGe crystal. Each output format consists of a class that inherits from a virtual base class (\emph{MGVOutputManager}) that contains the following main methods:
\begin{itemize}
\item \emph{BeginOfRunAction()} and \emph{EndOfRunAction()} methods that are executed at the beginning and end of a simulation run respectively. They are used to open and close data files, create file formats,
i.e. ROOT trees, create histograms, and allocate and deallocate data structures.
\item \emph{BeginOfEventAction()} and \emph{EndOfEventAction()} methods that are executed at the beginning and end of the events. They are used to clear arrays that store stepping information and to perform processing on event data and fill event histograms.
\item A \emph{SteppingAction()} method that is executed at the end of each step. It gathers all the information concerning the particle in the given track and adds it to a histogram, or accumulates it in a variable.
\end{itemize}
These classes, in turn, can be
inherited by classes that simulate and store detector responses and
save any relevant information.

\subsection{Materials}
\label{se:materials}

\textsc{Geant4}\ defines materials internally in terms of how they interact
with ionizing radiation, but \textsc{MaGe}{} adds information about radioactive contaminants and other properties
of interest to \textsc{Majorana}\ and \textsc{Gerda}.
{\sc MaGe} has the ability to read in all
relevant information about materials from a PostgreSQL database. This
is currently limited to quantities required by \textsc{Geant4}\ such as density, isotopic
abundance, etc. Once the \textsc{Majorana} or \textsc{Gerda}
detectors are constructed,
the materials used will be carefully assayed and
characterized. All this information will be saved in a database as
well. {\sc MaGe} can then use this information to include the measured
activities in the simulation on a component-by-component and time-dependent level,
reducing systematic uncertainties in sensitivity calculations. \\

\subsection{Example}
\label{se:simple_example}

\textsc{MaGe}\ is typically compiled into a single executable.
The user selects the particular instance of the components they require during run-times via macros saved as text files.
\textsc{MaGe}\ macros are based on the \textsc{Geant4}\ \emph{messenger} classes. A simple, self-explanatory example of such a macro is given below. The simulation consists of a simple block that is bombarded with a beam of neutrons. The scattered neutrons are analyzed by the output class and elastic cross-sections are computed. This example is used to verify the neutron cross-sections implemented by \textsc{Geant4}.

{\small
\begin{verbatim}
/MG/geometry/detector solidBlock
/MG/geometry/solidblock/material Hydrogen
/MG/geometry/solidblock/edgeLength 1.0 cm

/MG/eventaction/rootschema HPNeutronTest

/MG/generator/select SPS

/run/initialize
/run/beamOn 500000
\end{verbatim}
}

\section{Physics} \label{section:physics}

The physics list in {\sc MaGe} has been optimized for the reliable
simulation of the most common background
sources in $0\nu\beta\beta$-decay experiments. It was designed
according to the suggestions of the \textsc{{\sc Geant4}}
team~\cite{physics_list} and optimized for
low-background physics
applications~\cite{Bauer:2004an,Pandola:2007}.
The default  \textsc{MaGe}\ physics list is mainly based on the
Underground Physics advanced example which is distributed with {\sc
Geant4}~\cite{physics_list}.

By default, \textsc{MaGe}\ uses the Low-Energy models based on the Livermore data
libraries~\cite{livermore1,livermore2} for the description of the electromagnetic
interactions of electrons, $\gamma$-rays and ions. These models include atomic
effects (e.g. fluorescence and Doppler broadening) and can handle interactions of electrons, positrons and photons
with energies down to 250~eV. \\
For specialized applications, electromagnetic interactions of $\gamma$-rays,
electrons and ions can be simulated in \textsc{MaGe}\ by the so-called ``standard models'' provided
by \textsc{Geant4}~\cite{physics-manual}. These models are tuned to high-energy physics applications;
they are less precise in the low energy region and do not include atomic effects.
However, they are faster in terms of computing time.

\textsc{MaGe}\ uses the standard \textsc{Geant4}\
models for the electromagnetic interactions of muons and of positrons.
Furthermore, synchrotron radiation is included in the physics list for electrons and positrons.
The electromagnetic physics processes provided by \textsc{Geant4}\ for $\gamma$-rays and e$^{\pm}$
(both "standard" and "low-energy") have been systematically validated by the
\textsc{Geant4}\ Collaboration~\cite{Amako:2005} and by other groups~\cite{Poon:2005} at the
few-percent level. The precision of the electromagnetic models for muons is discussed
in ~\cite{Bogdavov:2006}.

Interactions of e$^{\pm}$ and $\gamma$-rays with nuclei are simulated on the basis of the
equivalent photon approximation. The reaction cross section is determined according
to a parameterization from experimental data for all incident energies from the hadron
production thresholds upwards. Different energy regimes are considered separately
in the parameterization, as described in Ref.~\cite{physics-manual}. Hadronic final
states are generated using a Chiral Invariant Phase Space (CHIPS) decay
model~\cite{chips3}.

For energies above 3.5~GeV, the final states of photo-nuclear reactions are generated
according to a theory-based parton-string model, called the Quark-Gluon String Precompound
model (QGSP). The model is composed of several components: the quark-gluon string
(QGS) part handles the formation of the initial strings in the initial collision~\cite{qgsp1};
string fragmentation into hadrons is handled by the quark-gluon string fragmentation
model~\cite{qgsp2}, while the pre-compound model~\cite{qgsp3} takes care of
the de-excitation of the residual nucleus.

Hadronic interaction of muons with nuclei are managed by the \texttt{G4MuNuclearInteraction}
model. Muons produce virtual photons which are in turn converted to pions which interact with
the nucleus using a model derived from
the \textsc{Gheisha} code~\cite{gheisha}. The physics list also includes
the capture process of $\mu^{-}$ by nuclei. The hadronic interactions included in the \textsc{MaGe}\ physics list handle elastic scattering,
inelastic scattering, capture (for neutrons, $\pi^{-}$ and K$^{-}$), fission (neutrons only) and
decay.

The elastic scattering of all long-lived hadrons is described by the \texttt{G4LElastic} model,
which is based on the \textsc{Gheisha} code~\cite{gheisha} and includes a
parameterization for cross section and final state.
Elastic scattering of neutrons from thermal energies to 20~MeV is simulated
according to the data-driven \texttt{G4NeutronHPElastic} model, which is based on the tabulated
cross section and final state data from the ENDF/B-VI database~\cite{endf1, endf2}.

Different energy regimes have been
considered in the \textsc{MaGe}\ physics list to simulate the inelastic interactions of long-lived
hadrons. Each energy regime has its own specialized model.
In particular:
\begin{enumerate}
\item theory-driven quark-gluon string precompound model (QGSP) for pions, kaons
and nucleons in the high-energy region, up to 100~TeV.
\item Low-energy parametrized (LEP) model for pions and nucleons with energies between 10 and
12~GeV and for kaons below 25~GeV. Such a model is applicable below about 30~GeV and is
derived from the \textsc{Gheisha} package~\cite{gheisha}. The cross section and the final state are
determined by parametrized functions which are fitted to experimental data.
\item Bertini (BERT) cascade model, based on a re-engineering of the
\textsc{INUCL} code~\cite{inucl}, to describe nucleon and $\pi$ interactions below energies of 10~GeV.
The model includes the Bertini intra-nuclear cascade model with
excitons~\cite{bert1, bert2}, a pre-equilibrium model~\cite{gri1,gri2}, a nucleus explosion
model~\cite{breakup} and an evaporation model~\cite{weiss}.
\item data-driven high-precision (HP) model for neutrons from thermal energies up to 20~MeV,
based on the ENDF/B-VI database.
\end{enumerate}
Such an ensemble of models for hadronic inelastic interactions is shortly labeled as
QGSP\_LEP\_BERT\_HP.\\

Alternative hadronic physics lists, which differ in the models used to describe the
inelastic interactions of nucleons, are available in \textsc{MaGe}\ and can be instantiated by
messenger commands. They are:
\begin{itemize}
\item QGSP\_LEP\_BIC\_HP, which employs the Binary cascade~\cite{physics-manual} model
(instead of the Bertini cascade model) for inelastic interactions of nucleons below
10~GeV. The model handles the intra-nuclear cascade as well as the remaining
fragment, which is treated by precompound and de-excitation models. Cross sections are
parametrized using experimental data. This physics list is similar to the one used in
Ref.~\cite{Lindote:09} for the simulation of muon-induced neutrons, the only difference
being that the list of Ref.~\cite{Lindote:09} does not include the LEP bridge between
the high-energy QGSP regime and the low-energy Binary cascade.
\item QGSP\_LEP\_HP, where parametrized LEP models are used for nucleons below 10~GeV, instead
of intra-nuclear cascade models. Such a list is much faster than the default one from
the point of view of the computation time.
\item QGSC\_BERT\_HP, which employs the quark-gluon string CHIPS model (QGSC) instead of
the quark-gluon string precompound model (QGSP) to simulate high-energy inelastic
interaction of nucleons (from 10~GeV to 100~TeV). The QGSC model differs by the QGSP
by the fact that the de-excitation of the residual nucleus is handled by a CHIPS
model~\cite{chips3,chips1,chips2},
rather than the pre-compound model. The physics list does not include
the LEP parametrized model. The QGSC model is also used to produce the final state
following photo-nuclear interactions above 3.5~GeV.
\end{itemize}

Multiple alternative lists for hadronic physics based on independent models are used
in \textsc{MaGe}\ mainly for testing purposes, namely to cross-check results and evaluate
systematic uncertainties from the simulations. Crucial background sources for many
underground experiments are due to high-energy interactions of cosmic-ray muons.
Background estimates (e.g. production of secondary neutrons or long-lived
unstable nuclei) often rely on the ability of simulation codes to model high-energy
hadronic and electromagnetic showers initiated by muon interactions. Fig.~\ref{fig:germanium_mc} shows
the neutron production yield by muons in Germanium
vs. the muon energy for the four hadronic physics lists provided by \textsc{MaGe}\ (one default and
three alternative). The simulations are run with \textsc{Geant4}\ 9.0.
Simulation results are very similar (within 15\%) in the full energy range,
indicating that details in the physics modeling of high-energy hadronic interactions do not have a large
impact on some parameters of interest for underground physics experiments.
Neutron yield in Ge obtained with \textsc{MaGe}\ (about $1.5\cdot10^{-3}$ neutrons per muon
per g/cm$^{2}$ at 280~GeV)
is consistent ($< 10\%$) with the results based on the \textsc{Fluka-1999}
code~\cite{fluka}, derived by the power law parameterization of Fig.~5 of
Ref.~\cite{Araujo:05}. \\
\begin{figure*}[tbh]
\centering
\includegraphics[width=0.70\textwidth]{plots/germanium_mc.pdf}
\caption{Neutron yield from muon-induced showers in metallic germanium.
\textsc{MaGe}\ had been run with the version 9.0 of \textsc{Geant4}.}\label{fig:germanium_mc}
\end{figure*}

The \textsc{Gerda}\ experiment uses a water Cherenkov muon veto and internal
shield consisting of liquid Ar. For this reason, \textsc{MaGe}\ takes advantage of
\textsc{Geant4}\ ability to simulate optical
photons. While the default {\sc MaGe} physics list does not include
interactions of optical photons these processes can be enabled during
runtime. The underlying models encompass scintillation light emission
(possibly with different light yields for electrons,
$\alpha$-particles and nuclei), Cherenkov light emission, absorption,
boundary processes, Rayleigh scattering and wavelength shifting.
If optical photon treatment is enabled, it is
necessary to specify all relevant optical properties of interfaces and
bulk materials (refraction index, absorption length, etc.) in the
geometry definition. \\

{\sc Geant4} tracks all simulated particles down to zero range, although
various options exist to manually limit step size, track length, time-of-flight,
and other parameters. Production cuts for
$\delta$-rays and for soft bremsstrahlung photons are
expressed in spatial ranges and are internally converted into energy
thresholds for the production of soft photons and $\delta$-rays in the
corresponding material. It is necessary to find a trade-off between
accuracy and computing time in most applications. Therefore {\sc MaGe}
provides three production cut realms: DarkMatter, DoubleBeta and
CosmicRays. The DarkMatter realm is used for high-precision
simulations, especially related to background studies for dark matter
applications and surface effects: the cuts for $\gamma$-rays and $e^{\pm}$
are 5~$\mu$m and 0.5~$\mu$m, respectively, corresponding to
a $\sim$1~keV energy threshold in metallic germanium. The DoubleBeta realm
({\sc MaGe} default) is suitable for signal and background studies
related to double-beta decay, i.e. in the MeV energy-region: the
range cut for $\delta$-ray production is relaxed to 0.1~mm,
corresponding to a 100~keV threshold in metallic germanium. The
CosmicRays realm is used for the simulation of extensive
electromagnetic showers induced by cosmic-ray muons. The
cut-per-region approach is used in this setup. Sensitive regions are
defined for which the production cuts are the same as for the
DoubleBeta realm. They are more relaxed everywhere else (5~cm for
$\gamma$-rays and 1~cm for $e^{\pm}$). By avoiding the
precise tracking of particles in the inactive detector components,
computing time is saved. \\

     {\sc MaGe} includes some provisions to improve agreement between
simulation and experimental results.
For instance, simulations do not account
for inefficient conversion of germanium nuclei recoil energy to
ionization energy, also known as quenching.  {\sc MaGe} contains output classes that simulate this conversion inefficiency using the parameterizations from~\cite{ljungvall}. \\

{\sc Geant4} performs simulations on an event-by-event basis, where
each event begins with the release of a particle from a generator, and
ends when the interactions of the primary particle and its secondaries
have finished.  When long-lived radioactive decays occur, a single
{\sc Geant4} event may span many simulated years.  Output classes in
{\sc MaGe} divide {\sc Geant4} events into intervals that span user-selectable
times.  The total energy deposited during specific time intervals can
be reported.  This information can be used to simulate the
effectiveness of timing cuts at removing backgrounds.  It can also
improve agreement between results of {\sc MaGe} simulations and
experimental data.  Simulated energy deposits occurring long after the
duration of an experiment can be excluded.  Simulated energy deposits
in close succession can be summed to approximate the pile-up due to the finite time resolution
of data acquisition hardware used in an experiment.

\section{Validation of the simulation}
\label{section:validation}

The {\sc Geant4} simulation toolkit is used in various applications of
modern physics. These range from simulations in high-energy particle
physics to astrophysics and medical science. In parallel to the
development of new simulation modules the verification of the
simulation code is an important task for developers and users. Several
modules have been developed to describe the interactions of low energy
photons, electrons and hadrons with matter. These
are of particular importance for applications such as {\sc MaGe} and
are tested within the two collaborations developing the software. In
the following, the current status of {\sc MaGe} validation efforts is
summarized. \\

The description of electromagnetic interactions in the energy region
up to several MeV was tested with high-purity germanium detector
systems. A reverse electrode coaxial germanium detector was operated
directly submerged in liquid nitrogen in the so-called cryoliquid-submersion
test stand~\cite{KKPHD}.  The crystal was exposed to the radioactive
sources $^{60}$Co, $^{228}$Th and $^{152}$Eu.  The comparison between
a simple \textsc{MaGe}~simulation and recorded data is shown in
Fig.\ref{fig:cryoliquidsubmersion_mc}. The biggest deviation was found
to be approximately 12\%. Most probably the decrease of
$Data/(MC+Bkg)$ with increasing energy is due to the fact that the
inactive material between source and detector, e.g. dewar walls, is
not precisely known.\\

\begin{figure}[tbh]
\centering
\includegraphics[height=0.48\textwidth,angle=90]{plots/DataOverMcPeaksAll.pdf}
\caption{Difference between number of events in characteristic photon peaks from $\rm ^{60}Co$, $\rm ^{152}Eu$ and $\rm ^{228}Th$ sources plus background and a simple~\textsc{MaGe}~simulation. The maximal deviation was found to be approx. 12\%.}
\label{fig:cryoliquidsubmersion_mc}
\end{figure}
The most extensive verification effort was performed with an 18-fold
segmented \textsc{Gerda} prototype detector. The segmented germanium
crystal was operated in vacuum and exposed to several radioactive
sources ($^{60}$Co, $^{228}$Th,
$^{152}$Eu)~\cite{Abt:2007rg,Abt:2007rf}.  A large fraction of the
emitted $\gamma$-rays deposit energy in more than one segment. This
feature allows these events to be distinguished from those which
deposit energy in relatively small volumes, such as
$0\nu\beta\beta$\ decays.  Such segmentation-based discrimination
between single- and multiple-site interactions was compared between
experiment and simulation, with deviations found on the
5\%~level. \\

\textsc{MaGe}\ has been used to simulate the response of a variety of low-background assay detectors in use by the \textsc{Majorana}\ collaboration. In general the simulation agreed with the data to a few percent and the largest discrepancies were ascribed to uncertainties in the geometry models that have been coded in the
simulation.
The detailed shape of the Compton continuum for a HPGe detector exposed to variety of sources was also studied in~\cite{evans}. The simulation and data compared favorably in this study.\\

The description of neutron interactions with Ge is probed by comparing data
from a measurement of an AmBe source with predictions from {\sc
MaGe}. The measurements have been performed with a {\sc Clover}
detector and with the 18-fold segmented detector cited
previously~\cite{siegfried_neutron}. At an energy of several
MeV, neutrons mostly interact through elastic and inelastic scattering
as well as neutron absorption. The measured energy spectra were
studied and photon lines from neutron interactions with the germanium
detector itself and the surrounding materials were
identified~\cite{Mei:2007zd}.
A few discrepancies in the \textsc{Geant4}\ simulation were identified:
\begin{itemize}
\item the 2223.0-keV peak from H(n,$\gamma$)D appears at 2224.6 keV
in {\sc Geant4} simulations (bug report \# 955)~\cite{G4bug};
\item meta-stable nuclear states are not produced by {\sc Geant4}
as a result of neutron interactions (bug report \# 956)~\cite{G4bug};
\item neutrons do not produce internal conversion electrons in
{\sc Geant4} (bug report \# 957)~\cite{G4bug}.
\end{itemize}
The first and third bugs are related to issues in the G4NDL nuclear data files distributed with
\textsc{Geant4}, and can be solved by editing or augmenting the data library.
The third bug also involved several coding errors, corrections for which were provided by
the \textsc{MaGe}\ group. Fixing the second bug requires functionality that the \textsc{Geant4}\ Collaboration does
not plan to implement at this time.
Nevertheless, a viable solution has been developed by the \textsc{MaGe}\
group. The work-around for the meta-stable problem consists in identifying
at run-time those neutron interactions that might produce a meta-stable nucleus. In this
case, the new nucleus track is discarded, since \textsc{Geant4}\ would always generate it in
the ground state, but its position and the parent neutron energy are logged in a file.
Later on, a new simulation job is launched which uses the position information from
the previous one and where the ratio between ground and meta-stable nuclear states
is set manually to the proper value.  \\

{\sc MaGe} was also used to study and verify the
simulation of spallation neutron production and propagation. At the
CERN NA55 experiment the neutron production from a 190~GeV muon beam
incident on different targets was measured. At the SLAC electron beam
dump experiment the neutron propagation through different thicknesses of concrete was
measured. Both experiments were simulated within the {\sc MaGe}
framework. The attenuation of the neutron propagation was found to be larger in the simulation
than measured in the SLAC experiment~\cite{mike_neutron}. A method to correct the neutron
over-attenuation in \textsc{MaGe}-based simulations was implemented in \textsc{MaGe}\ and is described
in~\cite{mike_neutron}. \\

It was also found that \textsc{MaGe}/\textsc{Geant4} underestimates
the neutron production from muon interactions measured by NA55, especially
in high-Z materials, by more than a factor of two~\cite{mike_neutron}.
Results obtained in the \textsc{MaGe} simulation of NA55 have been
compared with the \textsc{Geant4}- and \textsc{Fluka}-based~\cite{fluka}
Monte Carlo simulations of the same experiment performed in~\cite{Araujo:05}, and
found to be consistent.
The disagreement between Monte Carlo simulation codes and NA55 data for muon-induced
neutron production is discussed in detail in~\cite{Araujo:05} and \cite{ilias06}.
The muon-induced neutron yield in lead has been recently
re-measured at the Boulby underground laboratory~\cite{Araujo:08}. In this case,
in which the neutron energies are much lower, the
\textsc{Geant4}-based simulation is found to over-estimate the experimental neutron yield,
which is the opposite behavior with respect to the NA55 data.
\textsc{MaGe}\ results for muon-induced neutron production by muons in liquid scintillator
have been compared to the set of experimental data reported in Ref.~\cite{Araujo:05}.
Fig.~\ref{fig:scintillator_vs_data} displays the neutron production yield by muons derived
by the \textsc{MaGe}\ simulations
(using the alternative hadronic physics lists described in Sect.~\ref{section:physics})
vs. energy, superimposed with the experimental data~\cite{mudata} reported in~\cite{Araujo:05}
and the recent data from the KamLAND experiment~\cite{kamlanddata}.
The agreement between \textsc{MaGe}\ simulations and experimental data is typically
better than a factor of two. Results are consistent with those obtained in the
recent work Ref.~\cite{Lindote:09}. \\

\begin{figure*}[tbh]
\centering
\includegraphics[width=0.70\textwidth]{plots/scintillator_vs_data.pdf}
\caption{Neutron yield from muon-induced showers in liquid scintillator (C$_{10}$H$_{22}$).
Experimental data are from the Refs.~\cite{mudata} (reported
in~\cite{Araujo:05}) and from Ref.~\cite{kamlanddata}.
Monte Carlo results
(default physics list) are well parametrized by a scaling law $E^{0.685}$
(dashed curve). \textsc{Geant4}\ 9.0 had been used for this study.}\label{fig:scintillator_vs_data}
\end{figure*}

\section{Conclusions}
\label{section:conclusions}

We presented the {\sc MaGe} framework for simulating interactions in
neutrinoless double-beta decay experiments that utilize enriched HPGe
detectors. The benefits of {\sc MaGe} can be summarized as:
\begin{itemize}
\item Reliable Monte Carlo framework based on {\sc Geant4} for
low-background, low-energy experiments.
\item Ongoing tests of the code and validation of the physics processes.
\item Flexible geometry and physics code framework that emphasizes code
reuse and verification.
\item General purpose tools like surface
and volume sampling, custom isotope decay generators, etc.
\end{itemize}

In general there is good agreement between the \textsc{Geant4}\ simulation with the MaGe physics list
and the measurements of electromagnetic interactions with average
discrepancies of the order of (5-10)\%. Several problems have been
identified in the simulation of neutron interactions.
These problems have been reported to the {\sc Geant4} collaboration
and are under investigation by the \textsc{MaGe}\ developers.
\\

We anticipate that \textsc{MaGe}\ will form the foundation of the
simulation and analysis framework of the \textsc{Gerda}\ and \textsc{Majorana}\
experiments. This type of framework and its associated event generators and physics lists is also useful for other
low-background underground experiments, such as solar,
reactor and geological neutrino experiments, direct dark matter searches,
and other neutrinoless double-beta decay search.
These experiments share many detection techniques and background
issues in common with \textsc{Gerda} and \textsc{Majorana}.

\section*{Acknowledgments}

The authors would like to express their gratitude to all members of the
\textsc{Majorana} and \textsc{Gerda} collaborations, for their continuous
feedback and useful advices.
The authors would like to thank V.~Tretyak and C.~Tull for useful
discussions, and I.~Abt for useful
discussions and her efforts to help unify the Monte Carlo of both
collaborations. This work has been supported by the
\textsc{ILIAS} integrating activity (Contract No.RII3-CT-2004-506222)
as part of the EU FP6 programme, by BMBF under grant 05CD5VT1/8,
and by DFG under grant GRK 683.
This work was also supported by
Los Alamos National Laboratory's Laboratory-Directed Research and Development Program,
and by the Office of Science of the US Department of Energy
at the University of Washington under Contract No. DE-FG02-97ER41020,
at Lawrence Berkeley National Laboratory under Contract No. DE-AC02-05CH11231, at the University of North Carolina under Contract No. DE-FG02-97ER41041,
and at Pacific Northwest National Laboratory under Contract No. DE-AC06-76RLO1830.
This research used the Parallel Distributed Systems Facility at the
National Energy Research Scientific Computing Center, which is
supported by the Office of Science of the U.S. Department of Energy
under Contract No. DE-AC02-05CH11231.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Recently, cross sections for the production of $J/\psi$-pairs were measured
at the Tevatron \cite{Abazov:2014qba} and
the LHC \cite{CMS_jpsijpsi,ATLAS_jpsijpsi,Aaij:2016bqq,Aaij:2011yc}.
There remain a number of puzzles, especially with the CMS and ATLAS data.
Here the leading order of ${\cal{O}}(\alpha_S^4)$ (see e.g. \cite{JJ_kt,Baranov:2012re}) is clearly not
sufficient. The double parton scattering (DPS) contribution was claimed to be large
or even dominant in some corners of the phase space, when the rapidity distance $\Delta y$ between
two $J/\psi$ mesons is large. However the effective cross sections $\sigma_{\rm eff}$ found
from empirical analyses are about a factor $2.5$ smaller than the usually accepted
$\sigma_{\rm eff} = 15 \, \rm{mb}$. It is an open issue at the moment whether this points to
a nonuniversality of $\sigma_{\rm eff}$ or whether there are additional single parton scattering
mechanisms which can alleviate the tension.

The production of quarkonium pairs is interesting in a broader context. Here we wish to
consider production of pairs of $\chi_c$ mesons. This process is more difficult to measure
experimentally but interesting from the theoretical point of view.
A feed down to the double $J/\psi$ channel is interesting in the context of the puzzles
mentioned above.

The single-inclusive $\chi_c$ meson production was a topic of both experimental \cite{Aaij:2013dja,Chatrchyan:2012ub,ATLAS:2014ala}
and theoretical \cite{Hagler:2000dd,Kniehl:2006sk,Likhoded:2014kfa,Baranov:2015yea,Cisek:2016uxz} studies.
The cross section for single $\chi_c$ production is rather
large. The nonrelativistic perturbative QCD is the standard theoretical approach in this
context. In leading order the gluon fusion $g^* g^* \to \chi_c(J), J=0,1,2$ is the underlying production
mechanism. The $k_T$-factorization approach provides a reasonable description of
the experimental data \cite{Baranov:2015yea,Cisek:2016uxz}.

In the present letter we shall include the production of all combinations of $\chi_c$ meson
pair production. The cross section will be calculated in $k_T$-factorization approach using
newly derived off-shell matrix elements for the $g^* g^* \to \chi_c(i) \chi_c(j)$ process.

A first evaluation of the total cross section will be given. We also show some differential
distributions.

\subsection{The $p p \to \chi_c(J_1) \chi_c(J_2) X$ reaction, formalism}

It was shown in \cite{LS2015,LLP2016} that the $\chi_c J/\psi$ pair production
is possible only at $O(\alpha_s^5)$, while forbidden at $O(\alpha_s^4)$
due to $C$ parity conservation.
In contrast, the production of $\chi_c(J_1) \chi_c(J_2)$, see Fig.\ref{fig:diagram_chicchic},
is possible already at the $O(\alpha_s^4)$ order.

Of special importance for us is the fact that $\chi_c \chi_c$ states are
produced by the crossed-channel one-gluon exchange mechanism.
This implies that the production amplitudes are flat as a function
of $g^*g^*$ center of mass energy, which implies broad distributions
in the rapidity distance $\Delta y$ between the produced
$\chi_c$-mesons.

\begin{figure}[!h]
\begin{minipage}{0.32\textwidth}
\centerline{\includegraphics[width=1.0\textwidth]{diag_double_jspi_gamma.eps}}
\end{minipage}
   \caption{A diagrammatic representation of the leading order
mechanisms for $pp \to \chi_c(J_1) \chi_c(J_2) \to
(J/\psi+\gamma) (J/\psi+\gamma)$ reaction.
 }
\label{fig:diagram_chicchic}
\end{figure}

According to our knowledge this contribution was not discussed so far
in the literature. There was, however,
some calculations for $\chi_c \chi_b$ production \cite{protvino}.

We consider the gluon-gluon fusion mechanism shown diagramatically
in Fig.\ref{fig:diagram_chicchic}.
There are altogether six possible combinations of pair production
of $\chi_c(0), \chi_c(1), \chi_c(2)$ quarkonia.

In order to calculate the subprocess amplitudes,
we first turn to the $g^* g^* \to \chi_c(J)$ vertices.

\begin{figure}[!h]
\begin{minipage}{0.32\textwidth}
\centerline{\includegraphics[width=1.0\textwidth]{chic_ver.eps}}
\end{minipage}
   \caption{A diagrammatic representation of the
$g^* g^* \to \chi_c(\lambda)$ vertex being a building block of
corresponding $g^* g^* \to \chi_c(J_1) \chi_c(J_2)$ and discussed
in this section amplitudes.
 }
\label{fig:diag_gg_chic}
\end{figure}

\subsubsection{The $g^* g^* \to \chi_c(J)$ vertices}

The $g^* g^* \to \chi_c$ vertices with off-shell gluons
(see Fig.\ref{fig:diag_gg_chic}) are building blocks of the
elementary $g^* g^* \to \chi_c(J_1) \chi_c(J_2)$ amplitudes.

Here we follow the general rules of NRQCD as explained e.g. in
\cite{Guberina:1980dc,Cho:1995ce,Cho:1995vh}.
We restrict ourselves to the color singlet contribution and can write the
amplitude for the production of the $\chi_c(J)$ meson via
the fusion of two gluons as:
 \begin{eqnarray}
 V^{ab}_{\mu \nu}(J,J_z;q_1,q_2) = 4 \pi \alpha_S {{\rm Tr}[t^a t^b] \over \sqrt{N_c}} \,
 \sqrt{{2\over M}} \sum_{S_z,L_z} \int {d^4 k \over (2 \pi)^3} \delta(k^0 - {\vec{k}^2 \over M})
 \Psi_{1,L_z}(\vec{k})
 \nonumber \\
 \times \bra{1,S_z; 1,L_z}|\ket{J,J_z} \cdot {\rm{Tr}}[A_{\mu \nu} \Pi_{1,S_z}],
\label{eq:vertex}
 \end{eqnarray}
following closely the notation of
\cite{Hagler:2000dd,Pasechnik:2007hm,Pasechnik:2009bq,Pasechnik:2009qc},
where these vertices had been calculated for external reggeized gluons.
Below we will need the amplitudes (\ref{eq:vertex}) for arbitrary off-shell momenta
of gluons, not only the multiregge kinematics as in
\cite{Hagler:2000dd,Pasechnik:2007hm,Pasechnik:2009bq,Pasechnik:2009qc}.
There is however no additional difficulty related with this. As we concentrate on
the color-singlet mechanism, the three-gluon coupling does not enter and we really deal with a QED
problem.
Consequently the amplitudes (\ref{eq:vertex}) fulfill the QED-like gauge invariance conditions:
\begin{eqnarray}
  q_1^\mu V^{ab}_{\mu \nu}(J,J_z;q_1,q_2) = 0, q_2^\nu  V^{ab}_{\mu \nu}(J,J_z;q_1,q_2) = 0.
\end{eqnarray}
The calculation proceeds as follows.

The $g^* g^* \to Q \bar Q$ amplitude is (up to factors)
\begin{eqnarray}
 A_{\mu \nu} = \gamma_\mu {\hat p_Q - \hat q_1 + m_Q \over (p_Q - q_1)^2 - m_Q^2} \gamma_\nu +
 \gamma_\nu {\hat p_Q - \hat q_2 + m_Q \over (p_Q - q_2)^2 - m_Q^2} \gamma_\mu \, .
\end{eqnarray}
We parametrize
\begin{eqnarray}
 p_Q = {P \over 2 } + k \, , \, p_{\bar Q} = {P \over 2} - k \, ,
\end{eqnarray}
In spectroscopic notation, the $\chi_c$ mesons are $^{2S+1}L_J = ^3 P_J$ states,
where $J=0,1,2$. Therefore the spinorial part of the wavefunction is
an $S=1$ spin triplet state, and the relevant projector can be
written as
\begin{eqnarray}
 \Pi_{S=1,S_z} &=& {1 \over 2 \sqrt{2}  m_Q} \, \Big({\hat P \over 2} - \hat k - m_Q \Big) \hat \epsilon(S_z)
 \Big({\hat P \over 2} + \hat k + m_Q \Big).
\end{eqnarray}
Now, for the $P$-wave states, we should expand the product $A_{\mu \nu}\hat \Pi_{S=1,S_z}$ in
(\ref{eq:vertex}) to the first order in $k$. In fact the Taylor expansion for $P$-waves starts
from the term linear in $k$:
\begin{eqnarray}
  {\rm Tr}[A_{\mu \nu} \Pi_{1,S_z}] \rightarrow  k_\alpha \cdot {\partial \over \partial k_\alpha}
  {\rm Tr}[A_{\mu \nu} \Pi_{1,S_z}]\Big|_{k=0} \, .
\end{eqnarray}
Then, the integration over relative momentum $k$ reduces to the integral
\begin{eqnarray}
\int {d^3 \vec{k} \over (2 \pi)^3} k^\alpha \Psi_{1,L_z}(\vec{k}) =
-i \sqrt{3 \over 4 \pi} \, R'(0) \cdot \epsilon^\alpha(L_z) \, .
\end{eqnarray}
Here $R'(0)$ is the derivative of the radial wavefunction at the (spatial) origin.

For convenience, we introduce
\begin{eqnarray}
T_{\mu \nu} (q_1,q_2;J,J_z) \equiv
{\sqrt{2} M \over 8} \, \sum_{S_z,L_z} \bra{1,S_z; 1,L_z}|\ket{J,J_z} \epsilon_\alpha(L_z)
\cdot {\partial \over \partial k_\alpha}  {\rm Tr}[A_{\mu \nu} \Pi_{1,S_z}]\Big|_{k=0},
\end{eqnarray}
so that our  gluon-gluon fusion vertices take the form
\begin{eqnarray}
 V^{ab}_{\mu \nu}(J,J_z;q_1,q_2) = -i \, 4 \pi \alpha_S \,
 \delta^{ab} \, {2R'(0) \over \sqrt{\pi N_c M^3}} \, \sqrt{3}
 \cdot T_{\mu \nu}(J,J_z;q_1,q_2) \, ,
\end{eqnarray}
Performing the relevant Dirac-traces, we obtain
the explicit expressions for $T_{\mu \nu}(J,J_z;q_1,q_2)$:
\begin{enumerate}
 \item scalar, $J=0$:
\begin{eqnarray}
 T_{\mu \nu}(0,0;q_1,q_2) &=& {1 \over \sqrt{3}} \,  {M^2 \over (2 q_1 \cdot q_2)^2}
 \, \Big\{ g_{\mu \nu} \Big( 6(q_1 \cdot q_2) - q_1^2 -q_2 ^2 +
 {(q_2^2 - q_1^2)^2 \over M^2}\Big) \nonumber \\
 &+& q_{1 \mu} q_{2\nu} \, 2\Big( {q_1^2 + q_2^2 \over M^2} -1 \Big) \,
 + q_{2 \mu} q_{1 \nu} \, 2\Big( {q_1^2 + q_2^2 \over M^2} -3 \Big)
 \nonumber\\
 &+& q_{1 \mu} q_{1\nu} {4 q_2^2 \over M^2} + q_{2 \mu} q_{2 \nu} {4 q_1^2 \over M^2} \Big\}
 \label{eq:scalar}
\end{eqnarray}
 \item axial vector, $J=1$:
\begin{eqnarray}
 T_{\mu \nu} (1,J_z;q_1,q_2) &=& {i \over \sqrt{2}M} {1 \over (q_1 \cdot q_2)} \Big\{
  (q_1^2 - q_2^2)
 \epsilon_{\mu \nu \alpha \beta} (q_1+ q_2)^\alpha \epsilon^{\beta}(J_z)

 \nonumber \\
 &+& {q_1^2 + q_2^2 \over  (q_1 \cdot q_2)}   (a_\mu q_{1\nu} - a_\nu q_{2\mu})

 + 2 ( a_\nu q_{1\mu} - a_\mu q_{2\nu}) \Big\} \,

 \label{eq:axial}
\end{eqnarray}
with
\begin{eqnarray}
 a_\mu =  \epsilon_{\mu \rho \alpha \beta} q_1^\rho q_2^\alpha \epsilon^{\beta}(J_z) \, .
\end{eqnarray}
\item tensor, $J=2$:
\begin{eqnarray}
 T_{\mu \nu}(2,J_z;q_1,q_2) &=& {- M^2 \over (2q_1 \cdot q_2)^2}
 \Big\{ - g_{\mu \nu} (q_2 -q_1)^\alpha (q_2 - q_1)^\beta
 \epsilon_{\alpha \beta}(J_z)  + 4 (q_1 \cdot q_2) \epsilon_{\mu \nu}(J_z)
 \nonumber \\
 &+& 2 (q_2 -q_1)^\alpha \epsilon_{\alpha \nu}(J_z) q_{2\mu}
 - 2 (q_2 - q_1)^\alpha \epsilon_{\alpha \mu}(J_z) q_{1\nu} \Big\} \, ,
 \label{eq:tensor}
\end{eqnarray}
where
$\epsilon_{\mu \nu}(J_z) = \sum_{m_1,m_2} \bra{2,J_z}|\ket{1,m_1,1,m_2} \epsilon_\mu(m_1) \epsilon_\nu(m_2)$
is the polarization tensor of the $J=2$ state.
\end{enumerate}
Notice, that
\begin{eqnarray}
 2(q_1 \cdot q_2) = M^2 - q_1^2 - q_2^2 \, ,
\end{eqnarray}
and as gluons are always spacelike $q_i^2 <0$, the denominators of eqs (\ref{eq:scalar}
,\ref{eq:axial}, \ref{eq:tensor}) are always finite.

Besides the QED-like gauge invariance condition, these amplitudes also fulfill the Bose-symmetry
\footnote{Notice that it does not mean that $T_{\mu \nu}$ is a symmetric tensor, as the
results presented in \cite{protvino} (which violate the gauge invariance condition).}
\begin{eqnarray}
  T_{\mu \nu}(J,J_z;q_1,q_2) =   T_{\nu \mu}(J,J_z;q_2,q_1) \, .
\end{eqnarray}
A comment on the $J=1$ axial vector is in order.
Here the Landau-Yang theorem forbids the decay of the $\chi_c(1)$
into $\gamma \gamma$ or $gg$, and likewise its production
through fusion of on-shell photons or gluons.
Indeed, in the limit $q_1^2 \to 0, q_2^2 \to 0$,
we have
\begin{eqnarray}
 T_{\mu \nu} (1,J_z;q_1,q_2) \propto a_\nu q_{1\mu} - a_\mu q_{2\nu}  \, ,
\end{eqnarray}
which vanishes, when contracted with the polarization vectors
of on-shell photons/gluons
\begin{eqnarray}
 \epsilon_1^\mu \epsilon_2^\nu ( a_\nu q_{1\mu} - a_\mu q_{2\nu} ) = 0,
\end{eqnarray}
as required by the Landau-Yang theorem.
\subsubsection{The $g^* g^* \to \chi_c(J_1) \chi_c(J_2)$ amplitudes}
Now we wish to discuss the elementary
$g^* g^* \to \chi_c(J_1) \chi_c(J_2)$ amplitudes, which can be obtained from the
building blocks discussed above.

In all cases there are two diagrams ($t$ (left) and $u$ (right)
in Fig.\ref{fig:chicchic_amplitudes}).

\begin{figure}[!h]
\includegraphics[width=4cm]{diag_chic_chic_t.eps}
\includegraphics[width=4cm]{diag_chic_chic_u.eps}
   \caption{A diagrammatic representation of the
generic $g^* g^* \to \chi_c(J_1) \chi_c(J_2)$ $t$-channel (left)
and $u$-channel (right) amplitudes.
 }
\label{fig:chicchic_amplitudes}
\end{figure}

We can write the Feynman amplitudes corresponding to these diagrams as
\begin{eqnarray}
 {\cal M}^{ab}_{\mu \nu}(J_1,J_{1z},J_2,J_{2z}) &=&  V^{ac}_{\mu \alpha}(J_1,J_{1z};q_1,p_1-q_1)
 {-g^{\alpha \beta} \delta_{cd} \over \hat t} V^{db}_{\beta \nu}(J_2,J_{2z};p_2-q_2,q_2) \nonumber\\
 &+& V^{ac}_{\mu \alpha}(J_2,J_{2z};q_1,p_2-q_1)
 {-g^{\alpha \beta} \delta_{cd} \over \hat u} V^{db}_{\beta \nu}(J_1,J_{1z};p_1-q_2,q_1) \, ,
\label{eq:chi_chi}
 \end{eqnarray}
where $\hat t = (p_1 - q_1)^2 = (p_2-q_2)^2 \, , \, \hat u = (p_2 - q_1)^2 = (p_1 - q_2)^2$.
These amplitudes are infrared finite and gauge invariant.

To obtain the $k_T$-factorization amplitude one should contract (\ref{eq:chi_chi}) with the
polarization vectors of off-shell gluons
\begin{eqnarray}
 e_{1 \mu} = {q_{1T \mu} \over |\vec{q}_{1T}| } \, , \,
 e_{2 \nu} = {q_{2T \nu} \over |\vec{q}_{2T}| } \, .
\end{eqnarray}
Because of the QED-like Ward identities of the gluon fusion vertices,
these polarization vectors are equivalent to the more common Gribov's polarizations
$n^+_\mu, n^-_\nu$, for incoming gluons
in the high-energy kinematics $q_{1\mu} = q_1^+ n^+_\mu + q_{1T \mu}$,
 $q_{2\nu} = q_1^+ n^-_\nu + q_{2T \nu}$.

In the nonrelativistic QCD approach the cross section for
$\chi_c$ pair production is proportional to $|R'(0)|^4$.
The result is therefore extremely sensitive to the precise value
of the wave function derivative at the origin.
In our opinion the best estimate of the parameter can be obtained from:
\begin{equation}
\Gamma(\chi_c(0^+) \to \gamma \gamma) =
{27 e_c^4 \alpha_{\rm em} \over m_c^4} |R'(0)|^2 \; .
\end{equation}
From the experimental value of the diphoton decay width \cite{PDG}
one obtains for the $\chi_c$ P-wave function squared
\begin{equation}
|R'(0)|^2 = 0.042 \; {\rm GeV}^5.
\end{equation}

In the following the $\chi_c(J_1) \chi_c(J_2)$ cross section is
calculated within the $k_T$-factorization approach including
off-shell matrix elements for the $g^* g^* \to \chi_c(J_1) \chi_c(J_2)$
subprocess and modern unintegrated gluon distributions.

It is well known that about 30 \% of prompt single $J/\psi$ production
originates from radiative decays $\chi_c \to J/\psi + \gamma$
with branching fractions:
${\rm{Br}}(\chi_c(0) \to J/\psi \gamma) =  1.26 \pm 0.06 \% $,
${\rm{Br}}(\chi_c(1) \to J/\psi \gamma) = 33.9 \pm 1.2 \%$,
${\rm{Br}}(\chi_c(2) \to J/\psi \gamma) = 19.2 \pm 0.7 \%$ \cite{PDG}.
Obviously, regarding feed down into the $J/\psi J/\psi$ channel only
the $\chi_c(1) \chi_c(1)$,
$\chi_c(1) \chi_c(2)$ and $\chi_c(2) \chi_c(2)$ states could give potentially
important contributions. The details depend, however, on corresponding matrix elements
and cross sections for the $\chi_c(J_1) \chi_c(J_2)$ production.

The cross section for $p p \to \chi_c(J_1) \chi_c(J_2)$
is calculated in the $k_T$-factorization approach.
The corresponding differential cross section for the production
of $\chi_c(i) \chi_c(j)$ states, where $i$ and $j$ run through $0,1,2$
can be written as:
\begin{eqnarray}
&&\frac{d \sigma(p p \to \chi_c(i) \chi_c(j) X)}
{d y_1 d y_2 d^2 {\vec p}_{1T} d^2 \vec{p}_{2T}}
 =
\frac{1}{16 \pi^2 (x_1 x_2 s)^2} {1 \over 1+ \delta_{ij}} \int \frac{d^2 \vec{q}_{1T}}{\pi \vec{q}_{1T}^2}
\frac{d^2 \vec{q}_{2T}}{\pi \vec{q}_{2T}^2}
\overline{|{\cal M}_{g^{*} g^{*} \rightarrow \chi_c(i) \chi_c(j)}^{\rm{off-shell}}|^2}
\nonumber \\
&& \times \;\;
\delta^{(2)} \left( \vec{q}_{1T} + \vec{q}_{2T} - \vec{p}_{1T} - \vec{p}_{2T} \right)
{\cal F}(x_1,\vec{q}_{1T}^2,\mu_{F}^{2})
{\cal F}(x_2,\vec{q}_{2T}^2,\mu_{F}^{2}) \; .
\label{kt_fact_gg_jpsijpsi}
\end{eqnarray}
The unintegrated gluon distribution ${\cal F}(x_1,\vec{q}_{1T}^2,\mu_{F}^{2})$ is related to the
collinear one through
\begin{eqnarray}
x g(x,\mu_F^2) = \int^{\mu^2_F} {d \vec{q}^2_T \over\vec{q}^2_T  } \, {\cal F}(x,\vec{q}_{1T}^2,\mu_{F}^{2})\, ,
\end{eqnarray}
and the off-shell matrix element is obtained as
\begin{eqnarray}
 \overline{|{\cal M}_{g^{*} g^{*} \rightarrow \chi_c \chi_c}^{\rm{off-shell}}|^2}
 ={1 \over (N_c^2-1)^2} \sum_{a,b,J_{1z},J_{2z}}
 |e_{1\mu} e_{2\nu} {\cal{M}}^{ab}_{\mu \nu}(J_1,J_{1z},J_2,J_{2z})|^2
 \; .
\end{eqnarray}
The longitudinal momentum fractions $x_1$ and $x_2$ are calculated from $\chi_c$'s
transverse masses $m_{Ti} = \sqrt{m_c^2 + \vec{p}_{iT}^2}$ and rapidities:
\begin{eqnarray}
x_1 &=& {m_{T1} \over \sqrt{s} } e^{y_1} +  {m_{T2} \over \sqrt{s} } e^{y_2} \, , \nonumber \\
x_2 &=& {m_{T1} \over \sqrt{s} } e^{-y_1} +  {m_{T2} \over \sqrt{s} } e^{-y_2}  \, .
\end{eqnarray}
\subsection{Results for $\chi_c(J_1) \chi_c(J_2)$ production}

We start presentation of our results by showing integrated cross
sections.
As an example in Table \ref{table:cross_sections}
we show cross section in a broad range of $\chi_c$ rapidities.
We used an unintegrated gluon distribution constructed from the KMR
prescription \cite{Kimber:2001sc} based on the MSTW2008 collinear NLO
gluon distribution \cite{MSTW08}.
For the renormalization scales $\mu^2_{r1},\mu^2_{r2}$  of the running coupling
and factorization scales $\mu^2_{F1},\mu^2_{F2}$ entering the unintegrated
gluon distribution, we choose
\begin{eqnarray}
\mu_{r1}^2 &=& \mu^2_{F1} = {\rm max}\{m_c^2, |\vec{q_{1T}}|^2\} \, ,\nonumber \\
\mu_{r2}^2 &=& \mu^2_{F2} = {\rm max}\{m_c^2, |\vec{q_{2T}}|^2\} \ ,
\end{eqnarray}
where these scales refer to the running coupling/gluon distribution coupling
to gluon $q_1$ or $q_2$ respectively.
We refrain from a detailed analysis of dependence on the factorization scale,
the distributions shown below simply serve to get an impression of the
salient features of the production mechanism.
A more detailed analysis, including theoretical errors will be given in a future work
\cite{paperJPsi}, where we will address the feeddown into the $J/\psi J/\psi$ channel.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
             &  $\chi_c(0)$  &  $\chi_c(1)$  &  $\chi_c(2)$  \\
\hline
$\chi_c(0)$  &  1.32       &   1.71         &   4.24       \\
$\chi_c(1)$  &  ....       &   0.84         &   2.88       \\
$\chi_c(2)$  &  ....       &   ....         &   3.45     \\
\hline
\end{tabular}
\end{center}
\caption{Cross sections in nb for production of different combinations
	of $\chi_c(J_1) \chi_c(J_2)$ dimeson states
	for -8 $< y_1, y_2 <$ 8 at $\sqrt{s}$ = 8 TeV.
	The numbers are obtained in the $k_T$-factorization approach.
	We used an unintegrated gluon distribution constructed from the KMR
	prescription \cite{Kimber:2001sc} based on the MSTW2008 collinear NLO
	gluon distribution \cite{MSTW08}.
	In all cases the gauge invariant matrix elements discussed in the
	present paper were used.}
\label{table:cross_sections}
\end{table}

There are six independent cross sections related to the different spin
combinations (see Table \ref{table:cross_sections}).
We see that the cross sections for different spin combinations
are of the same order of magnitude.

In Fig.\ref{fig:dsig_dy} we show rapidity distributions for
$\chi_c$ mesons for different pair combinations.
In the left panel we show:
$\chi_c(0) \chi_c(0)$ (solid line), $\chi_c(1) \chi_c(1)$ (dahed line)
and $\chi_c(2) \chi_c(2)$ (dotted line).
In the other panels we show distributions for:
$\chi_c(0) \chi_c(1)$ (solid line), $\chi_c(0) \chi_c(2)$ (dashed line)
and $\chi_c(1) \chi_c(2)$ (dotted line).
In the upper keft panel the distribution of the first listed quarkonium is shown,
while the distributions of the second listed quarkonium
are shown in the lower panel. Evidently for the nonidentical quarkonia the
distribution of the first and second meson are not the same.

\begin{figure}[!h]

	\includegraphics[width=.4\textwidth]{dsig_dy1_chic_ii.eps}
	\includegraphics[width=.4\textwidth]{dsig_dy1_chic_ij.eps}
	\includegraphics[width=.4\textwidth]{dsig_dy2_chic_ij.eps}
	\caption{Rapidity distributions of quarkonia for different spin
		combinations.}
	\label{fig:dsig_dy}
\end{figure}

In Fig.\ref{fig:dsig_dpt} we show similar distributions in quarkonia
transverse momenta. The distributions for $\chi_c(1)$ quarkonia are
less steep than those for the other mesons. This may have important
consequences for large transverse momenta, also for $J/\psi$ pair production
(CDF, ATLAS, CMS), but goes beyond the scope of the present letter.

\begin{figure}[!h]
\includegraphics[width=.4\textwidth]{dsig_dpt1_chic_ii.eps}
\includegraphics[width=.4\textwidth]{dsig_dpt1_chic_ij.eps}
\includegraphics[width=.4\textwidth]{dsig_dpt2_chic_ij.eps}
\caption{Transverse momentum distributions of quarkonia for different
spin combinations.}
\label{fig:dsig_dpt}
\end{figure}

The exchange of gluons leads to broad distributions in the difference
of rapidities $\Delta y$ of the two quarkonia, as shown in
Fig.\ref{fig:dsig_dydiff}. All final states have in common also a rather
deep dip at $\Delta y =0$. Therefore the $\chi_c$ pair production will be
potentially important rather for experimental setups that cover a large range in
rapidities.

\begin{figure}[!h]
\includegraphics[width=.4\textwidth]{dsig_deltay_chic_ii.eps}
\includegraphics[width=.4\textwidth]{dsig_deltay_chic_ij.eps}
\caption{Distributions in the rapidity separation between $\chi_c$'s
for different spin combinations.}
\label{fig:dsig_dydiff}
\end{figure}

In calculations based on collinear gluon distributions, the two $\chi_c$
mesons are produced back-to-back at the lowest order.
This is not so in the $k_T$-factorization approach discussed here.
In Fig.\ref{fig:dsig_dptsum} we show distributions of the transverse momentum
of the meson pair, $p_{T,sum}$.
The distribution for the $\chi_c(1) \chi_c(1)$ extends to large pair
transverse momenta, which is related to the corresponding vertex structure.

\begin{figure}[!h]
\includegraphics[width=.4\textwidth]{dsig_dptsum_chic_ii.eps}
\includegraphics[width=.4\textwidth]{dsig_dptsum_chic_ij.eps}
\caption{Distributions in the transverse momentum of quarkonium pairs
for different spin combinations.}
\label{fig:dsig_dptsum}
\end{figure}

The $\chi_c$ mesons radiatively decay into $J/\psi$ mesons.
The double feed down leads to a new contribution to the $J/\psi J/\psi$
channel. The direct $J/\psi J/\psi$ contribution is more than order
of magnitude larger than the feed-down contribution.
However, the $\chi_c \chi_c$ contribution has its own specificity.
In Fig.\ref{fig:dsig_dydiff_jpsijpsi} we show distribution in rapidity
difference for all $\chi_c \chi_c$ contributions weighted
by branching fractions into $J/\psi$ channel (solid line)
compared to the standard direct $J/\psi J/\psi$ contribution (dashed
line). At large rapidity difference the feed-down contribution dominates
over the contribution of the standard mechanism.
Here we assumed that the $J/\psi$'s from the decay will be collinear to
their parent $\chi_c$'s. How important is the feed-down contribution for
different experimental situations will be discussed elsewhere \cite{paperJPsi}.

\begin{figure}[!h]
\includegraphics[width=8cm]{dsig_deltay_jpsijpsi.eps}
\caption{Distributions in the rapidity difference
between two $J/\psi$ (dashed line) and for the sum
over all $\chi_c \chi_c$ combinations multiplied by
combined branching fractions.}
\label{fig:dsig_dydiff_jpsijpsi}
\end{figure}

\section{Conclusions}

We have made a first exploratory study of $\chi_c$ pair production
in proton-proton collisions.
The $g^* g^* \to \chi_c(J_i) \chi_c(J_j)$ amplitudes for off-shell gluons
and different spin combinations $J_i, J_j = 0,1,2$ are calculated based
on $g^* g^* \to \chi_c(J)$ verticies calculated within the color-singlet
nonrelativistic pQCD approach. In this approach the vertices are proportional to
the derivative of the spatial wave function at the origin $|R'(0)|$.
The value of this quantity can be obtained from models of the quarkonia
states. Here it has been obtained from the $\chi_c(0) \to \gamma \gamma$
branching fraction which was measured experimentally.

We have performed calculations within the $k_T$-factorization approach
for the $p p \to \chi_c \chi_c X$ process at $\sqrt{s} = 8 \, \rm{TeV}$ using
Kimber-Martin-Ryskin \cite{Kimber:2001sc}
type unintegrated gluon distribution based on the MSTW2008
\cite{MSTW08} collinear gluons.

We have found that the cross sections for different combinations
of $\chi_c$ quarkonia are of a similar size.
The integrated cross sections for different channels are of the order
of a few nb.
This is of the same order of magnitude as the cross section for
$J/\psi$ pair production.
This means that a feedown from the double $\chi_c$ decays
$\chi_c \to J/\psi \gamma$ leads to extra nonnegligible contribution
which has to be included in the total prompt production of two $J/\psi$
mesons. Due to specific branching fractions the $\chi_c(1) \chi_c(1)$,
$\chi_c(1) \chi_c(2)$ and $\chi_c(2) \chi_c(2)$ channels are
the dominant ones. The other three contributions can be safly neglected.

The $\chi_c \chi_c$ contribution to the $J/\psi J/\psi$ final state
is interesting but goes beyond the scope of the present analysis and
will be studied in detail in future dedicated analyses.

The salient feature of the $t$ and $u$-channel gluon exchange mechanism
are the broad distributions in rapidity difference $\Delta y$ between $\chi_c$ mesons.
This is to be contrasted with the narrow $\Delta y$ distribution of $J/\psi$ pairs
at leading order.
A feed-down from double $\chi_c$ production to the double $J/\psi$
channel is therefore expected to be important at large $\Delta y$ and may
mimic the kinematical behaviour of double parton scattering mechanisms.

\vspace{1cm}

{\bf Acknowledgments}

We would like to thank Sergey Baranov for discussions and comments on the
manuscript.
This study was partially supported by the Polish National Science Center
grant DEC-2014/15/B/ST2/02528 and by the Center for Innovation and
Transfer of Natural Sciences and Engineering Knowledge in
Rzesz{\'o}w.