#  A Peters cycle at the end of the cosmic ray spectrum?

Marco Stein Muzio  0000-0003-4615-5529<EMAIL_ADDRESS>Department of Astronomy
and Astrophysics, Pennsylvania State University, University Park, PA 16802,
USA  Department of Physics, Pennsylvania State University, University Park, PA
16802, USA  Institute of Gravitation and the Cosmos, Center for Multi-
Messenger Astrophysics, Pennsylvania State University, University Park, PA
16802, USA  Luis A. Anchordoqui  0000-0003-1463-7136 Department of Physics and
Astronomy, Lehman College, City University of New York, NY 10468, USA
Department of Physics, Graduate Center, City University of New York, NY 10016,
USA  Department of Astrophysics, American Museum of Natural History, NY 10024,
USA  Michael Unger  0000-0002-7651-0272 Institut für Astroteilchenphysik,
Karlsruher Institut für Technologie, Karlsruhe 76344, Germany  Institutt for
fysikk, Norwegian University of Science and Technology (NTNU), Trondheim,
Norway

###### Abstract

We investigate the degree to which current ultrahigh energy cosmic ray
observations above the ankle support a common maximum rigidity for all nuclei,
often called a Peters cycle, over alternative scenarios for the cosmic ray
spectra escaping sources. We show that a Peters cycle is not generally
supported by the data when compared with these alternatives. We explore the
observational signatures of non-Peters cycle scenarios, and the opportunities
to explore both ultrahigh energy cosmic ray source conditions, as well as,
physics beyond the Standard model they present.

##  I Introduction

One of the most challenging open questions regarding the origin of ultrahigh
energy cosmic rays (UHECRs) deals with the relative maximum energies of the
spectra escaping the source for different nuclei. A commonly-used simplifying
assumption is that the cosmic ray source spectra trace a Peters cycle, in
which the maximum cosmic ray energy scales linearly with $Z$, i.e., with the
charge of the UHECR in units of the proton charge [1] (see e.g. [2, 3, 4] for
fits of UHECR data with a Peters cycle at the source). A Peters cycle arises
naturally for acceleration processes which depend on magnetic fields to
confine cosmic rays to the acceleration region. This limits the maximum
rigidity of the accelerator to $R_{\mathrm{max}}\lesssim BL$, where $B$ and
$L$ are the magnetic field strength and size of the accelerating region. This
condition, often called the Hillas criterion [5], results in a common maximum
rigidity among all nuclei which are accelerated.

Current UHECR data, though, is not compatible with a pure Peters cycle from a
single source population across the entire spectrum [6, 7] above $10^{18}$ eV.
In particular, current spectrum and composition data require at least one of
two possibilities in order to be reconciled. The first possibility is that the
UHECR flux below the ankle is dominantly produced by a second class of
extragalactic sources (see [7]). However, the number of different source
populations and their relative variance cannot be too large [8]. Another
possibility is that an alternative scaling of maximum energies is required to
explain the full UHECR spectrum and composition data (see, e.g., [9, 10, 11,
12]). Such alternative scalings are a natural assumption for models which
consider the possibility that UHECRs escape their source only after suffering
significant energy losses. A specific illustration of this in the context of
gamma-ray bursts (GRBs) is presented in [13].

However, since the CR escape time naturally decreases with rigidity, it is
possible that the highest energy CRs suffer minimal energy losses and escape
their sources with a Peters cycle intact. Therefore, here we consider whether
the UHECR data above the ankle, in particular above $10^{18.8}$ eV, can still
be reconciled with a pure Peters cycle, rather than one modified by energy
losses.

Energy losses are typically parametrized in terms of $Z$ and the UHECR baryon
number $A$. For example, the energy loss rate of synchrotron and curvature
radiative processes scales as $Z^{4}/A^{2}$ and $Z^{2}$, respectively, see
e.g. [14], leading to different scalings for the maximum energy [15, 16, 17].
For instance, when considering a diffusive acceleration mechanism, synchrotron
radiation limits the cosmic ray maximum energy so that it scales as $E_{\rm
max}\propto A^{4}/Z^{4}$. On the other hand, for one-shot acceleration
processes, synchrotron radiation constrains the maximum energy to scale as
$E_{\rm max}\propto A^{2}/Z^{3/2}$, whereas curvature radiation yields $E_{\rm
max}\propto A/Z^{1/4}$. Photonuclear interactions with the thermal radiation
fields and hadronic interactions with the ambient gas are driven by the
scattering cross section, which scales with $A$. Even consideration of a
single particle species at injection can provide a complex multi-species
source spectra after propagation through the source environment [9, 12, 18].

The maximum energy of cosmic rays may even be independent of both $Z$ and $A$.
A specific example of this type of framework is the dark dimension scenario,
which naturally addresses the cosmological hierarchy problem by adding one
mesoscopic dimension of micron scale [19]. Since within this scenario physics
becomes strongly coupled to gravity around $10^{10}~{}{\rm GeV}$, it was
recently conjectured that new universal energy losses deep-rooted within the
dark dimension could control the cosmic ray maximum energy [19, 20, 21].

In this paper we take a pragmatic approach to investigate whether existing
data favor any of the above mentioned scenarios over a pure Peters cycle above
the ankle. Using data from the Pierre Auger Observatory, we carry out a
statistical analysis to study the degree to which the observed spectrum and
nuclear composition constrain the source emission spectra.

The layout of the paper is as follows. In Sec. II we provide an overview of
the data and we introduce the framework for a data-driven analysis to model
the maximum cosmic ray energy and constrain the source spectra. In Sec. III we
discuss the results for benchmark models. We find that a Peters cycle is not
the most favorable model. In Sec. IV we provide an astrophysical
interpretation of our results and we discuss possible multimessenger
connections and implications. Our conclusions are collected in Sec. V.

##  II Model

We adopt the working assumption that the all source spectra can be described
by

$\displaystyle\frac{dN(A)}{dE}\varpropto E^{-\gamma}\,e^{-E/E^{A}_{{\rm
max}}}\,,$ (1)

with

$E_{{\rm max}}^{A}=E_{0}\,Z^{\alpha}A^{\beta},$ (2)

where $E_{0}$ is the proton maximum energy of the system, and the spectral
index is common to all nuclei. Note that the case with $\alpha=1$ and
$\beta=0$ corresponds to the Peters cycle, whereas the case with $\alpha=0$
and $\beta=0$ stands for the case in which the spectrum is dominated by some
universal energy loss which is independent of the cosmic ray Lorentz boost.

With this model for the spectra escaping the source we consider five mass
groups representing $p$, He, CNO, Si, and Fe, whose relative abundances are
adjusted to obtain the best-fit. The mass groups are propagated through the
cosmic microwave background (CMB) and extragalactic background light (EBL)
using propagation matrices built from CRPropa3 [22]. The observed spectrum of
CRs is fit to the spectrum [23] and composition [24] of the Pierre Auger
Observatory (Auger) by minimizing the value of a combined $\chi^{2}$,

$\displaystyle\chi^{2}=\displaystyle\sum_{i}\frac{(J_{i}-\hat{J}_{i})^{2}}{\sigma_{J,i}^{2}}+\sum_{i}\frac{(\mu_{i}-\hat{\mu}_{i})^{2}}{\sigma_{\mu,i}^{2}}+\sum_{i}\frac{(V_{i}-\hat{V}_{i})^{2}}{\sigma_{V,i}^{2}}~{},$
(3)

where $J$ is the UHECR flux and $\mu$ and $V$ are the mean and variance of the
distribution of the logarithm of cosmic-ray mass number, $\ln{A}$. Quantities
with a hat are the model prediction and the subscript $i$ denotes the energy.
In particular, to interpret the air shower data provided by Auger we must
adopt a hadronic interaction model. We consider two hadronic interaction
models, Sibyll2.3d [25] and Epos-LHC [26], to interpret the depth of shower
maximum, $X_{\mathrm{max}}$, data in terms of $\ln{A}$.

We perform a fit to the Auger data [23, 24] above $10^{18.8}$ eV so as to
directly address the question of whether a pure Peters cycle is preferred
above the ankle. At lower energies, we assume either an additional source
population contributes to the spectrum or that energy loss processes have
significantly distorted the original Peters cycle. Given our free parameters
(the proton maximum energy $E_{0}$, spectral index $\gamma$, and $4$
parameters controlling $5$ mass group fractions) this choice of fit range
leaves $N_{\mathrm{dof}}=29$. Here we focus on a benchmark set of systematic
data shifts, which provide the best-fit to the Auger data overall: shifting
the energy scale by $\mathrm{dlg}E=+0.1$ and shifting the $\langle
X_{\mathrm{max}}\rangle$ by $-1\sigma_{X}$, since these gave the best overall
fit of the shifts explored. We consider the sensitivity of our results to
other systematic shifts of the data in Appendix A.

##  III Results

After optimizing the model parameters for a particular $(\alpha,\beta)$
combination we calculate the relative goodness-of-fit, $N_{\sigma}$, compared
to a Peters cycle in units of sigma. To calculate $N_{\sigma}$ we apply Wilks’
theorem to convert the $\Delta\chi^{2}$ between an alternative scenario and a
Peters cycle to a p-value, given that alternative scenarios have $\Delta
N_{\mathrm{dof}}=2$. We calculate $\Delta\chi^{2}$ as

$\displaystyle\Delta\chi^{2}\equiv
S^{-1}\sqrt{\left\lvert\chi_{\alpha,\beta}^{2}-\chi^{2}_{\mathrm{Peters}}\right\rvert}$
(4)

where
$S=\smash[b]{\left(\min(\chi_{\alpha,\beta}^{2},\chi^{2}_{\mathrm{Peters}})/N_{\mathrm{dof}})\right)}^{1/2}$
is a scale factor introduced to enlarge the uncertainties to account for a
$\chi^{2}_{\mathrm{min}}/N_{\mathrm{dof}}>1$ [27]. Additionally, we assign
$N_{\sigma}$ a sign equal to
$\operatorname{sgn}\left(\chi_{\alpha,\beta}^{2}-\chi^{2}_{\mathrm{Peters}}\right)$
to encode whether the fit has improved or worsened compared to a Peters cycle.
In other words, negative values of $N_{\sigma}$ represent the statistical
significance at which one can reject the null hypothesis of a pure Peters
cycle in favor of the alternative scenario $(\alpha,\beta)$.

(a)  (b)

Figure 1: Change in quality of fit to the UHECR spectrum and composition
relative to a Peters cycle. We consider (a) Sibyll2.3d and (b) Epos-LHC with
data shifted by $\mathrm{dlg}E=0.1$ and $-1\sigma_{X}$. The family of
scenarios with $\alpha+\beta=0$ are indicated by the white dashed line. The
Peters cycle (PC) and a number of alternative scenarios (green dots) are
highlighted: a photodisintegration-limited spectrum (PD), a synchrotron-
limited diffusion accelerated spectrum (SDA), a synchrotron-limited one-shot
accelerated spectrum (S1A), a curvature radiation-limited one-shot accelerated
spectrum (C1A), and a universal energy loss spectrum (UEL).

Figure 1 shows $N_{\sigma}$ for Sibyll2.3d and Epos-LHC using our benchmark
systematic shifts for the data. For reference, the Peters cycle (PC) and a
number of alternative scenarios are highlighted: a photodisintegration-limited
spectrum (PD), a synchrotron-limited diffusion accelerated spectrum (SDA), a
synchrotron-limited one-shot accelerated spectrum (S1A), a curvature
radiation-limited one-shot accelerated spectrum (C1A), and a universal energy
loss spectrum (UEL).

The value of $\alpha$ and $\beta$ change the relative position of nuclei in
the spectrum emerging from the source (as illustrated in Appendix B). This in
turn changes the relative position of mass groups at Earth which affects the
quality of fit to the CR composition data. Therefore, one expects that points
in the $\alpha-\beta$ plane with similar ratios of maximum energy between
nuclei to produce similar quality fits. In practice this is realized due to
the approximate degeneracy between $A$ and $Z$, since $A\simeq 2Z$ for stable
nuclei with the exception of protons where $A=Z$ (though these fall below our
fit range in some cases, including the standard Peters cycle scenario). In
particular, the ratio between maximum energies for two nuclei $A$ and
$A^{\prime}$ will be constant along lines of constant $\alpha+\beta$, since:

$\displaystyle\frac{E_{\mathrm{max}}^{A}}{E_{\mathrm{max}}^{A^{\prime}}}=\left(\frac{A}{A^{\prime}}\right)^{\alpha+\beta}~{}.$
(5)

By contrast, the ratio between the maximum energies of a nucleus $A$ and
protons will be constant along lines of constant
$(1-\log_{A}{2})\alpha+\beta$, since:

$\displaystyle\frac{E_{\mathrm{max}}^{A}}{E_{\mathrm{max}}^{p}}=2^{-\alpha}A^{\alpha+\beta}~{}.$
(6)

Equations (5) and (6) have two important consequences. First, in the case that
the escaping proton flux is not significant relative to other mass groups,
different acceleration scenarios fall into families which have a common
$\alpha+\beta$. Within these families, one can expect that the observed UHECR
spectrum and composition are nearly indistinguishable despite the
fundamentally different processes responsible for them. Second, in the case
that a substantial proton population escapes the source, the proton spectrum
can peak at arbitrarily high energies compared to the peak energies of nuclei
irrespective of the sign of $\alpha+\beta$. As an example, this means that
even for scenarios in the Peters cycle family $\alpha+\beta=1$, nuclei remain
ordered in terms of their peak energy but protons can peak below helium, above
iron, or between the two.

As can be seen from Fig. 1, alternative scenarios are favored over a simple
Peters cycle. For Epos-LHC alternative scenarios can achieve $N_{\sigma}<-4$,
indicating that these scenarios could be used to reject the Peters cycle
hypothesis with strong statistical significance. Notably the region with the
most negative value of $N_{\sigma}$ has a slope deviating from $-1$,
indicating that the position of protons relative to heavier nuclei is
responsible for the improvement in fit quality. In particular, this region
satisfies $\beta\simeq 0.4-0.8\alpha$ and we can infer that it is the ratio of
the maximum proton energy to the maximum energy of $A\simeq 32$ which drives
the fit.

It is clear from Fig. 1 that the global minimum is outside the plotted range
(which was driven by the range of $\alpha$ and $\beta$ values of the
alternative scenarios we considered). To explore how far outside of this range
the global minimum lies, we performed a 1-D scan along the line giving the
best-fit for Epos-LHC: $\beta\simeq 0.4-0.8\alpha$. The results of this scan
are shown in Fig. 2. Independent of the hadronic interaction model, the best-
fit is found for $\alpha\simeq 6.75$ and $\beta\simeq-5$. Currently, we are
unaware of any scenarios which could produce such an exotic dependence on mass
and charge.

Figure 2: Reduced $\chi^{2}$ along the line $\beta=0.4-0.8\alpha$ for
Sibyll2.3d and Epos-LHC. In both cases the minimum appears around
$\alpha\simeq 6.75$ and $\beta\simeq=-5$.

Figure 3 shows the best-fit spectrum and composition-related observables under
the Peters cycle and alternative scenarios highlighted in Fig. 1. The global
best-fit (BF) scenario found from the 1-D scan discussed above is also plotted
in Fig. 3. While predictions for the UHECR spectrum are very similar between
all the scenarios considered, there are differences in their predictions for
the UHECR composition. It is worthwhile noting that the models presented here
explain nearly all the flux in the two data points below our fit range,
requiring a fine-tuned transition to a second, extragalactic source
population. This may hint that the UHECR spectrum is best explained by a
energy-loss-modified spectrum, because the nucleons created during energy
losses in the source environment can populate the flux just below the ankle,
see [9].

As is hinted at by Fig. 1 and (5), the best-fit spectrum and composition for
models with similar $\alpha+\beta$ are nearly indistinguishable from each
other if no significant proton component exists in the escaping spectrum. This
explains the stark similarity both between a Peters cycle and
photodisintegration-limited scenario ($\alpha+\beta=1$), as well as, between a
synchrotron-limited diffusion accelerated scenario and a universal energy loss
scenario ($\alpha+\beta=0$). Owing to its having $\alpha+\beta=0.75$, the
curvature radiation-limited one-shot accelerated scenario falls between the
Peters cycle/photodisintegration-limited scenario and the synchrotron-limited
one-shot accelerated scenario ($\alpha+\beta=0.5$). While some of these
alternative scenarios give slightly better or worse fits compared to the
Peters cycle, Fig. 3 makes clear the difficulty in distinguishing between
these scenarios in a statistically significant way.

This pattern is clearly broken by the best-fit scenario, which is most similar
to the $\alpha+\beta=0$ family of scenarios despite its having
$\alpha+\beta\simeq 1.75$. This is due to its significant escaping proton flux
and large value of $\alpha$, so that this proton flux is at high energies
relative to nuclei. This can be seen explicitly in Figs. 6f and 7f. While this
feature may not be easily accessible through the high-level information
provided by the UHECR spectrum, $\langle X_{\mathrm{max}}\rangle$, and
$\sigma(X_{\mathrm{max}})$ data, its other signatures will be discussed in the
following Section.

(a) Fit using Sibyll2.3d for the relation between mass and shower maximum,
$X_{\text{max}}$.

(b) Fit using Epos-LHC for the relation between mass and shower maximum,
$X_{\text{max}}$.

Figure 3: Best-fit spectra (top panels) and composition-related observables
(bottom panels) assuming a Peters cycle and other alternative scenarios
(colored lines). Data points are the Auger 2021 spectrum [23] and 2019
composition data [24] shifted by $\mathrm{dlg}E=+0.1$ and $-1\sigma_{X}$.
Predictions for $\langle X_{\mathrm{max}}\rangle$ and
$\sigma(X_{\mathrm{max}})$ are made assuming (a) Sibyll2.3d and (b) Epos-LHC.
Precitions for pure proton and iron compositions for each hadronic interaction
model are shown for reference (solid gray lines). Models were fit to the data
above $10^{18.8}$ eV (solid points). Data points below the fit range are shown
for reference (open points). The different models are: Peters cycle (PC), a
photodisintegration-limited spectrum (PD), a synchrotron-limited diffusion
accelerated spectrum (SDA), a synchrotron-limited one-shot accelerated
spectrum (S1A), a curvature radiation-limited one-shot accelerated spectrum
(C1A), and a universal energy loss spectrum (UEL).

##  IV Signatures of alternative scenarios

Given the difficulty of distinguishing a Peters cycle from alternative
scenarios using high-level information like the UHECR spectrum and composition
it is worthwhile to explore other signatures which might provide a smoking gun
to the dominant process in the universe. In principle, accurate measurements
of the spectra of each mass group in the UHECR spectrum would directly access
information about $\alpha$ and $\beta$, but this level of mass separation is
difficult above the ankle since it would require separation between nuclei of
a single mass group. If data below the ankle is also included, such a mass
separation is possible [28] but it is not possible to distinguish whether the
resulting maximum energy scalings are due to an alternative scenario or a
superposition of source populations. However, other signals may provide
different means to distinguish between a Peters cycle and alternative
scenarios.

First, let us consider the case where there is a substantial proton flux
escaping the source. For small values of $\alpha$, this component will peak at
low energies relative to nuclei – a factor of $2^{-\alpha}$ below the energy-
per-nucleon of nuclei for the Peters cycle family of scenarios. At such low
energies, this proton component likely will be deep into the spectrum’s
Galactic-to-extragalactic transition. While this lower energy than expected
proton component compared to a Peters cycle would be an indication that an
alternative scenario is at work, such a signal may be difficult to distinguish
from both the protons from photodisintegrated nuclei or a second source
population below the ankle.

For large values of $\alpha$, the proton component will peak at higher than
expected energies – a factor of $2^{\alpha}$ above the energy-per-nucleon of
nuclei for the Peters cycle family of scenarios. This would imply that a
significant proton component exists in the spectrum at peaking above heavier
nuclei, which can result in an increase in $\sigma(X_{\mathrm{max}})$ at high
energies (see Fig. 3). Additionally, for large enough values of $\alpha$, this
proton component will extend beyond the GZK threshold and, therefore, produce
a flux of cosmogenic neutrinos at ${\sim}$ EeV energies. Good measurements of
the proton fraction throughout the spectrum or observation of cosmogenic
neutrinos will probe whether such a component exists and measurement of its
peak energy will constrain the value of $\alpha$.

Second, let us consider the case where there is no substantial proton flux
escaping the source. In this case, a proton flux will still arrive at Earth
due to UHECR photodisintegration interactions with the CMB and extragalactic
background light (EBL). These interations preserve the energy-per-nucleon of
the primary CR, so that:

$\displaystyle E_{A\mathrm{PD}}^{p}=\frac{E^{A}}{A}~{}.$ (7)

This implies that the peak in the spectrum of photodisintegrated protons from
primary CRs of mass $A$ will be given by:

$\displaystyle
E_{A\mathrm{PD,max}}^{p}=2^{-\alpha}E_{0}A^{\alpha+\beta-1}~{}.$ (8)

In the case of the Peters cycle family of scenarios, where $\alpha+\beta=1$,
this implies that the spectra of photodisintegrated proton peak at an energy
$2^{-\alpha}E_{0}$ irrespective of their parent CR’s mass $A$. Moreover, the
peak energy of any primary proton spectrum will be the same as the
photodisintegrated spectrum’s up to a factor of $2^{-\alpha}$. Distinguishing
a Peters cycle from alternative scenarios in its family is very difficult and
can only be probed by through this factor of $2^{-\alpha}$ difference in the
peak energies of the primary proton component and photodisintegrated proton
component.

However, this is not true for alternative families of scenarios where
$\alpha+\beta\neq 1$, and the ratio between the peaks in the spectra of
photodisintegrated protons from CRs of mass $A$ and mass $A^{\prime}$ will be
given by:

$\frac{E^{p}_{A\mathrm{PD,max}}}{E^{p}_{A^{\prime}\mathrm{PD,max}}}=\left(\frac{A}{A^{\prime}}\right)^{\alpha+\beta-1}~{}.$
(9)

Therefore, a signature of alternative scenarios outside the Peters cycle
family is a multi-peaked proton component in the UHECR spectrum. While it may
be difficult to resolve these separate peaks, this signature also implies an
extended proton spectrum throughout the entire UHECR spectrum, which is
unexpected for scenarios in the Peters cycle family. This extended proton
component may be more easily detected than resolve the peak energy of spectrum
of each mass group, and can result in larger-than-expected values of
$\sigma(X_{\mathrm{max}})$ at high energies.

An important caveat, however, is that an alternative scenario to a Peters
cycle may be difficult to distinguish from a second UHECR source population
which produces a substantial high-energy proton flux (this possibility has
been explored in a number of studies, including [29, 12, 30, 31]). Such a
source population could mimic alternative scenarios by producing cosmogenic
neutrinos or a larger than expected proton fraction through the spectrum. If
the fluxes of heavier nuclei produced by this source population are not large
enough to be detected, it may be impossible to distinguish between the two.

##  V Conclusions

In this study we explored the degree to which observations of the ultrahigh
energy cosmic ray spectrum and composition above the ankle favor a Peters
cycle over alternative scenarios for spectra escaping sources. Such
alternatives are motivated by energy loss and beyond the Standard Model
processes which imprint particular scalings between the maximum energy of
nuclei and their mass and charge. We find that alternative scenarios explain
the UHECR data above the ankle better than a Peters cycle, regardless of the
hadronic interaction model or systematic data shifts assumed. This result
raises the possibility that a Peters cycle is not realized in any energy range
of the observed UHECR spectrum.

We investigated the observational signatures which might be used to further
discriminate alternative scenarios from a Peters cycle. These include
unexpected scalings of the peak energy of different mass groups at Earth, a
substantial GZK neutrino flux, an extended proton component across the
spectrum’s full energy range, or an unexpected proton component at the highest
energies (which could result in a large $\sigma(X_{\mathrm{max}})$ and
therefore may be constrained by deep neural network-based $X_{\mathrm{max}}$
measurements [32]). However, some of these signatures are difficult to
distinguish from an additional UHECR source population below the ankle or one
producing a substantial high-energy proton flux above it.

We emphasize that alternative scenarios to a Peters cycle represent an
exciting observational opportunity. In particular, further constraints on
these scenarios will not only significantly reduce the theoretical
uncertainties in UHECR modeling and open a window to the conditions inside
UHECR sources, but they will also enable UHECR data to directly constrain new
physics processes.

## Acknowledgments

The work of L.A.A. is supported by the U.S. National Science Foundation (NSF
Grant PHY-2112527). The research of M.S.M. is supported by the NSF MPS-Ascend
Postdoctoral Award #2138121.

## References

  * Peters [1961] B. Peters, Il Nuovo Cimento 22, 800 (1961).
  * Allard  _et al._ [2007] D. Allard, E. Parizot, and A. V. Olinto, Astropart. Phys. 27, 61 (2007), arXiv:astro-ph/0512345 .
  * Taylor [2014] A. M. Taylor, Astropart. Phys. 54, 48 (2014), arXiv:1401.0199 [astro-ph.HE] .
  * Aab  _et al._ [2017] A. Aab _et al._ (Pierre Auger), JCAP 04, 038, [Erratum: JCAP 03, E02 (2018)], arXiv:1612.07155 [astro-ph.HE] .
  * Hillas [1984] A. M. Hillas, Ann. Rev. Astron. Astrophys. 22, 425 (1984).
  * Luce  _et al._ [2022] Q. Luce, S. Marafico, J. Biteau, A. Condorelli, and O. Deligny, Astrophys. J. 936, 62 (2022), arXiv:2207.08092 [astro-ph.HE] .
  * Halim  _et al._ [2023] A. A. Halim _et al._ (Pierre Auger), JCAP 05, 024, arXiv:2211.02857 [astro-ph.HE] .
  * Ehlert  _et al._ [2023a] D. Ehlert, F. Oikonomou, and M. Unger, Phys. Rev. D 107, 103045 (2023a), arXiv:2207.10691 [astro-ph.HE] .
  * Unger  _et al._ [2015] M. Unger, G. R. Farrar, and L. A. Anchordoqui, Phys. Rev. D 92, 123001 (2015), arXiv:1505.02153 [astro-ph.HE] .
  * Globus  _et al._ [2015] N. Globus, D. Allard, and E. Parizot, Phys. Rev. D 92, 021302 (2015), arXiv:1505.01377 [astro-ph.HE] .
  * Kachelrieß  _et al._ [2017] M. Kachelrieß, O. Kalashev, S. Ostapchenko, and D. V. Semikoz, Phys. Rev. D 96, 083006 (2017), arXiv:1704.06893 [astro-ph.HE] .
  * Muzio  _et al._ [2019] M. S. Muzio, M. Unger, and G. R. Farrar, Phys. Rev. D 100, 103008 (2019), arXiv:1906.06233 [astro-ph.HE] .
  * Biehl  _et al._ [2018] D. Biehl, D. Boncioli, A. Fedynitch, and W. Winter, Astron. Astrophys. 611, A101 (2018), arXiv:1705.08909 [astro-ph.HE] .
  * Anchordoqui [2019] L. A. Anchordoqui, Phys. Rept. 801, 1 (2019), arXiv:1807.09645 [astro-ph.HE] .
  * Medvedev [2003] M. V. Medvedev, Phys. Rev. E 67, 045401 (2003), arXiv:astro-ph/0303271 .
  * Ptitsyna and Troitsky [2010] K. V. Ptitsyna and S. V. Troitsky, Phys. Usp. 53, 691 (2010), arXiv:0808.0367 [astro-ph] .
  * Rieger [2011] F. M. Rieger, Int. J. Mod. Phys. D 20, 1547 (2011), arXiv:1107.2119 [astro-ph.CO] .
  * Muzio  _et al._ [2022] M. S. Muzio, G. R. Farrar, and M. Unger, Phys. Rev. D 105, 023022 (2022), arXiv:2108.05512 [astro-ph.HE] .
  * Montero  _et al._ [2023] M. Montero, C. Vafa, and I. Valenzuela, JHEP 02, 022, arXiv:2205.12293 [hep-th] .
  * Anchordoqui [2022] L. A. Anchordoqui, Phys. Rev. D 106, 116022 (2022), arXiv:2205.13931 [hep-ph] .
  * Noble  _et al._ [2023] N. T. Noble, J. F. Soriano, and L. A. Anchordoqui, Phys. Dark Univ. 42, 101278 (2023), arXiv:2306.03666 [hep-ph] .
  * Alves Batista  _et al._ [2016] R. Alves Batista, A. Dundovic, M. Erdmann, K.-H. Kampert, D. Kuempel, G. Müller, G. Sigl, A. van Vliet, D. Walz, and T. Winchen, JCAP 05, 038, arXiv:1603.07142 [astro-ph.IM] .
  * Abreu  _et al._ [2021] P. Abreu _et al._ (Pierre Auger), Eur. Phys. J. C 81, 966 (2021), arXiv:2109.13400 [astro-ph.HE] .
  * Yushkov [2020] A. Yushkov (Auger), PoS ICRC2019, 482 (2020).
  * Riehn  _et al._ [2020] F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev, Phys. Rev. D 102, 063002 (2020), arXiv:1912.03300 [hep-ph] .
  * Pierog  _et al._ [2015] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, Phys. Rev. C 92, 034906 (2015), arXiv:1306.0121 [hep-ph] .
  * Rosenfeld [1975] A. H. Rosenfeld, Ann. Rev. Nucl. Part. Sci. 25, 555 (1975).
  * Abdul Halim  _et al._ [2023a] A. Abdul Halim _et al._ (Pierre Auger), PoS ICRC2023, 438 (2023a).
  * van Vliet  _et al._ [2019] A. van Vliet, R. Alves Batista, and J. R. Hörandel, Phys. Rev. D 100, 021302 (2019), arXiv:1901.01899 [astro-ph.HE] .
  * Muzio  _et al._ [2023] M. S. Muzio, M. Unger, and S. Wissel, Phys. Rev. D 107, 103030 (2023), arXiv:2303.04170 [astro-ph.HE] .
  * Ehlert  _et al._ [2023b] D. Ehlert, A. van Vliet, F. Oikonomou, and W. Winter, (2023b), arXiv:2304.07321 [astro-ph.HE] .
  * Abdul Halim  _et al._ [2023b] A. Abdul Halim _et al._ (Pierre Auger), PoS ICRC2023, 278 (2023b).

##  Appendix A Results for systematic shifts of data

Figures 4 and 5 show the effect of systematically shifting the spectrum and
composition data under both hadronic interaction models considered. We
consider all combinations shifting the energy scale by
$\mathrm{dlg}E=-0.1/0/+0.1$ and shifting the $\langle X_{\mathrm{max}}\rangle$
by $-1/0/+1\sigma_{X}$. As can be seen, these shifts have a non-trivial impact
on the particular regions of parameter space which best describe the data.
However, generically it is the case that no matter the systematic shifts on
consider alternative scenarios describe the data as well or better than a
Peters cycle.

If one has a strong theoretical prior towards a Peters cycle or a particular
alternative scenario, this will strongly effect the systematic data shifts one
favors. Importantly, we note that the overall quality of fit (i.e. the value
of the $\chi^{2}$) depends strongly on the systematic shifts chose regardless
of the $(\alpha,\beta)$ value one considers.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

(i)

Figure 4: Change in quality of fit to the UHECR spectrum and composition
relative to a Peters cycle when considering different systematic shifts of the
data. We consider Sibyll2.3d shifting the data by $\mathrm{dlg}E=-0.1/0/+0.1$
(left/middle/right columns) and $-1/0/+1\sigma_{X}$ (top/middle/bottom
rows).The Peters cycle (PC) and a number of alternative scenarios (green dots)
are highlighted: a photodisintegration-limited spectrum (PD), a synchrotron-
limited diffusion accelerated spectrum (SDA), a synchrotron-limited one-shot
accelerated spectrum (S1A), a curvature radiation-limited one-shot accelerated
spectrum (C1A), and a universal energy loss spectrum (UEL).

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

(i)

Figure 5: Same as Fig. 4 for Epos-LHC.

##  Appendix B Alternative scenario escaping spectra

Figures 6 and 7 show the best-fit escaping spectra for alternative scenarios
compared to a Peters cycle for our benchmark systematic data shifts. These
plots underscore the similarity of scenarios belonging to the same
$\alpha+\beta$ family: Peters cycle & photodisintegration ($\alpha+\beta=1$),
and synchrotron-limited diffusion acceleration & universal energy losses
($\alpha+\beta=0$).

The best-fit scenario ($\alpha\simeq 6.75$, $\beta\simeq-5$) also demonstrates
how scenarios with a significant proton component and large $\alpha$ result in
a proton component which does not obey the normal mass ordering expected in a
Peters cycle. In this case, the proton component peaks at similar energies as
the CNO component, resulting in a peak in the $\sigma(X_{\mathrm{max}})$ at
high energies (see Fig. 3).

(a)  (b)

(c)  (d)

(e)  (f)

Figure 6: Best-fit escaping spectra for non-Peters cycle scenarios (solid
lines) compared to those for a Peters cycle (dashed lines). Spectra are broken
down by mass group (colored lines). All plots are for Sibyll2.3d with data
shifted by $\mathrm{dlg}E=+0.1$ and $-1\sigma_{X}$.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 7: Same as Fig. 6 but for Epos-LHC.


[a]Simone Bacchio

# Lyncs-API: a Python API for Lattice QCD applications

Jacob Finkenrath  Christodoulos Stylianou

###### Abstract

We present Lyncs-API, a Python API for Lattice QCD applications currently
under development. Lyncs aims to bring several widely used libraries for
Lattice QCD under a common framework. Lyncs flexibly links to libraries for
CPUs and GPUs in a way that can accommodate additional computing architectures
as these arise, achieving performance-portability for the calculations while
maintaining the same high-level workflow. Lyncs distributes calculations using
Dask and mpi4py, with bindings to the libraries automatically generated by
cppyy. While Lyncs is designed to allow linking to multiple libraries, we
focus on a set of targeted packages that include DDalphaAMG, tmLQCD, QUDA and
c-lime. More libraries will be added in the future. We also develop generic-
purpose tools for facilitating the usage of Python in Lattice QCD and HPC in
general. The project is open-source, community-oriented and available on
Github.

##  1 Introduction

Three keywords in High-Performance Computing (HPC) are Performance,
Portability and Productivity, known as the three “P”. Performance is commonly
the main focus of Lattice QCD applications, which achieve scalability and
high-performance on major computing architectures. Portability and
productivity, instead, are new emerging requirements, usually sacrificed at
the cost of performance and fast development, but recently widely chased in
response to the new challenges arising in HPC. In this work, we introduce the
Lyncs-API, a new community effort whose goal is to achieve the three “P” by
developing a Python framework for Lattice QCD calculations. In particular, we
aim to portability and productivity by designing a framework that,
respectively, (i) can achieve high-performance on the largest HPC systems and
(ii) is flexible, user-friendly, with a clean API and can be easily integrated
with other tools as well.

###  1.1 Related works in Lattice QCD

Recently, significant efforts towards portability are being made by developers
of various libraries for Lattice QCD. For example, the QUDA library [1] is
implementing various backends for extending its employment on GPUs made by
other manufacturers than NVIDIA, such as AMD or Intel [2]. And the Grid
library, whose development focused on achieving portability on CPU
architectures [3], has been also implementing GPU support [4]. There are a
handful of more examples where the lattice QCD community is either extending
or developing new software for adapting to the diversity of computing
architectures available nowadays. The Lyncs-API is in its way different from
these efforts by not implementing directly Lattice QCD kernels but relying on
the availability of libraries for providing high-coverage on tools,
portability and performance. Indeed, one of its main goals is to bring several
widely used libraries for Lattice QCD under a common framework and involve a
large community.

On the other hand, productivity is probably the least satisfied feature by
Lattice QCD software and it is also the most challenging to achieve. Indeed,
most libraries often lack of user-friendliness, flexibility and backwards-
compatibility making them difficult to be integrated in applications and also
interact between each other. The Lyncs-API, instead, has productivity as its
main focus. It is the reason behind the choice of Python as high-level
programming language and most of its features are inspired by the impressive
achievements of the Python community in this direction. A compatible effort is
being made by the Grid Python Toolkit (GPT) [5] which provides a Python
interface and Python bindings for Grid. We foresee the usage of GPT in the
future for interfacing to the Grid library.

###  1.2 Overview of related Python packages

In recent years, Python has become the most used programming language in data
science. This is thanks to its strong community and many successful
applications, e.g. in machine learning. The employment of Python in HPC,
instead, is limited but growing and there is a number of tools meant to
support it. In the following we will provide an brief overview of key packages
used by the Lyncs-API, which facilitate the usage of Python in HPC.

Numpy

[6] is the most well-known Python package. It offers performant
implementations of multi-dimensional array operations and it defines a
standard employed by many other packages.

CuPy

[7] is an implementation of Numpy-compatible multi-dimensional array
operations on GPUs using CUDA. It supports multiple devices and Just-In-Time
(JIT) compilation of kernels. Experimentally, it also supports AMD GPUs using
ROCm.

cppyy

[8] is an automatic, run-time, Python-C++ bindings generator, for calling
C/C++ from Python and Python from C++. Libraries are required to be compiled
as shared objects and bindings are generated automatically by parsing the
header files. The run-time generation enables detailed specialization for
higher performance, lazy loading, Python-side cross-inheritance and callbacks
for working with C++ frameworks (e.g. template instantiation).

mpi4py

[9] (MPI for Python) provides Python bindings for the Message Passing
Interface (MPI) standard with an object oriented interface resembling MPI-2
C++. Since version 3.1.0, it supports GPUDirect RDMA (GDR) communications with
GPU-aware MPI implementations.

h5py

[10] implements a Python interface for the HDF5 library and file format. It
efficiently supports all features of the HDF5 standard, including parallel IO.

Dask

[11] provides a Python API for data- and task-parallelism. It supports
distributed Numpy arrays via domain decomposition and uses a centralized
scheduler for managing multiple workers. Another focus of the package is
remote data representation via futures.

Numba

[12] is an open source JIT compiler that translates a subset of Python and
NumPy code into performant machine code. It supports various CPU and GPU
architectures and, respectively, Numpy and CuPy arrays. It also supports
caching or Ahead-of-Time (AOT) compilation for avoiding repeated overheads
during run-time.

##  2 Overview of the Lyncs-API

In the following we provide an overview of the Lyncs-API, focusing on its
scope, design and structure. The package as a whole is currently under
development and some features described below might be missing or incomplete.
Others instead have reached a stable state. In this section we will not
discuss implementation status but we will provide a time-agnostic description
of the project valid in the future. Implementation details are given in the
following section.

###  2.1 Statement of purpose

The Lyncs-API aims to provide a complete framework for performing Lattice QCD
calculations using the Python programming language. The project is community-
oriented, open source and welcomes contributions. Its motivations and purpose
are the following:

  * •

To create a productive, user-friendly and ease-to-use framework for Lattice
QCD calculations. Python is chosen as top-level programming language for
implementing the API and it is used for driving and deploying calculations.
All computational kernels, instead, are executed by compiled software and
libraries developed in C/C++.

  * •

To achieve performance and portability on HPC computing architectures by
interfacing to major libraries for Lattice QCD. Many libraries are targeted in
order to reach high-coverage on tools and performance-portability. Moreover
many implementations of the same kernels allow for crosschecks and performance
benchmarks. Interfaces to open-source and free-to-use libraries are welcome
and the project aims to involve the whole community.

  * •

To target multitasking, modular computing and resiliency. The high-level
framework is designed to drive calculation on vast partitions of HPC systems
beyond scalability of single kernels. This is done via multitasking as a mean
to seamlessly increase scalability of applications, to deploy computations
simultaneously on partitions with different architectures, i.e. modular
computing, and to allow for calculations to proceed in the event of node
failure having more nodes at disposal and data duplicates, i.e. resiliency.

  * •

To follow high-standards in community software development. We employ version
control, unit-testing targeting high-coverage, backwards compatibility,
linting, code formatting, issues tracking and updated documentation. We
publish on GitHub and use GitHub Actions for Continuous Integration and
Deployment (CI/CD). Packages are distributed on the Python Packaging Index
(PyPI) and can be installed via pip. Code coverage is tracked in the
implementation and new tests are implemented for new functionalities.

###  2.2 Design and structure

The Lyncs-API is designed upon a modular structure in order to separate
concerns of different modules and maintain their implementation and deployment
independent from each other as much as possible. All modules together aim to
provide a complete Python framework for Lattice QCD calculations but, on the
other hand, they can be used separately within other Python projects allowing
for easy code reuse. The modules are developed and deployed independently. We
use version control for tracking changes in the software stack and update
setup requirements appropriately to ensure that packages point to the correct
versions of their dependencies. All packages are collected on the GitHub
organization111https://github.com/Lyncs-API with description and documentation
in the respective README files and overview on the online
webpage222https://lyncs-api.github.io.

The current structure of packages in the Lyncs-API can be divided in three
categories: (i) generic purpose tools, (ii) interfaces to Lattice QCD
libraries, and (iii) high-level tools for Lattice QCD calculations. In the
first category we find a variety of generic purposes packages that could be
also used outside Lattice QCD applications. In the second category we have a
number of interfaces to Lattice QCD libraries providing a “Pythonization” of
their functions/classes resembling as close as possible the library itself.
And in the last category we find high-level modules for lattice QCD
calculations that provide end-user functionalities based on the low-level
packages in the previous categories. More details on the package structure are
given in the following sub-sections.

####  2.2.1 Generic purpose tools

Part of the Lyncs-API are a number of packages that provide generic purpose
tools whose development has been needed in constructing the API itself.
Briefly, we have (i) lyncs_setuptools, which provides setup-level
functionalities and developer tools including support for CMake; (ii)
lyncs_utils, which is a collection of small and independent functionalities
that recur widely in the implementation of all packages; (iii) lyncs_cppyy,
which provides a wrapping layer on top of cppyy for simplifying its usage when
interfacing to libraries (more details in Section 3.1.1); (iv) lyncs_mpi with
tools for MPI using mpi4py and in combination with Dask ; and lastly we have
two major packages, lyncs_io and lyncs_field, described in the following.

lyncs_io aims to provide an all-inclusive framework for IO in Python. It
responds to the necessity of missing tools for parallel IO and user-
friendliness in various file formats. Therefore it implements a very simple
interface with two main functionalities, load and save; it supports various
file formats and, where possible, parallel IO for domain decomposed arrays.
Additionally, it also implements file formats that are common in Lattice QCD
as, for instance, the Lattice QCD Interchange Message Encapsulation (LIME).
More details are given in Section 3.3.

lyncs_field333currently under development provides a set of meta-classes for
abstraction of properties of multi-dimensional arrays and tensor fields.
Indeed, an issue when interfacing to various libraries is the variety of
formats and conventions used for storing the data. Especially in Lattice QCD,
many choices can be made in ordering the axis, whether using a checkerboard
splitting, how to store the halos, etc. With lyncs_field one can easily track
information on the properties of arrays and seamlessly convert from a standard
to another. Numpy and CuPy are the default backends for array’s operations on
CPU and GPU architectures, respectively, but also any Numpy-compatible backend
is supported.

####  2.2.2 Interfaces to Lattice QCD libraries

A major task of the Lyncs-API is to provide Python interfaces to a wide range
of libraries. Within the scope of the Lyncs-API is to target both, cutting-
edge software packages as well as legacy software, in order to achieve on one
side high-coverage on tools and performance-portability as well on the other
to perform crosschecks and benchmarks ensuring continuity in the developments.
Currently the number of interfaces is limited to lyncs_DDalphaAMG, lyncs_quda,
lyncs_tmLQCD and lyncs_clime that interface respectively to DDalphaAMG [13],
QUDA [14], tmLQCD [15] and c-lime444https://usqcd-software.github.io/c-lime/.
The interfaces follow as closely as possible the structure of the libraries in
order to simplify their implementation and to remain familiar to users of the
library itself. Major differences may appear when libraries are not object-
oriented. We expect in the future to increase significantly the list of
targeted libraries while enlarging the community. More details on the
interfaces are given in Section 3.1 with examples in Section 3.1.2.

####  2.2.3 High-level tools for Lattice QCD

Apart from interfaces to libraries, all other tools for Lattice QCD are
expected to be on the top level only. Here we find generic APIs, abstraction
of common workflows and applications. For instance, we have implemented tools
for Hybrid Monte-Carlo (HMC) in lyncs_hmc and we aim to provide a complete
end-user API in lyncs . Unfortunately, this is currently in an early
development stage and we do not have a complete list of packages yet.

##  3 Implementation details

As typical when developing high-level software, the design of the Lyncs-API
proceeded top-down while its implementation is proceeding bottom-up. For this
reason, not all features described in the previous section are available, but,
on the other hand, the modular structure of the Lyncs-API allows for most of
the developed components to be used as stand-alone packages. Hereafter we
provide implementation details on key components with examples.

###  3.1 Interfaces to libraries

All packages that interface to libraries implement the following procedure:

  1. 1.

Library’s source-code is cloned during installation from its online repository
targeting a chosen release (or commit). We fix the latter, because most
libraries do not ensure backward compatibility and/or stability. Update of the
library’s release is done on a regular basis after respectively updating the
interface and testing all functionalities.

  2. 2.

Patches to the source-code are applied where needed and they are distributed
within the interface repository. Patches either fix installation or execution
issues, e.g. enabling compilation of the library as shared-object or exposing
more features in the interface.

  3. 3.

The above items, compilation and installation of the library are managed by
CMake. Indeed, libraries are not compiled AOT and are not distributed as
Python wheels for optimization purposes. The installation of the package,
including the execution of CMake, follows the default Python procedure with
setup script in setup.py that can be installed via pip.

  4. 4.

Python bindings are automatically generated via cppyy using our wrapper
lyncs_cppyy, which simplifies its usage with some additional features (see
3.1.1).

  5. 5.

A uniform and well documented Python interface is created for key features of
the library resembling its structure. Data conversion and input validation are
also part of this step, using typical Python features in data management and
error handling.

  6. 6.

Unit-testing and high-coverage of the interfaced functionalities is reached.
Functions are tested for many lattice sizes, parallel partitioning, data
formats, floating precision, etc. sometimes finding issues and providing
feedbacks for the developers of the library itself.

####  3.1.1 Automatic bindings with cppyy

Giving more details on the fourth step, cppyy automatically generates Python
bindings for C/C++ libraries by parsing header files and loading the compiled
library at run-time as a shared object. By collecting information on the
signature of functions and classes in the headers, it can then allocate
appropriately memory, manage data, call functions, convert Python objects to
inputs and vice versa for outputs. This process is based on Cling [16], an
interactive C++ interpreter, built on the top of LLVM and Clang libraries with
specific Python functionalities provided by cppyy. The package also supports
JIT compilation of C++ templates and headers-only implementations. The
standard API of cppyy provides two functions for this purpose include and
load_library. In lyncs_cppyy we use an object-oriented approach representing
the library as a class. For example,

⬇

from lyncs_cppyy import Lib

lib = Lib(header="quda.h", library="quda", namespace="quda")

the above can be used to load e.g. quda and access its functionalities as
attributes/methods of the lib object. Additional features compared to cppyy is
the usage of a default namespace, lazy instantiation (loading the library at
its first use), definition and access to macros. Additional remarks on using
cppyy compared other standard approaches are

  * •

Resiliency to changes: since interfaces are not compiled AOT, all
functionalities not affected by changes are operative and usable. Any issue
will appear at run-time. This is on one side an advantage when transiting from
a version of a library to another, but on the other interfaces require full
coverage for catching any broken features.

  * •

Error tolerance: in connection to the above, cppyy is error tolerant and
allows for the usage of “try/except” when calling functions or instantiating
classes. This is useful for supporting multiple versions of libraries and
easily addressing changes. Moreover, cppyy can also convert abortive signals,
e.g. segmentation faults, into Python exceptions avoiding to interrupt
execution and easily debug software.

  * •

Class cross-inheritance: another very useful feature is cross-inheritance of
C++ classes in Python allowing the latter to be called by C++. Indeed, method
overloading and virtual method overriding are fully supported.

####  3.1.2 Examples

Hereafter we provide two examples for computing the inverse of the Dirac
operator using Lyncs’ interfaces.

⬇

import lyncs_quda as quda

from mpi4py import MPI

# Creating a Cart. comm via MPI

lattice = [4, 4, 4, 4]

procs = [2, 2, 1, 1]

comm = MPI.COMM_WORLD

comm = comm.Create_cart(procs)

# Creating a random gauge field

gauge = quda.gauge(lattice=lattice,

comm=comm, device="GPU")

gauge.random()

plaq = gauge.plaquette

print("Plaquette:", plaq)

# Creating a random spinor

rhs = quda.spinor_like(gauge)

rhs.random()

# Inverting using CG

dirac = gauge.Dirac(kappa=0.125)

mat = dirac.MMdag

sol = mat.solve(rhs, inv_type="CG")

res = mat(sol) - rhs

⬇

from lyncs_DDalphaAMG import Solver

from lyncs_mpi import Client

# Creating a Cart. comm via Dask

lattice = [4, 4, 4, 4]

procs = [2, 2, 1, 1]

client = Client()

comm = client.create_cart(procs)

# Instantiating DDalphaAMG

solver = Solver(lattice=lattice,

comm=comm, kappa=0.125)

# Reading a random configurations

conf = solver.read_configuration(

"test/conf.random")

plaq = solver.set_configuration(

conf)

print("Plaquette:", plaq)

# Inverting using multigrid

rhs = solver.random()

sol = solver.solve(rhs)

res = solver.D(sol) - rhs

res = res.compute()

The first example (left panel) uses QUDA via lyncs_quda and the second uses
DDalphaAMG via lyncs_DDalphaAMG. Notation in the two implementation differ
slightly since the interfaces follow closesly the structure of the libraries.
The interface to QUDA is based on its C++ implementation (not the C interface)
and the interface to DDalphaAMG is based on the header
DDalphaAMG.h555https://github.com/sbacchio/DDalphaAMG/blob/master/src/DDalphaAMG.h.
Additionally the first example parallelizes directly via MPI while the second
uses Dask at the higher level which then uses MPI. The two approaches are
interchangeable and they have a very similar initialization (i.e. need to
exchange how the communicator is created in the two examples). On the other
hand significant differences in the execution apply as outlined in the
following section.

###  3.2 Parallel computing with Dask and MPI

Beside standard data-parallelism via domain decomposition using MPI, the
Lyncs-API also employs Dask for more advanced features, e.g. task parallelism.
In Lyncs the usage of Dask is optional, made compatible with MPI and ease to
use (see above example). The core of Dask is a centralized scheduler that
receives instructions from a client and assigns tasks to workers. Within the
duties of the scheduler are knowledge on the data location, properties and
status of workers, management and transfer of data, submission of tasks and
dependencies resolution. On top of that, Dask provides high-level interfaces
like distributed Numpy arrays via domain decomposition which is of interest
for Lattice QCD and it supports all features needed including halos. In
lyncs_mpi we have developed tools for bridging the gap between Dask and MPI,
e.g. grouping workers within a communicator and easy integration in classes.

Going into details of the above examples, while in the left panel the
computation proceeds as expected in a MPI-distributed application, in the
right panel most of the calculations are delayed. Starting from the creation
of the communicator, nothing is actually computed until “.compute()” is
called. This means that all steps are collected in a computational graph,
which is then passed to the scheduler and computed by the first group of
“workers” available.

###  3.3 Simplified parallel IO in Python

Another outcome of the Lyncs-API is a high-level tool for parallel I/O in
Python within lyncs_io. The package has a general purpose and implements a
simple API for loading and saving distributed data using domain decomposition
via either Dask or an MPI Cartesian communicator. The objectives of the
package are

  * •

Seamlessly IO: we have tried to simplify the IO as much as possible, without
having to open files, directly specify the format or using more than one line
of code. Three main functions are provided save(obj, filename) for saving,
obj=load(filename) for loading and obj=head(filename) for retrieving metadata
without loading data. More arguments are available as outlined in the
following.

  * •

Many formats supported. The file format can be either deduced from the
filename’s extension or with the option format passed as an argument to the
above functions. The structure of the package is flexible enough to easily
accommodate new/customized file formats as these arise. For instance we
support file formats typical in lattice QCD like LIME.

  * •

Omission of extension When saving, if the extension is omitted, the optimal
file format is deduced from the data type and the extension is added to the
filename. When loading, any extension is considered, i.e. filename.*, and if
only one match is available, the file is loaded.

  * •

Support for archives. In case of archives, e.g. HDF5, zip etc., the content
can be accessed directly by specifying it in the path. For instance using
directory/file.h5/content, directory/file.h5 is the file path, and the
remaining is content to be accessed that will be searched inside the file
(this is inspired by h5py ).

  * •

Support for Parallel IO. Where possible and implemented, parallel IO is
supported. This is enabled either via MPI providing a Cartesian communicator
with the option comm, or via Dask with the option chunks.

##  4 Conclusions and future plans

In this contribution we have presented the pillars of the Lyncs-API, a new
software effort for Lattice QCD applications using the Python programming
language. We have described its design, purpose and long term vision. We have
provided details on the implementation of its low-level building blocks, such
as interfaces to libraries, data- and task-parallelism and I/O. We have
motivated the choice of Python as a tool for easily addressing new
computational challenges in HPC and Lattice QCD, such as the integration of
machine learning techniques, usage of task parallelism and achieving the
three-Ps towards Exascale supercomputing. The development of the Lyncs-API is
strongly motivated by these and wants to be an alternative to more
traditional/low-level approaches.

In the short term, our plan is to complete the interface to QUDA, lyncs_quda,
closely following the developments of the library towards 2.0, and to
implement lyncs_field which will be the connecting layer between all libraries
allowing for easy exchange and reordering of data. This will be followed by
the implementation of more high-level tools and first large-scale
applications. A major target of our development is to achieve scalability in
HMC simulations using modular computing and multi-tasking. Concluding, with
this presentation we wanted to raise awareness on the project and welcome
contributions and/or feedback.

## Acknowledgments

This project has received funding under PRACE-6IP, Grant agreement ID: 823767,
Project name: LyNcs. LyNcs is one of 10 applications supported by PRACE-6IP,
WP8 “Forward Looking Software Solutions”. S.B., J.F. and C.S. have received
funding under this project. The authors would also like to thank the other
members of the LyNcs project for the stimulating collaboration.

## References

  * [1] M.A. Clark, R. Babich, K. Barros, R.C. Brower and C. Rebbi, _Solving Lattice QCD systems of equations using mixed precision solvers on GPUs_ , _Comput. Phys. Commun._ 181 (2010) 1517 [0911.3191].
  * [2] K. Clark, _Preparing for quda 2.0_ , in _Lattice2021_ , https://indi.to/SgQDB.
  * [3] P. Boyle, A. Yamaguchi, G. Cossu and A. Portelli, _Grid: A next generation data parallel C++ QCD library_ , 1512.03487.
  * [4] P. Boyle and M. Lin, _Performance portability strategies for grid c++ expression templates_ , .
  * [5] C. Lehner, _Grid python toolkit (gpt)_ , in _Lattice2021_ , https://indi.to/XXwcm.
  * [6] 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.
  * [7] R. Okuta, Y. Unno, D. Nishino, S. Hido and Crissman, _Cupy : A numpy-compatible library for nvidia gpu calculations_ , 2017.
  * [8] W.T. Lavrijsen and A. Dutta, _High-performance python-c++ bindings with pypy and cling_ , in _2016 6th Workshop on Python for High-Performance and Scientific Computing (PyHPC)_ , pp. 27–35, IEEE, 2016.
  * [9] L. Dalcín, R. Paz and M. Storti, _Mpi for python_ , _Journal of Parallel and Distributed Computing_ 65 (2005) 1108.
  * [10] A. Collette, _Python and HDF5: unlocking scientific data_ , " O’Reilly Media, Inc." (2013).
  * [11] Dask Development Team, _Dask: Library for dynamic task scheduling_ , 2016\.
  * [12] S.K. Lam, A. Pitrou and S. Seibert, _Numba: A llvm-based python jit compiler_ , in _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_ , LLVM ’15, (New York, NY, USA), Association for Computing Machinery, 2015, DOI.
  * [13] S. Bacchio, C. Alexandrou, J. Finkenrath, A. Frommer, K. Kahl and M. Rottmann, _DDalphaAMG for Twisted Mass Fermions_ , _PoS_ LATTICE2016 (2016) 259 [1611.01034].
  * [14] K. Clark, M. Wagner, D. Howarth, E. Weinberg, A. Strelchenko, J. Tu et al., _lattice/quda: Quda v1.1.0_ , Oct., 2021.  10.5281/zenodo.5610079.
  * [15] K. Jansen and C. Urbach, _tmLQCD: A Program suite to simulate Wilson Twisted mass Lattice QCD_ , _Comput. Phys. Commun._ 180 (2009) 2717 [0905.3331].
  * [16] V. Vasilev, P. Canal, A. Naumann and P. Russo, _Cling – the new interactive interpreter for ROOT 6_ , _Journal of Physics: Conference Series_ 396 (2012) 052071.